---
title: "From CFA to a defensible measurement model"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{From CFA to a defensible measurement model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
library(nomologR)

knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.width = 7,
  fig.height = 4.8
)
```

# Why measurement evidence is staged

A confirmatory factor model, a reliability coefficient, and a discriminant-validity
statistic answer different questions. Treating them as interchangeable checkboxes
can produce misleading conclusions. `nomologR` therefore keeps the workflow staged:

1. `nomo_cfa()` — **Does the proposed measurement structure reproduce the data reasonably, and where is there model strain?**
2. `nomo_reliability()` — **How consistently does the intended score represent its target under that measurement model?**
3. `nomo_validity()` — **Do indicators converge on their intended constructs, and are theoretically distinct constructs empirically distinguishable?**

The sequence is evidence-guided rather than pass/fail. A reference value can trigger
inspection, but it never automatically deletes an indicator, merges constructs, or
establishes construct validity.

# A healthy two-factor model

We begin with a known two-factor population in which both factors are well measured
and moderately correlated.

```{r healthy-data}
set.seed(2026)
n <- 500

f1 <- rnorm(n)
f2 <- .30 * f1 + sqrt(1 - .30^2) * rnorm(n)

healthy <- data.frame(
  A1 = .82 * f1 + rnorm(n, sd = .55),
  A2 = .79 * f1 + rnorm(n, sd = .58),
  A3 = .76 * f1 + rnorm(n, sd = .62),
  A4 = .80 * f1 + rnorm(n, sd = .57),
  B1 = .83 * f2 + rnorm(n, sd = .54),
  B2 = .78 * f2 + rnorm(n, sd = .60),
  B3 = .75 * f2 + rnorm(n, sd = .63),
  B4 = .81 * f2 + rnorm(n, sd = .56)
)

model <- '
  F1 =~ A1 + A2 + A3 + A4
  F2 =~ B1 + B2 + B3 + B4
'
```

```{r healthy-fit}
cfa <- nomo_cfa(model, data = healthy)

# Modest resampling here keeps vignette build time reasonable.
# Use more resamples for final reporting.
rel <- nomo_reliability(
  cfa,
  ci = "bootstrap",
  ci_boot = 100,
  ci_seed = 2026
)
val <- nomo_validity(cfa, htmt = "both")
```

The compact summaries are intentionally different from the raw engine output. They
surface the evidence most relevant to the user's next decision without repeating the
same full loading table in several places.

```{r healthy-tables}
summary(rel)
summary(val)
```

The figures also avoid duplication. Item-level loading evidence belongs to the CFA;
reliability gets a coefficient plot; convergent validity gets an AVE plot; and
construct separation gets an HTMT-family plot.

```{r healthy-loading-plot}
plot(cfa, type = "loadings")
```

```{r healthy-reliability-plot}
plot(rel)
```

```{r healthy-ave-plot}
plot(val, type = "ave")
```

```{r healthy-discriminant-plot}
plot(val, type = "discriminant")
```

A healthy result contributes several complementary pieces of evidence. It does not
turn the measurement model into a permanently "validated" object. Validity remains
an accumulating argument tied to construct interpretation, population, setting, and
intended use.

# Good global fit does not guarantee good measurement

This is a critical teaching case. The following one-factor model is correctly
specified, but its indicators are weak measures of the latent variable.

```{r weak-data}
set.seed(2027)
n <- 650
f <- rnorm(n)

weak <- data.frame(
  W1 = .30 * f + rnorm(n, sd = .95),
  W2 = .34 * f + rnorm(n, sd = .94),
  W3 = .28 * f + rnorm(n, sd = .96),
  W4 = .32 * f + rnorm(n, sd = .95),
  W5 = .31 * f + rnorm(n, sd = .95)
)

weak_cfa <- nomo_cfa(
  'Weak =~ W1 + W2 + W3 + W4 + W5',
  data = weak
)
weak_rel <- nomo_reliability(weak_cfa)
weak_val <- nomo_validity(weak_cfa, htmt = "none")
```

```{r weak-fit}
weak_cfa$fit_evidence
```

Because the one-factor structure is the correct covariance model, global fit can look
quite good. That does **not** imply that the items strongly represent the construct or
that the resulting score is precise.

```{r weak-measurement}
summary(weak_rel)
summary(weak_val)
```

This is why `nomologR` does not collapse measurement quality into a single global-fit
verdict. Weak standardized loadings, low AVE, and weak model-based reliability remain
important even when CFI, RMSEA, or SRMR look favorable.

# Strong reliability does not guarantee construct separation

The reverse problem is also possible. Two constructs can each be measured very
consistently while remaining difficult to distinguish empirically.

```{r overlap-data}
set.seed(2028)
n <- 700
f1 <- rnorm(n)
f2 <- .94 * f1 + sqrt(1 - .94^2) * rnorm(n)

overlap <- data.frame(
  A1 = .86 * f1 + rnorm(n, sd = .48),
  A2 = .83 * f1 + rnorm(n, sd = .51),
  A3 = .84 * f1 + rnorm(n, sd = .50),
  B1 = .86 * f2 + rnorm(n, sd = .48),
  B2 = .83 * f2 + rnorm(n, sd = .51),
  B3 = .84 * f2 + rnorm(n, sd = .50)
)

overlap_model <- '
  F1 =~ A1 + A2 + A3
  F2 =~ B1 + B2 + B3
'

overlap_cfa <- nomo_cfa(overlap_model, data = overlap)
overlap_rel <- nomo_reliability(overlap_cfa)
overlap_val <- nomo_validity(overlap_cfa, htmt = "both")
```

```{r overlap-summary}
summary(overlap_rel)
summary(overlap_val)
```

Here, strong omega and AVE do not erase a high latent correlation or HTMT2 value.
The package therefore recommends investigating theoretical distinctiveness, item
content, cross-construct overlap, and the intended construct boundary. It does not
automatically merge factors or delete indicators.

```{r overlap-plot}
plot(overlap_val, type = "discriminant")
```

# Ordered indicators: keep the score scale explicit

For ordered Likert-type indicators, `nomo_cfa()` uses an ordered-data CFA workflow
when the indicators are declared ordered. `nomo_reliability()` then distinguishes the
observed ordinal-score reliability estimand from the hypothetical latent-response
scale.

```{r ordinal-example}
set.seed(2029)
n <- 500
f <- rnorm(n)
z <- data.frame(
  I1 = .82 * f + rnorm(n, sd = .60),
  I2 = .79 * f + rnorm(n, sd = .62),
  I3 = .76 * f + rnorm(n, sd = .65),
  I4 = .80 * f + rnorm(n, sd = .60)
)

ordinal <- as.data.frame(lapply(z, function(x) {
  ordered(cut(x, breaks = c(-Inf, -.8, -.25, .25, .8, Inf), labels = 1:5))
}))

ordinal_cfa <- nomo_cfa(
  'F =~ I1 + I2 + I3 + I4',
  data = ordinal,
  ordered = names(ordinal)
)
ordinal_rel <- nomo_reliability(
  ordinal_cfa,
  ordinal_scale = TRUE,
  include_alpha = TRUE
)
summary(ordinal_rel)
```

Observed-scale omega is available for the practical ordinal composite. Observed-scale
alpha is deliberately reported as unavailable in this ordered-CFA estimand rather than
silently replacing it with a different definition. The point is not to maximize the
number of coefficients reported; it is to keep the score and estimand clear.

# Continuing the exploratory example

The [exploratory walkthrough](exploratory-workflow.html) flagged a cross-loading
item (`a5`) and a weak item (`b5`) in `nomo_demo_continuous`. Suppose the
researcher prespecifies the intended two-factor simple structure with all ten
items and fits it as a CFA:

```{r demo-cfa}
demo_model <- nomo_model(list(
  A = c("a1", "a2", "a3", "a4", "a5"),
  B = c("b1", "b2", "b3", "b4", "b5")
))

demo_cfa <- nomo_cfa(demo_model, data = nomo_demo_continuous)
demo_cfa$fit_evidence
head(demo_cfa$top_modification_indices, 5)
```

Three lessons carry over from the exploratory stage:

- **Modification indices are quarantined.** The largest index points toward the
  cross-loading that the data-generating model contains, but `nomologR` does not
  add it. Allowing a cross-loading, revising the item, or evaluating a model
  without it is a theoretical decision; freeing parameters because an index is
  large is capitalization on chance (MacCallum, Roznowski, & Necowitz, 1992).
- **The same sample is not independent confirmation.** This CFA uses the data
  that suggested the structure. `nomo_split()` or a new sample provides a
  stronger test.
- **Missing data are handled explicitly.** By default lavaan uses complete cases
  here; `missing = "fiml"` requests full-information maximum likelihood for
  continuous indicators, and the decision is recorded. Whether the results
  depend on that choice is a separate question, taken up
  [below](#does-a-result-depend-on-how-missing-data-were-handled).

# Comparing models with a recorded rationale

`nomo_compare()` evaluates competing measurement models side by side. It
requires a rationale, checks whether the models are nested, reports the
difference test that matches the estimator when they are, and never selects a
model automatically.

## A cross-loading suggested by the results

The largest modification index above pointed to a loading of `a5` on factor B.
Because that idea came from the results rather than from theory stated in
advance, the comparison is labeled post hoc:

```{r compare-cross}
cross_cfa <- nomo_cfa(
  "A =~ a1 + a2 + a3 + a4 + a5\nB =~ b1 + b2 + b3 + b4 + b5 + a5",
  data = nomo_demo_continuous
)

cross_comparison <- nomo_compare(
  simple_structure = demo_cfa,
  cross_loading = cross_cfa,
  rationale = "The largest modification index suggested that a5 also reflects factor B.",
  origin = "post_hoc"
)
cross_comparison
```

```{r compare-values, include=FALSE}
fmt_p <- function(p) {
  if (!is.finite(p)) return("unavailable")
  if (p < .001) "< .001" else paste("=", formatC(p, format = "f", digits = 3))
}
cc <- cross_comparison$comparisons
```

Estimating the cross-loading reduces misfit (chi-square difference =
`r sprintf("%.1f", cc$chisq_diff)`, df = `r cc$df_diff`, p `r fmt_p(cc$p_value)`;
CFI changes by `r sprintf("%+.3f", cc$delta_cfi)`). That is exactly the feature
built into these simulated data, but with real data the evidence alone does not
settle the question. Does the wording of `a5` plausibly reflect both
constructs? Because the comparison is post hoc, the decision log records it as
such and recommends confirming the retained model in independent data.

## Is a weak item needed?

`b5` loads weakly on B. Dropping `b5` from the model changes the data being
modeled, so a model without the column cannot be tested against the full model.
Instead, keep `b5` and fix its loading to zero:

```{r compare-b5}
b5_zero_cfa <- nomo_cfa(
  "A =~ a1 + a2 + a3 + a4 + a5\nB =~ b1 + b2 + b3 + b4 + 0*b5",
  data = nomo_demo_continuous
)

b5_comparison <- nomo_compare(
  full = demo_cfa,
  b5_loading_zero = b5_zero_cfa,
  rationale = paste(
    "Evaluate whether the weakly loading item b5 contributes to factor B",
    "before deciding whether to keep it."
  )
)
summary(b5_comparison)
```

```{r b5-values, include=FALSE}
bc <- b5_comparison$comparisons
b5_loading <- b5_comparison$loadings$full[b5_comparison$loadings$item == "b5"]
```

Fixing the loading to zero worsens fit (chi-square difference =
`r sprintf("%.1f", bc$chisq_diff)`, df = `r bc$df_diff`, p `r fmt_p(bc$p_value)`),
so `b5` is statistically related to B. Its standardized loading, however, is
only `r sprintf("%.2f", b5_loading)`. With `r b5_comparison$models$n_used[[1]]`
cases, even a weak relation is detectable; statistical detectability is not the
same as an adequate indicator. The decision still depends on whether `b5` covers
content the construct needs.

```{r compare-loadings-plot}
plot(b5_comparison, type = "loadings")
```

The summary's measurement evidence carries a note: omega for B in the
zero-loading model still includes `b5` in the composite, so it does not
describe a shortened four-item scale. To see that scale's reliability, fit the
model without `b5`. Because the two models contain different observed
variables, `nomo_compare()` reports this comparison descriptively, without a
difference test or information criteria:

```{r compare-shortened}
shortened_cfa <- nomo_cfa(
  "A =~ a1 + a2 + a3 + a4 + a5\nB =~ b1 + b2 + b3 + b4",
  data = nomo_demo_continuous
)

shortened <- nomo_compare(
  full = demo_cfa,
  four_item_b = shortened_cfa,
  rationale = "Describe the reliability of a four-item B scale alongside the full scale."
)

nomo_table(shortened, "comparisons")[, c("model", "relation", "test_available", "ic_available")]

omega_rows <- nomo_table(shortened, "evidence")
omega_rows[omega_rows$metric == "omega", c("model", "construct", "estimate")]
```

Together these results separate three questions that are easy to blur: whether
an item is statistically related to its factor, whether it is a strong
indicator, and how the score's reliability changes without it. None of them
alone decides whether the item stays.

# Does a result depend on how missing data were handled?

`demo_cfa` used lavaan's default, listwise deletion, which analyses only the
cases observed on every item. `nomo_missing()` refits the same model under
alternative strategies and reports what changes. For continuous indicators it
compares listwise deletion with full-information maximum likelihood (FIML),
which uses every observed response:

```{r missing-demo}
demo_missing <- nomo_missing(demo_cfa, data = nomo_demo_continuous)
demo_missing
```

The two strategies make different assumptions. Listwise deletion requires data
missing completely at random (MCAR). FIML requires the weaker condition that
data are missing at random (MAR), meaning that whether a value is missing may
depend on other observed values but not on the missing value itself. In Enders
and Bandalos's (2001) simulations, all of these methods were unbiased under
MCAR, and FIML was the most efficient. Under MAR, FIML alone remained unbiased
across the parameters they examined.

The missing values in `nomo_demo_continuous` were introduced completely at
random, so theory predicts agreement, and that is what the comparison shows.
Listwise deletion discards
`r demo_missing$strategies$n_discarded[demo_missing$strategies$strategy == "listwise"]`
cases. No estimate moves by as much as half a standard error. That is the size
Schafer and Graham (2002) treat as practically important for bias, because
beyond it confidence intervals noticeably lose coverage. Reliability barely
moves either:

```{r missing-demo-reliability}
nomo_table(demo_missing, "reliability")
```

## When the strategies disagree

Differences appear when data are missing at random but not completely at
random. In this simulation the population factor correlation is .50. Most cases
with high scores on factor A are missing every B item, so whether B is missing
depends on the observed A items:

```{r missing-mar}
set.seed(32)
n <- 4000
f <- matrix(rnorm(n * 2), n, 2) %*% chol(matrix(c(1, .5, .5, 1), 2))
items <- cbind(
  f[, 1] %o% rep(.7, 4) + matrix(rnorm(n * 4, sd = sqrt(1 - .7^2)), n),
  f[, 2] %o% rep(.7, 4) + matrix(rnorm(n * 4, sd = sqrt(1 - .7^2)), n)
)
mar_data <- as.data.frame(items)
names(mar_data) <- c(paste0("a", 1:4), paste0("b", 1:4))

high_a <- rowMeans(mar_data[, paste0("a", 1:4)]) > 0
drop <- high_a & runif(n) < .90
mar_data[drop, paste0("b", 1:4)] <- NA

mar_model <- "A =~ a1 + a2 + a3 + a4\nB =~ b1 + b2 + b3 + b4"
mar_missing <- nomo_missing(
  nomo_cfa(mar_model, data = mar_data),
  data = mar_data,
  reliability = FALSE
)

mar_estimates <- nomo_table(mar_missing, "estimates")
mar_estimates[mar_estimates$type == "factor_correlation",
              c("strategy", "estimate", "se", "difference_in_se")]
```

```{r missing-mar-values, include=FALSE}
mar_r <- mar_estimates[mar_estimates$type == "factor_correlation", ]
mar_listwise <- mar_r[mar_r$strategy == "listwise", ]
mar_fiml <- mar_r[mar_r$strategy == "ml", ]
se_from_truth <- function(row) sprintf("%.1f", abs(row$estimate - .50) / row$se)
```

FIML's estimate, `r sprintf("%.3f", mar_fiml$estimate)`, is
`r se_from_truth(mar_fiml)` of its standard errors from the population value,
within sampling error. Listwise deletion's,
`r sprintf("%.3f", mar_listwise$estimate)`, is `r se_from_truth(mar_listwise)`
standard errors below it. The complete cases under-represent high scorers on A,
and that restricted range attenuates the correlation. The decision log flags
the difference and says what it can and cannot mean:

```{r missing-mar-log}
mar_log <- nomo_table(mar_missing, "decision_log")
mar_log$recommendation[mar_log$metric == "estimate_difference"]
```

The flag carries two qualifications:

- The difference estimates listwise deletion's bias only if the data are MAR
  and the model is correct, because only then is FIML consistent.
- The two strategies analyze different cases, so part of any difference is
  sampling variability. The more cases listwise deletion discards, the larger
  that part becomes.

## What the comparison cannot tell you

Whether data are MAR cannot, in general, be tested from the data at hand
(Schafer & Graham, 2002). Agreement between strategies therefore shows only
that a result does not depend on the choice between them. It does not show
that either strategy is unbiased: when data are missing not at random, neither
is guaranteed to be. The strategy belongs in the analysis plan, decided in
advance, and `nomo_missing()` reports how much the results depend on it.

For ordered indicators estimated with WLSMV, lavaan does not offer FIML, so
`nomo_missing()` compares listwise with pairwise deletion. Both require MCAR,
so a difference between them shows that the result depends on the choice, but
it cannot be attributed to either strategy.

Mean substitution is not offered. Replacing each missing value with the item
mean shrinks variances and distorts covariances. Schafer and Graham (2002) show
that even under MCAR it narrows confidence intervals: with a quarter of values
missing, a nominal 95% interval for a mean covers the true value only about 86%
of the time. Multiple imputation, which Schafer and Graham recommend alongside
maximum likelihood, is a candidate for a later release.

# Reading the evidence as an argument

A useful measurement conclusion is rarely "all cutoffs passed." A stronger conclusion
identifies what each analysis contributes and what remains uncertain. The workflow is
therefore designed around four questions:

- **Structure:** Is the hypothesized factor model a reasonable representation of the data?
- **Precision:** Is the intended score sufficiently reliable for its use, conditional on that model?
- **Convergence:** Do indicators substantially represent the intended construct?
- **Separation:** Are constructs that theory distinguishes also distinguishable empirically?

The resulting decision logs preserve observations, references, recommendations, and
researcher decisions without automatically changing the scale or model.

# Research basis

The reliability workflow follows the literature recommending model-based omega for
congeneric measurement and cautioning against mechanical use of coefficient alpha,
including Dunn, Baguley, and Brunsden (2014), Flora (2020), and Bell, Chalmers, and
Flora (2024). Ordered-score reliability follows the distinction emphasized by Green
and Yang (2009) and implemented in `semTools`.

The validity workflow treats AVE as convergent evidence rather than reliability,
following Fornell and Larcker (1981). For construct separation, HTMT-family evidence
is prioritized because simulation research shows that older criteria can miss overlap
(Henseler, Ringle, & Sarstedt, 2015). HTMT2 is emphasized for congeneric measurement
because it relaxes the original HTMT tau-equivalence assumption (Roemer, Schuberth, &
Henseler, 2021). Fornell-Larcker output remains available only as optional historical/
supporting information.

The missing-data comparison follows Enders and Bandalos (2001), whose
simulations compared FIML with listwise and pairwise deletion under MCAR and
MAR, and Schafer and Graham (2002), who review the older methods and explain
why MAR cannot be tested from the data at hand.

The [research basis](research-basis.html) article places each of these methods
between historical and contemporary practice, and
`nomo_methods(stage = c("cfa", "compare", "reliability", "validity"))` returns
the same information as a table, with full references.