Overview

The goal of this exercise is for you to become more familiar with fitting basic structural equation models using the lavaan R package. We will start first with a very simple dataset consisting of SAT and ACT scores collected in 700 individuals. To learn additional details, see: ?psych::sat.act.

Here is a quick description of the data:

##             n   mean     sd median min max  skew kurtosis
## sex*      700   1.65   0.48      2   1   2 -0.61    -1.62
## education 700   3.16   1.43      3   0   5 -0.68    -0.07
## age       700  25.59   9.50     22  13  65  1.64     2.42
## ACT       700  28.55   4.82     29   3  36 -0.66     0.53
## SATV      700 612.23 112.90    620 200 800 -0.64     0.33
## SATQ      687 610.22 115.64    620 200 800 -0.59    -0.02

As you can see, there are SAT verbal scores (SATV), quantitative scores (SATQ), and overall ACT scores for 700 people (a bit of missingness on SATQ).

Test a single-factor of achievement

What if we believe that scores on the three achievement tests reflect something like an overall measure of academic aptitude/achievement? How would you fit a one-factor CFA to these data?

#one-factor CFA here

Standardized estimates

#How can you print out the standardized factor loadings here?

Fit diagnosis

How well does the model fit according to global indices including CFA, SRMR, RMSEA, and \(\chi^2\)? If you know SEM: is the model overidentified, underidentified, or just-identified? How does this qualify your interpretation of fit?

#Print out and evaluate global model fit

Predictors of achivement

Assuming that SAT and ACT scores were recorded first (i.e., SAT and ACT came before post-secondary education), does higher academic achievement (our latent factor) predict higher ultimate educational level?

Additional questions

Do the lavaan model warnings worry you?

How strong is the association (think correlation or standardized regression coefficient)?

What about missing data? How is it being handled, and could you do better?

#modifications and details here

Sex differences in achievement

Are there significant mean differences in achievement scores between men and women? If so, interpret the direction of the effect.

#recode to make the direction and name clearer in the output
sat.act$female <- as.numeric(sat.act$sex=="female")

#fit model allowing for sex differences

Bonus question: moderation by sex

If you like, examine whether the association between achievement scores and ultimate educational attainment differs between men and women. Note that this is best done as a multiple-groups SEM using the group= argument in lavaan (specifically, when calling the sem function). Also note that you may wish to control which parameters are free to differ between sexes and which must be equal using group.equal. See ?lavOptions for details. Also examine the modification indices to see if model fit could be improved.

A more complex model

Next, let’s look at a more detailed dataset that allows for more complex modeling of latent structure. We’ll just be dealing with measurement models here – that is, variants of CFA. This is the iqitems dataset from the psych package, which has questions from different subscales of an IQ test.

describe(iqitems) %>% select(-vars, -trimmed, -mad, -se, -range)
##              n mean   sd median min max  skew kurtosis
## reason.4  1523 3.39 1.24      4   0   6 -1.28     1.25
## reason.16 1524 3.39 1.14      4   0   6 -1.61     1.70
## reason.17 1523 3.80 1.34      4   0   6 -1.17     2.08
## reason.19 1523 4.79 1.84      6   0   6 -1.28     0.24
## letter.7  1524 4.95 1.70      6   0   6 -1.72     1.98
## letter.33 1523 2.79 1.25      3   0   6 -0.06     0.64
## letter.34 1523 3.38 1.30      4   0   6 -1.15     0.65
## letter.58 1525 3.26 1.52      4   0   6 -0.71    -0.61
## matrix.45 1523 4.14 1.36      5   0   6 -1.44     1.71
## matrix.46 1524 2.47 1.36      2   0   6  0.97     0.56
## matrix.47 1523 2.59 1.37      2   0   6  1.01     0.63
## matrix.55 1524 3.69 1.56      4   0   6 -0.31    -0.30
## rotate.3  1523 4.65 2.16      4   0   8 -0.04    -0.64
## rotate.4  1523 4.77 2.47      4   0   8 -0.18    -1.32
## rotate.6  1523 4.40 2.56      5   0   8 -0.26    -1.26
## rotate.8  1524 4.63 2.40      4   0   8 -0.14    -1.25

Note that the names of the items denote the respective IQ subscale from which they came. Thus, we are interested in whether we can first corroborate that there are four components of intelligence: mental rotation, matrix reasoning, letter sequences, and basic verbal reasoning. There are four items from each domain.

ggcorrplot(lavCor(iqitems))

Fit the a priori four-factor CFA

First, fit the four-factor model based on the expected constructs, printing out the standardized loading in the summary.

Diagnosing problems

Does the warning provided by lavaan give you reason for concern? If so, what’s the diagnosis? (Hint: use a variant of inspect aka lavInspect.)

How well does the model fit?

How would you interpret the quality of global fit based on fit indices?

Fit a simpler model

Based on the pattern of observed correlations and the results of the initial four-factor model, conceptualize and fit a simpler model with fewer factors.

Does this model yield a better fit? Hint: you can use the anova() function to compare alternative factor models using a likelihood ratio test.

Finding an even better model

Residuals

Look at the residual correlations of the model. What correlations are being poorly estimated? Does this inform your assessment of how to improve the model?

Modification indices

Do the modification indices suggest any plausible parameters that are omitted from the model? Alternatively, do the modification indices suggest problems with particular items?

Can you improve the model further? Consider:

  • dropping items that aren’t loading well on any factor
  • shifting items to another factor if they don’t load where we expected
  • aggregating items that seem highly and uniquely overlapping

Bonus question: fit a hierarchical factor model

In principle, there may be a general factor of intelligence that explains the correlations among the subscales here. Fit a hierarchical factor model in which ‘g’ is a superordinate factor that explains the lower-order factors. Does this model fit significantly worse than a model in which all lower-order factors simply correlate, but there is no ‘g’?

(Noting that the four-factor model is bad.) BIC is indifferent between the alternatives, the LRT says model fit is significantly worse for hierarchical (but we have a lot of data, so a lot of power to reject the null). AIC is 11 points lower for the four-factor model. Altogether, this is equivocal support for a four-factor model. We could debate which one is ‘better,’ but there is not a clear winner.