--- title: "Scoring a measurement model with nomo_scores()" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Scoring a measurement model with nomo_scores()} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} library(nomologR) knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE ) ``` After a measurement model is established, researchers usually need numbers to carry into later analyses. `nomo_scores()` produces them, and reports what they are and are not. Scoring is a modeling decision. `nomo_scores()` documents the decision; it does not make it. ## A model to score ```{r fit} items <- paste0("x", 1:8) model <- " Agency =~ x1 + x2 + x3 + x4 Persistence =~ x5 + x6 + x7 + x8 " set.seed(33) loadings <- cbind( Agency = c(.80, .75, .70, .65, 0, 0, 0, 0), Persistence = c(0, 0, 0, 0, .60, .55, .50, .45) ) phi <- matrix(c(1, .5, .5, 1), 2) sigma <- loadings %*% phi %*% t(loadings) diag(sigma) <- 1 dimnames(sigma) <- list(items, items) dat <- as.data.frame(matrix(rnorm(600 * 8), 600, 8) %*% chol(sigma)) names(dat) <- items fit <- lavaan::cfa(model, data = dat, std.lv = TRUE) ``` ## Adding items is a model, not arithmetic The most common way to score a scale is to add its items. McNeish and Wolf (2020) show that this is not a model-free calculation: adding items assumes a *parallel* model, in which the unstandardized loadings and the residual variances are equal across items. That assumption needs the same justification as any other measurement model. `nomo_scores()` fits that constrained model and compares it with the model you supplied: ```{r sum} summed <- nomo_scores(fit, method = "sum") summed ``` The loadings here were generated to differ, so the constraints are rejected. That does not forbid a sum score. It establishes that the items are not interchangeable in the way adding them assumes, so the choice needs a reason beyond convenience, and the diagnostics say what it costs. ## What a score is, and is not Grice (2001) evaluates scores on three criteria, all reported here and all computed from the fitted model rather than from the scores themselves: ```{r diagnostics} nomo_table(summed, "diagnostics") ``` `validity` is the correlation between a score and the factor it estimates. Gorsuch (1983) recommended at least .80, and above .90 where scores are used as substitutes for the factors themselves. For `method = "regression"` this quantity is the factor determinacy coefficient that [`nomo_hierarchical()`](hierarchical-models.html) reports, because that method maximizes it. `univocality` is a score's correlation with the factors it does *not* represent. A score that is not univocal cannot be treated as though it measured its own factor alone. `correlational_accuracy` is the one to read before using scores in a later analysis. It is the difference between the correlations among the scores and the correlations among the factors they stand in for. ## Scoring changes the answer Compare what each method implies about the relationship between the two constructs: ```{r comparison} comparison <- lapply(c("sum", "regression", "bartlett"), function(m) { s <- nomo_scores(fit, method = m) data.frame( method = m, factor_correlation = s$factor_correlations[1, 2], score_correlation = s$score_correlations[1, 2], discrepancy = s$diagnostics$correlational_accuracy[[1]] ) }) do.call(rbind, comparison) ``` The factor correlation is the same in every row, because it is a property of the model. The score correlation is not. Unit-weighted and Bartlett scores understate the relationship here, and regression scores overstate it. The direction is not a constant that could be corrected for: it depends on the scoring method and on the model. Grice reports the discrepancy running one way for orthogonal factors and the other for oblique ones. So `nomologR` reports it and leaves it visible rather than adjusting for it. The practical consequence is worth stating plainly. A relationship estimated from scores carries this discrepancy as bias. Where the question can be asked of the latent variables instead, asking it of scores replaces an unbiased answer with a biased one. ## One design recovers a regression There is one established exception, and its conditions are specific. For a linear regression of one factor on others, Skrondal and Laake (2001) proved that the regression coefficients are estimated consistently when * the predictors are scored with the regression method and the outcome with the Bartlett method, and * each block is scored from a measurement model of its own, not from one model containing both. Here is the design applied to the data above, with Agency predicting Persistence, next to the two ways of getting it wrong: ```{r skrondal-laake} agency_model <- lavaan::cfa("Agency =~ x1 + x2 + x3 + x4", data = dat, std.lv = TRUE) persistence_model <- lavaan::cfa("Persistence =~ x5 + x6 + x7 + x8", data = dat, std.lv = TRUE) slope <- function(outcome, predictor) unname(coef(lm(outcome ~ predictor))[2]) agency <- nomo_scores(agency_model, method = "regression")$scores$Agency persistence_bartlett <- nomo_scores(persistence_model, method = "bartlett")$scores$Persistence persistence_regression <- nomo_scores(persistence_model, method = "regression")$scores$Persistence joint_agency <- nomo_scores(fit, method = "regression")$scores$Agency joint_persistence <- nomo_scores(fit, method = "bartlett")$scores$Persistence data.frame( design = c( "Latent variables (the target)", "Regression then Bartlett, separate models", "Regression scores for both, separate models", "Regression then Bartlett, one joint model" ), slope = c( lavaan::lavInspect(fit, "cov.lv")["Agency", "Persistence"], slope(persistence_bartlett, agency), slope(persistence_regression, agency), slope(joint_persistence, joint_agency) ) ) ``` Both factors have a variance of 1, so the target slope is their covariance in the fitted model. The design lands next to it. Using the regression method for both blocks falls well short, and scoring both factors from the one two-factor model overshoots, although it uses the right method for each. Each condition matters. Three limits come with the result. It is about regression coefficients, not the correlations `correlational_accuracy` describes. The standard errors `lm()` reports treat the scores as observed data, and Skrondal and Laake note that corrected standard errors and confidence intervals may require resampling. And it does not extend to nonlinear models, for which they caution that factor score regression appears to perform very badly. `nomologR` does not apply the design for you, because it depends on which factor is the outcome, which is part of the research question. Where the question can be asked of the latent variables themselves, that remains the simpler route. ## Choosing a method ```{r methods} nomo_table(nomo_scores(fit, method = "regression"), "diagnostics") ``` With one factor, regression and Bartlett scores differ only by a scaling constant and rank people identically, so the choice between them matters only when there is more than one factor. Reliability coefficients follow the weighting. Coefficient alpha is the reliability of a unit-weighted scale; coefficient H, reported by [`nomo_hierarchical()`](hierarchical-models.html), belongs to optimally weighted scores. Reporting one for the other describes a scale that was not used. ## What is never done The data you supply is never modified, and scores are returned only for the cases the model actually used, so a case dropped from the fit does not acquire a score. No threshold here is applied as a rule. Gorsuch's recommendations appear where a value falls below them, labeled as his, in the same way fixed fit-index cutoffs are treated elsewhere in the package. ## References Gorsuch, R. L. (1983). *Factor analysis* (2nd ed.). Lawrence Erlbaum. Grice, J. W. (2001). Computing and evaluating factor scores. *Psychological Methods, 6*(4), 430–450. McNeish, D., & Wolf, M. G. (2020). Thinking twice about sum scores. *Behavior Research Methods, 52*(6), 2287–2305. Skrondal, A., & Laake, P. (2001). Regression among factor scores. *Psychometrika, 66*(4), 563–575.