--- title: "EFA with ordinal and missing data" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{EFA with ordinal and missing data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.align = "center" ) # The multiple-imputation section uses mice to create the imputations. It is only # suggested by the package, so the imputation chunks are evaluated only when it is # installed. mice_ok <- requireNamespace("mice", quietly = TRUE) ``` Two features of real data complicate an exploratory factor analysis: items are often **ordinal** (a handful of Likert categories rather than a continuous scale), and some responses are usually **missing**. `EFAtools` handles both within the ordinary `efa_fit()` workflow, without switching packages. This vignette shows how. It assumes familiarity with the basic workflow covered in the [EFAtools](EFAtools.html) vignette and focuses on what changes for ordinal and incomplete data. ```{r} library(EFAtools) ``` So that the examples are self-contained and reproducible, we generate the data with `efa_simulate()` from a known three-factor population (18 indicators, six per factor, with moderately correlated factors), using fixed seeds throughout. ```{r} Lambda <- population_models$loadings$baseline # 18 x 3 loading pattern Phi <- population_models$phis_3$moderate # moderate factor intercorrelations ``` ## Ordinal Data Rating-scale items are not continuous: they take a few ordered values, and a Pearson correlation between two such items underestimates the association between the underlying constructs. The polychoric correlation instead estimates the correlation of the continuous latent variables assumed to underlie the observed categories, and pairing it with a categorical estimator removes the bias that treating the items as continuous introduces. We draw 400 responses on a four-category scale. Because the latent data are normal, cutting them at the standard-normal category thresholds already leaves the population *polychoric* correlation of the discretised data equal to the target correlation; `match = "polychoric"` records that this is what we are after. ```{r} d_ord <- efa_simulate(N = 400, Lambda = Lambda, Phi = Phi, categories = 4, match = "polychoric", seed = 2024)$data d_ord[1:5, 1:6] ``` ### Screening Ordinal Data `efa_screen()` reports, among its diagnostics, how many response categories each item has and whether the data are multivariate normal — the two things that decide whether an ordinal treatment is worthwhile. ```{r, warning = FALSE} efa_screen(d_ord, seed = 42) ``` The sampling adequacy (KMO) and sphericity checks confirm the data are factorable. The telling parts are the multivariate-normality section and the recommendations: Mardia's kurtosis and the Henze-Zirkler test reject normality, and the recommendations flag that every item has fewer than five response categories. Together, these point to the same conclusion. With few categories and non-normal data, a polychoric correlation with a categorical estimator such as DWLS is less biased than normal-theory maximum likelihood. The normal-theory standard errors and fit indices are also better replaced by robust (sandwich) versions. ### Polychoric Correlations, DWLS, and Robust Standard Errors `efa_fit()` computes the polychoric correlation when `cor_method = "poly"` and fits it with diagonally weighted least squares when `estimator = "DWLS"` — the estimator recommended for ordinal data because it accounts for how precisely each polychoric correlation is estimated. Requesting `se = "sandwich"` adds robust standard errors and a scaled (Satorra-Bentler) chi-square that stay valid under the non-normality these data show. ```{r} efa_poly <- efa_fit(d_ord, n_factors = 3, cor_method = "poly", estimator = "dwls", rotation = "oblimin", se = "sandwich") efa_poly ``` The pattern matrix recovers the three factors cleanly (six indicators each), and the model fit reports a **scaled** chi-square with its CFI, TLI, and RMSEA. Because the chi-square is a scaled statistic, the AIC and BIC (which are defined on the unscaled likelihood discrepancy) are left `NA`. For binary items, `cor_method = "tetra"` computes tetrachoric correlations and runs the same DWLS and sandwich machinery. The robust standard errors accompany each estimated quantity; for the rotated loadings, for example: ```{r} round(efa_poly$SE$rot_loadings, 3) ``` The matching confidence intervals live in `efa_poly$CI`, and `summary(efa_poly)` prints them as a labelled table alongside the model diagnostics. ### Why Not Just Treat the Items as Continuous? To see what the ordinal treatment buys, fit the same data as if they were continuous — a Pearson correlation with maximum likelihood — and compare the rotated loadings with `efa_compare()`. ```{r} efa_cont <- efa_fit(d_ord, n_factors = 3, cor_method = "pearson", estimator = "ML", rotation = "oblimin") cmp <- efa_compare(efa_poly$rot_loadings, efa_cont$rot_loadings, x_labels = c("Polychoric / DWLS", "Pearson / ML")) cmp plot(cmp) ``` The two solutions agree on the structure, but the polychoric loadings are systematically a little larger: treating the items as continuous attenuates the loadings, because the Pearson correlation understates the latent associations. With only four categories here the gap is modest, but it widens as the number of categories drops (it is largest for binary items) and as the category thresholds grow more asymmetric (skewed items). This is why a polychoric or tetrachoric treatment is preferable for genuinely ordinal items with few categories. The polychoric route does make its own demands, though: it assumes a normal latent variable underlies each item, and it needs an adequate sample size and reasonably populated response-category combinations. When categories are very sparse (rare responses, small samples), the polychoric asymptotic covariance behind the DWLS weights and the robust standard errors becomes unreliable, and collapsing rare categories can help. ## Missing Data When some responses are missing, dropping every incomplete case (listwise deletion) wastes data and can bias the results unless the values are missing completely at random. `EFAtools` offers two principled alternatives that assume only that the data are missing at random (MAR): a single-analysis route via full-information maximum likelihood, and a multiple-imputation route via `efa_mi()`. We simulate 250 continuous cases with about 15% of values missing at random, where each item's missingness depends on another item's value. By default, `efa_simulate()` would make each item's missingness depend on another item that is itself partly missing — a kind of missingness the two routes below are not designed for. Rather than failing outright, they run with a residual bias that grows as `missing_prop` and `missing_strength` increase. `missing_vars` and `missing_predictor` avoid that. Here the first nine items carry the missing values and each is driven by one of the last nine, which stay complete. Every predictor is therefore fully observed, satisfying the *ignorable* MAR assumption the two routes below rely on -- not just MAR, but MAR with no incomplete predictor left unaccounted for. ```{r} d_miss <- efa_simulate(N = 250, Lambda = Lambda, Phi = Phi, missing = "MAR", missing_prop = 0.15, missing_vars = 1:9, missing_predictor = 10:18, seed = 2024)$data round(mean(is.na(d_miss)), 3) # overall proportion missing round(colMeans(is.na(d_miss)), 3) # holed items only ``` ### Two-Stage Full-Information Maximum Likelihood With `cor_method = "fiml"`, `efa_fit()` estimates the saturated mean and covariance from all the observed data by an EM algorithm (assuming the data are MAR) and analyses the resulting correlation — a single fit that uses every case rather than only the complete ones. The model fit is reported as corrected two-stage (Satorra-Bentler) statistics. ```{r} efa_fiml <- efa_fit(d_miss, n_factors = 3, cor_method = "fiml", estimator = "ml", rotation = "oblimin") efa_fiml ``` The solution again recovers the three factors, and the printout records that the correlation was obtained by two-stage FIML. Standard errors are available here too: for `estimator = "ML"` or `"ULS"`, `se = "information"` or `"sandwich"` return the corrected two-stage standard errors, and `se = "np-boot"` works with any estimator. FIML estimates the model from every case, but it does not fill in the missing values. Any step that needs complete rows is therefore still complete-case. `efa_scores(d_miss, f = efa_fiml)`, for example, scores only the complete cases and warns that the rest are `NA`: ```{r} sum(complete.cases(d_miss)) # cases a score can be formed for ``` Use the `efa_mi()` route below if you want case-level scores for every respondent. Each imputed dataset is complete, so scoring it with the pooled loadings covers all 250 cases. The correction itself needs a well-behaved saturated covariance. When that cannot be formed — typically in a small sample with a high proportion of missing values, or with near-collinear items — `efa_fit()` warns and keeps the plain two-stage likelihood-ratio statistic instead of discarding the test. Such a fit is never presented as corrected: the printed chi-square line is labelled **uncorrected**, and the accompanying *p*-value, CFI, TLI, and RMSEA should be read as indicative only. ### Multiple Imputation with `efa_mi()` The alternative is to impute the missing values several times, fit each completed dataset, and pool the results. `EFAtools` does not impute the data itself — use a dedicated tool such as the [mice](https://CRAN.R-project.org/package=mice) package — but `efa_mi()` takes the list of completed datasets and does the factor-analytic pooling. Here we create five imputations with mice (a Bayesian linear-regression model, appropriate for these continuous items) and collect them into a list. ```{r, eval = mice_ok} imp <- mice::mice(as.data.frame(d_miss), m = 5, method = "norm", printFlag = FALSE, seed = 123) dat_list <- lapply(seq_len(imp$m), function(i) mice::complete(imp, i)) ``` `efa_mi()` fits the same `efa_fit()` model to each imputed dataset — the extraction, rotation, and standard-error options are passed through `...` — aligns the solutions to a common factor space (rotation is only identified up to reflection and permutation, so the imputations must be matched before averaging), and pools them. ```{r, eval = mice_ok} efa_pooled <- efa_mi(dat_list, n_factors = 3, estimator = "ml", rotation = "oblimin") efa_pooled ``` The pooled loadings recover the three factors. Point estimates are averaged across the imputations after alignment. The model chi-square, and the AIC/BIC derived from it, are pooled with the same rule as RMSEA, but applied to a different discrepancy scale — which is why the printout labels the pooled chi-square **D2-pooled**, and why you should not expect it to reconcile by hand with the pooled RMSEA. The incremental CFI and TLI are instead averaged across the per-imputation fits. Requesting standard errors in the call (for example `se = "information"` or `se = "np-boot"`) additionally pools them with Rubin's rules, so the between-imputation variability inflates the pooled standard errors. Because multiple imputation propagates the extra uncertainty from the missing data, its pooled fit statistics are not directly comparable with the single FIML fit above; read them together with the per-imputation fits stored in the returned object. Which route to prefer is largely practical. FIML is a single, efficient fit and is the simpler default when the analysis model is the whole story. Multiple imputation is more flexible when the imputation model should draw on auxiliary variables not in the factor model, or when the same imputations feed several downstream analyses. ## Ordinal *and* Missing Data Questionnaire data are usually both: a few response categories *and* some unanswered items. The two treatments above do not simply combine, because the ordinal machinery needs complete cases exactly where the FIML route needs continuous ones. We draw the same four-category items as before, now with 8% of the responses missing completely at random. ```{r} d_ord_miss <- efa_simulate(N = 300, Lambda = Lambda, Phi = Phi, categories = 4, match = "polychoric", missing = "MCAR", missing_prop = 0.08, seed = 2024)$data round(mean(is.na(d_ord_miss)), 3) # overall proportion missing sum(complete.cases(d_ord_miss)) # respondents who answered every item ``` Eight percent per item is mild, but spread over 18 items it leaves only a fifth of the sample complete. That number, not the 300 rows supplied, is what the ordinal route has to work with. ### The Effective Sample of a Polychoric DWLS Fit Polychoric correlations can be estimated pair by pair, using every respondent who answered a given pair of items. Their asymptotic covariance cannot: the DWLS weights and the sandwich standard errors describe *one* set of estimates from *one* set of cases, so whenever they are requested the correlations and the covariance are both computed on the listwise-complete rows. `efa_fit()` announces that override rather than applying it silently. ```{r} efa_ord_miss <- efa_fit(d_ord_miss, n_factors = 3, cor_method = "poly", estimator = "dwls", rotation = "oblimin") efa_ord_miss$settings$N # cases the fit is actually based on ``` Every loading, fit index, standard error, and confidence interval of that fit rests on the complete cases. `efa_screen()` names both denominators for the same reason — its multivariate-normality and outlier results state how many complete cases they used out of how many rows were supplied — so the reduction is visible before a model is fitted. Two adjustments avoid it, each with a cost: - Ask for no asymptotic covariance. With `estimator = "ULS"` or `"ML"` and no sandwich standard errors, nothing needs the weight matrix, so the polychoric matrix stays pairwise-complete and every case contributes to the pairs it answered. In exchange you give up the DWLS weighting and the robust standard errors, and the reported `N` is the nominal row count while each correlation rests on its own smaller subset — so fit statistics and analytic standard errors treat the data as more complete than they are. - Switch to FIML — but only as a continuous analysis. `cor_method = "fiml"` estimates a saturated *normal* covariance and analyses the resulting Pearson-type correlation; there is no polychoric asymptotic covariance in it, so it cannot supply DWLS with weights. Combining the two is refused rather than silently approximated: ```{r, error = TRUE} efa_fit(d_ord_miss, n_factors = 3, cor_method = "fiml", estimator = "dwls") ``` ### Keeping Both the Ordinal Treatment and the Cases Multiple imputation is the route that keeps both. Each imputed dataset is complete, so the polychoric correlations and their asymptotic covariance are estimated on all 300 cases, and the imputation uncertainty is carried into the pooled results. Since `efa_mi()` forwards its arguments to `efa_fit()`, this is the ordinal fit above applied to each imputation — here with predictive mean matching, which resamples observed values and so respects the response scale. ```{r, eval = mice_ok} imp_ord <- mice::mice(as.data.frame(d_ord_miss), m = 5, method = "pmm", printFlag = FALSE, seed = 123) dat_ord <- lapply(seq_len(imp_ord$m), function(i) mice::complete(imp_ord, i)) efa_ord_pooled <- efa_mi(dat_ord, n_factors = 3, cor_method = "poly", estimator = "dwls", rotation = "oblimin") efa_ord_pooled ``` The pooled pattern matrix is estimated from all 300 respondents rather than the 62 complete ones — compare it with `efa_ord_miss$rot_loadings` — and the between-imputation variability enters any pooled standard errors requested in the call. As always with imputation, the imputation model has to be defensible: for ordinal items that means a method returning values on the observed scale (predictive mean matching, or an ordinal model such as `method = "polr"`), and enough respondents per item to estimate it. So, for ordinal items with missing values: use `efa_mi()` when the missingness is more than incidental and the ordinal treatment matters; fit `cor_method = "poly"` with `estimator = "DWLS"` directly when the complete-case sample is still comfortably large; and use FIML when the items have enough categories to be treated as continuous in the first place. ## Where to Next This vignette covered the ordinal and missing-data extensions of the workflow. For the core analysis — screening, factor retention, extraction and rotation, and the post-processing tools — see the [EFAtools](EFAtools.html) vignette, and the individual help pages for the statistical details and references. Run `browseVignettes("EFAtools")` for the vignettes installed with the package, or visit the [package website](https://mdsteiner.github.io/EFAtools/) for the full set of articles.