Parcels and second-order CFA

Parcel construction and second-order factors with lavaan

Author
Affiliation

George Michaelides

Norwich Business School, University of East Anglia

Published

July 7, 2026

Learning objectives

By the end of this tutorial you will be able to:

  • explain what a parcel is and why researchers use them
  • construct parcels from individual items using base R
  • fit a CFA on parceled indicators instead of individual items
  • specify and fit a second-order factor model in lavaan
  • decide when parceling is appropriate and when it is not

Setup

Setup code (you have seen this before) 
library(haven)
library(lavaan)
library(tidySEM)

What is a parcel?

A parcel is an aggregate — usually a mean or sum — of two or more items, used as a single indicator of a latent variable. Instead of specifying a factor with six individual items as indicators, you might group those six items into three parcels of two items each, and specify the factor with three parcel indicators.

Why parcel?

Researchers parcel for several reasons:

  1. Fewer indicators per factor. A model with 30 individual items has many parameters to estimate and many degrees of freedom. Parceling the same 30 items into 10 parcels reduces both, often producing a model that converges more easily and fits better.

  2. Reduced measurement error. Averaging two or more items reduces the random error in each individual item, so parcels are more reliable than individual items.

  3. Continuous-ish indicators. A mean of two 5-point Likert items can take values like 3.5, which behaves more like a continuous variable than an individual item with only 5 possible values. This makes the normality assumption of ML estimation more plausible.

Why not parcel?

Parceling is controversial:

  1. Loss of item-level information. Once you average items into a parcel, you cannot test whether individual items load well on the factor. A parcel can mask a badly-functioning item.

  2. Parcels are not interpretable. A parcel made of con1 and con10 does not correspond to any single questionnaire item, so you cannot say what a high score on that parcel “means” in the way you can for an individual item.

  3. It can hide a multidimensional structure. If your six items actually measure two different sub-factors, parceling them into three parcels of two items each can hide this — the parcels will all appear to load on a single factor even though the underlying structure is two-dimensional.

WarningParcel only when justified

Parcel only when (a) you have a theory that says the items measure a single unidimensional construct, and (b) you have checked this with an item-level CFA or EFA first. If an item-level CFA fits well, parceling is a reasonable simplification for a larger SEM. If an item-level CFA fits poorly, parceling will mask the problem rather than solve it.

Parcel construction

There are two main strategies for constructing parcels, both of which appear in the SEM literature.

Domain-representative parcels

Each parcel contains items from across the sub-domains of the construct. If an EFA suggests the items split into two sub-factors, you build parcels that each contain one item from each sub-factor. This ensures each parcel represents the full breadth of the construct.

Homogeneous parcels

Each parcel contains items from a single sub-domain. This is simpler but produces parcels that may reflect sub-factors rather than the overall construct.

Building parcels in R

We use base R — rowMeans() is the cleanest way to compute parcel scores from a data frame. No dplyr needed.

We will parcel the six job-control items (con1, con3, con6, con8, con9, con10) into three parcels of two items each, using a domain-representative strategy. The pairing is arbitrary for this teaching example; in real use you would pair items based on an EFA or theoretical reasoning.

ps$p1 <- (ps$con9 + ps$con10) / 2
ps$p2 <- (ps$con1 + ps$con8)  / 2
ps$p3 <- (ps$con3 + ps$con6)  / 2

Each parcel is the mean of two items. The new columns p1, p2, p3 are added to the ps data frame and can be used as indicators in a CFA.

NoterowMeans() for larger parcels

For parcels of three or more items, rowMeans() is cleaner than summing and dividing manually:

ps$p1 <- rowMeans(ps[, c("con9", "con10", "con1")], na.rm = TRUE)

The na.rm = TRUE argument ensures that a participant with one missing item still gets a parcel score from the remaining items. Be careful: if all items are missing, rowMeans() returns NaN, not NA. You may want to recode NaN to NA after parceling.

A CFA on parceled indicators

Now we fit a one-factor CFA on the three job-control parcels:

model_parcel <- '
  con =~ p1 + p2 + p3
'
fit_parcel <- cfa(model_parcel, data = ps, estimator = "MLR")
summary(fit_parcel, standardized = TRUE, fit.measures = TRUE)
## lavaan 0.6-21 ended normally after 17 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         6
## 
##   Number of observations                           343
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                                 0.000       0.000
##   Degrees of freedom                                 0           0
## 
## Model Test Baseline Model:
## 
##   Test statistic                               389.047     346.953
##   Degrees of freedom                                 3           3
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  1.121
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.000       1.000
##                                                                   
##   Robust Comparative Fit Index (CFI)                            NA
##   Robust Tucker-Lewis Index (TLI)                               NA
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1017.748   -1017.748
##   Loglikelihood unrestricted model (H1)      -1017.748   -1017.748
##                                                                   
##   Akaike (AIC)                                2047.496    2047.496
##   Bayesian (BIC)                              2070.522    2070.522
##   Sample-size adjusted Bayesian (SABIC)       2051.489    2051.489
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000          NA
##   90 Percent confidence interval - lower         0.000          NA
##   90 Percent confidence interval - upper         0.000          NA
##   P-value H_0: RMSEA <= 0.050                       NA          NA
##   P-value H_0: RMSEA >= 0.080                       NA          NA
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.000
##   P-value H_0: Robust RMSEA <= 0.050                            NA
##   P-value H_0: Robust RMSEA >= 0.080                            NA
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.000       0.000
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   con =~                                                                
##     p1                1.000                               0.690    0.773
##     p2                0.803    0.066   12.254    0.000    0.554    0.797
##     p3                0.902    0.072   12.602    0.000    0.622    0.795
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .p1                0.320    0.038    8.443    0.000    0.320    0.402
##    .p2                0.177    0.021    8.599    0.000    0.177    0.366
##    .p3                0.226    0.029    7.863    0.000    0.226    0.369
##     con               0.475    0.061    7.734    0.000    1.000    1.000

Reading the output

The structure is identical to the CFA we fitted in Tutorial 07 — only the indicators have changed. The key differences:

  • Fewer parameters. A one-factor CFA with three indicators has 6 free parameters (2 loadings + 1 latent variance + 3 residual variances) and \(\text{df} = 6 - 6 = 0\). The model is just-identified — it fits perfectly by construction.
  • Higher loadings. Parcels are more reliable than individual items, so the standardised loadings on the factor are usually higher — often above .80.

The just-identification problem

A one-factor CFA with three parcels is just-identified — it has zero degrees of freedom and fits perfectly. This means we cannot test whether the factor model is correct. The parcels will always appear to load on a single factor, even if the underlying items have a multidimensional structure.

This is the central tension of parceling: you gain a simpler, better-fitting model, but you lose the ability to test the measurement model. For this reason, parceling is most defensible when:

  1. You have already confirmed the item-level factor structure in a separate CFA, and
  2. You are using the parcels in a larger SEM (Tutorial 10) where the measurement model is not the primary focus.
TipTry it

Parcel the four problem-solving demands items (prob1, prob2, prob3, prob4) into two parcels of two items each. Fit a one-factor CFA on the two parcels. What happens to the degrees of freedom? Why?

ps$prob_p1 <- (ps$prob1 + ps$prob2) / 2
ps$prob_p2 <- (ps$prob3 + ps$prob4) / 2

model_prob <- '
  prob =~ prob_p1 + prob_p2
'
fit_prob <- cfa(model_prob, data = ps, estimator = "MLR")
fitMeasures(fit_prob, c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))
##  chisq     df pvalue    cfi    tli  rmsea   srmr 
##     NA     -1     NA     NA     NA      0      0

A one-factor CFA with two indicators has \(q = 2\), so \(q(q+1)/2 = 3\), and \(t = 3\) (1 loading + 1 latent variance + 2 residual variances), giving \(\text{df} = 0\). The model is just-identified. This is the same problem as a one-factor CFA with two individual items — parceling has not helped. To get an overidentified model, you need at least three parcels per factor.

Second-order factor models

A second-order factor is a latent variable whose indicators are themselves latent variables. Instead of items loading directly on a single factor, items load on first-order factors, and those first-order factors load on a second-order factor.

When to use a second-order factor

Second-order factors are used when theory suggests that several related but distinct constructs are all manifestations of a single higher-order construct. A classic example is burnout, measured by two first-order factors — emotional exhaustion and cynicism — which are both manifestations of the higher-order construct “burnout”.

In our PS dataset, we can model affect as a second-order factor measured by two first-order factors: positive affect (PA) and negative affect (neg). The theory is that PA and neg are distinct but related dimensions of an underlying “affect” construct.

Specifying a second-order factor in lavaan

The syntax extends naturally. We define the first-order factors as usual, then define the second-order factor as measured by the first-order factors:

model_2nd <- '
  PA     =~ pa1 + pa2 + pa3
  neg    =~ na1 + na2 + na3
  affect =~ PA + neg
'
fit_2nd <- cfa(model_2nd, data = ps, estimator = "MLR")

The line affect =~ PA + neg says “the second-order factor affect is measured by the first-order factors PA and neg”. The first-order factors are now endogenous — they are predicted by the second-order factor — and their residual variances (\(\zeta\), disturbances) are estimated.

Identification of second-order factors

The identification rules for second-order factors are the same as for first-order factors:

  • The second-order factor’s scale must be set (by the marker-variable method or std.lv = TRUE)
  • Each second-order factor should have at least three first-order indicators

Our affect second-order factor has only two first-order indicators (PA and neg). This makes the second-order portion just-identified — the relationship between affect, PA, and neg fits perfectly. The overall model’s degrees of freedom come entirely from the first-order measurement model.

To have an overidentified second-order portion, you would need at least three first-order factors loading on the second-order factor — for example, if you had PA, neg, and a third factor like “activation” loading on affect.

NoteTwo first-order factors is common in practice

Many published second-order models have only two first-order factors per second-order factor. This is acceptable when the second-order factor is not the primary focus of the analysis and the overall model is overidentified through the first-order portion. Just be aware that the second-order portion itself is not testable with only two indicators.

Reading the output

summary(fit_2nd, standardized = TRUE, fit.measures = TRUE)
## lavaan 0.6-21 ended normally after 23 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        14
## 
##                                                   Used       Total
##   Number of observations                           330         343
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                                 9.870       8.801
##   Degrees of freedom                                 7           7
##   P-value (Chi-square)                           0.196       0.267
##   Scaling correction factor                                  1.121
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              1009.810     961.135
##   Degrees of freedom                                15          15
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  1.051
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.997       0.998
##   Tucker-Lewis Index (TLI)                       0.994       0.996
##                                                                   
##   Robust Comparative Fit Index (CFI)                         0.998
##   Robust Tucker-Lewis Index (TLI)                            0.996
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -2610.126   -2610.126
##   Scaling correction factor                                  0.885
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -2605.191   -2605.191
##   Scaling correction factor                                  0.964
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                                5248.253    5248.253
##   Bayesian (BIC)                              5301.440    5301.440
##   Sample-size adjusted Bayesian (SABIC)       5257.032    5257.032
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.035       0.028
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.082       0.074
##   P-value H_0: RMSEA <= 0.050                    0.641       0.734
##   P-value H_0: RMSEA >= 0.080                    0.057       0.029
##                                                                   
##   Robust RMSEA                                               0.030
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.081
##   P-value H_0: Robust RMSEA <= 0.050                         0.680
##   P-value H_0: Robust RMSEA >= 0.080                         0.056
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.015       0.015
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   PA =~                                                                 
##     pa1               1.000                               0.973    0.829
##     pa2               1.000    0.061   16.280    0.000    0.972    0.764
##     pa3               1.086    0.067   16.105    0.000    1.056    0.856
##   neg =~                                                                
##     na1               1.000                               0.948    0.848
##     na2               0.935    0.058   16.071    0.000    0.887    0.819
##     na3               0.983    0.056   17.668    0.000    0.932    0.836
##   affect =~                                                             
##     PA                1.000                                  NA       NA
##     neg               0.802    0.097    8.280    0.000       NA       NA
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .pa1               0.431    0.053    8.118    0.000    0.431    0.313
##    .pa2               0.676    0.071    9.560    0.000    0.676    0.417
##    .pa3               0.406    0.063    6.480    0.000    0.406    0.267
##    .na1               0.350    0.046    7.659    0.000    0.350    0.280
##    .na2               0.386    0.050    7.798    0.000    0.386    0.329
##    .na3               0.375    0.043    8.677    0.000    0.375    0.302
##    .PA                1.447    0.120   12.110    0.000    1.530    1.530
##    .neg               1.222    0.148    8.245    0.000    1.359    1.359
##     affect           -0.501    0.023  -21.818    0.000       NA       NA

The output now has an additional layer:

  • First-order loadings (PA =~ pa1, pa2, pa3; neg =~ na1, na2, na3) — identical to the two-factor CFA from Tutorial 07
  • Second-order loadings (affect =~ PA, neg) — the loadings of the first-order factors on the second-order factor
  • First-order residual variances (\(\zeta\) for PA and neg) — the variance in each first-order factor not explained by the second-order factor

The second-order loadings tell you how strongly each first-order factor reflects the higher-order construct. The first-order residual variances tell you how much of each first-order factor is not explained by the second-order factor — a large residual variance means the first-order factor has substantial unique variance beyond the higher-order construct.

Fit indices

fitMeasures(fit_2nd, c("chisq", "df", "pvalue",
                       "cfi", "tli",
                       "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                       "srmr"))
##          chisq             df         pvalue            cfi            tli 
##          9.870          7.000          0.196          0.997          0.994 
##          rmsea rmsea.ci.lower rmsea.ci.upper           srmr 
##          0.035          0.000          0.082          0.015

Compare these to the fit indices from the two-factor CFA in Tutorial 07. The second-order model has the same number of observed variables (6) but fewer free parameters — the covariance between PA and neg has been replaced by a structural relationship through affect. This gives the second-order model more degrees of freedom. If the second-order structure is correct, the fit should be nearly as good as the first-order model. If the second-order structure is wrong (PA and neg do not share a common higher-order cause), the fit will be worse.

Comparing nested models

We can formally compare the second-order model to the first-order model with a likelihood-ratio test. Because both models were fitted with MLR, we use lavTestLRT() rather than anova(). First, we refit the first-order two-factor CFA:

model_1st <- '
  PA  =~ pa1 + pa2 + pa3
  neg =~ na1 + na2 + na3
  PA  ~~ neg
'
fit_1st <- cfa(model_1st, data = ps, estimator = "MLR")
lavTestLRT(fit_1st, fit_2nd)
## 
## Chi-Squared Difference Test
## 
##         Df    AIC    BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_2nd  7 5248.3 5301.4  9.87                                    
## fit_1st  8 5246.3 5295.6  9.87                0       1

The output tests whether the second-order model fits significantly worse than the first-order model. A non-significant result means the second-order model fits as well as the first-order model — the more parsimonious second-order structure is supported. A significant result means the second-order restriction worsens fit, and the first-order model is preferred.

WarninglavTestLRT() for MLR models

When comparing nested models fitted with MLR, do not use anova() — it does not apply the Satorra–Bentler scaling correction. Use lavTestLRT() instead, which by default uses the scaled difference test. We cover this in more detail in Tutorial 13.

Drawing the second-order model

tidySEM::graph_sem(fit_2nd, std = TRUE)

Path diagram for the second-order factor model of affect. The top circle is the second-order factor ‘affect’; the middle circles are the first-order factors PA and neg; the bottom squares are the observed indicators.

The diagram shows three layers:

  1. Top: the second-order factor affect (circle)
  2. Middle: the two first-order factors PA and neg (circles), with arrows from affect — these are the second-order loadings
  3. Bottom: the six observed indicators (squares), with arrows from their respective first-order factors — these are the first-order loadings

When to use parcels and second-order factors

Technique When to use When to avoid
Item-level CFA When testing a measurement model is the primary goal When the model is too large to converge
Parcels When the item-level structure is confirmed and you need a simpler model for a larger SEM When you have not verified unidimensionality; when item-level information matters
Second-order factors When theory posits a hierarchical structure with distinct sub-factors under a common construct When the first-order factors are not theoretically related; when you have fewer than three first-order factors

Common errors

WarningParcels with NaN instead of NA

rowMeans(..., na.rm = TRUE) returns NaN when all items in a parcel are missing. lavaan treats NaN differently from NA, which can cause problems. After parceling, replace NaN with NA:

ps$p1[is.nan(ps$p1)] <- NA

Second-order model with two first-order factors is just-identified

If your second-order factor has only two first-order indicators, the second-order portion is just-identified. The overall model may still be overidentified through the first-order portion, but the second-order structure itself is not testable. Do not claim that the second-order model “fits well” based solely on the overall fit indices — the second-order portion could be replaced by any other just-identified relationship without changing the fit.

Negative residual variances in second-order models

Second-order models sometimes produce Heywood cases — negative residual variances for the first-order factors. This often means the second-order structure is not appropriate for the data; the first-order factors are not sufficiently related to share a common higher-order cause.

Exercises

  1. Parcel the six job-control items into three parcels of two items each (use the pairing from the tutorial). Fit a one-factor CFA on the parcels. Then fit a one-factor CFA on the six individual items. Compare the standardised loadings — are the parcel loadings higher? Why?

  2. Create a second-order factor model where affect is measured by PA and neg, and PA and neg are each measured by their three items. Compare this model’s fit indices to the first-order two-factor CFA from Tutorial 07. Does the second-order restriction worsen fit?

  3. Use lavTestLRT() to formally compare the first-order and second-order models from Exercise 2. Is the second-order model significantly worse than the first-order model?

  4. A researcher wants to parcel 12 items into 4 parcels of 3 items each. They plan to fit a one-factor CFA on the 4 parcels. How many degrees of freedom will this model have? Is it overidentified?

Exercise 1

ps$p1 <- (ps$con9 + ps$con10) / 2
ps$p2 <- (ps$con1 + ps$con8)  / 2
ps$p3 <- (ps$con3 + ps$con6)  / 2

fit_parcel <- cfa("con =~ p1 + p2 + p3", data = ps, estimator = "MLR")
fit_item   <- cfa("con =~ con1 + con3 + con6 + con8 + con9 + con10",
                  data = ps, estimator = "MLR")

summary(fit_parcel, standardized = TRUE)
## lavaan 0.6-21 ended normally after 17 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         6
## 
##   Number of observations                           343
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                                 0.000       0.000
##   Degrees of freedom                                 0           0
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   con =~                                                                
##     p1                1.000                               0.690    0.773
##     p2                0.803    0.066   12.254    0.000    0.554    0.797
##     p3                0.902    0.072   12.602    0.000    0.622    0.795
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .p1                0.320    0.038    8.443    0.000    0.320    0.402
##    .p2                0.177    0.021    8.599    0.000    0.177    0.366
##    .p3                0.226    0.029    7.863    0.000    0.226    0.369
##     con               0.475    0.061    7.734    0.000    1.000    1.000
summary(fit_item, standardized = TRUE)
## lavaan 0.6-21 ended normally after 24 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        12
## 
##   Number of observations                           343
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                                16.447      15.854
##   Degrees of freedom                                 9           9
##   P-value (Chi-square)                           0.058       0.070
##   Scaling correction factor                                  1.037
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   con =~                                                                
##     con1              1.000                               0.518    0.644
##     con3              1.238    0.136    9.128    0.000    0.641    0.697
##     con6              1.209    0.123    9.848    0.000    0.626    0.690
##     con8              1.104    0.104   10.595    0.000    0.572    0.691
##     con9              1.348    0.131   10.328    0.000    0.698    0.731
##     con10             1.293    0.139    9.326    0.000    0.670    0.576
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .con1              0.379    0.038   10.073    0.000    0.379    0.585
##    .con3              0.435    0.042   10.384    0.000    0.435    0.514
##    .con6              0.430    0.045    9.575    0.000    0.430    0.523
##    .con8              0.358    0.030   11.832    0.000    0.358    0.522
##    .con9              0.426    0.043    9.844    0.000    0.426    0.466
##    .con10             0.905    0.068   13.345    0.000    0.905    0.669
##     con               0.268    0.041    6.488    0.000    1.000    1.000

The parcel loadings should be higher than most of the individual-item loadings because averaging two items reduces measurement error, making the parcel a more reliable indicator of the factor.

Exercise 2

model_2nd <- '
  PA     =~ pa1 + pa2 + pa3
  neg    =~ na1 + na2 + na3
  affect =~ PA + neg
'
fit_2nd <- cfa(model_2nd, data = ps, estimator = "MLR")
fitMeasures(fit_2nd, c("chisq", "df", "pvalue", "cfi", "tli",
                       "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "srmr"))
##          chisq             df         pvalue            cfi            tli 
##          9.870          7.000          0.196          0.997          0.994 
##          rmsea rmsea.ci.lower rmsea.ci.upper           srmr 
##          0.035          0.000          0.082          0.015

Compare these to the Tutorial 07 fit indices (\(\chi^2(8) = 8.79\), CFI = 1.00, RMSEA = 0.02, SRMR = 0.01). The second-order model should have slightly more degrees of freedom and very similar fit if the second-order structure is appropriate.

Exercise 3

fit_1st <- cfa("PA =~ pa1 + pa2 + pa3
               neg =~ na1 + na2 + na3
               PA ~~ neg", data = ps, estimator = "MLR")
lavTestLRT(fit_1st, fit_2nd)
## 
## Chi-Squared Difference Test
## 
##         Df    AIC    BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_2nd  7 5248.3 5301.4  9.87                                    
## fit_1st  8 5246.3 5295.6  9.87                0       1

A non-significant result means the second-order model fits as well as the first-order model — the more parsimonious second-order structure is supported.

Exercise 4

With \(q = 4\) parcels, the covariance matrix has \(4 \times 5 / 2 = 10\) elements. Free parameters: 3 loadings (first fixed to 1) + 1 latent variance + 4 residual variances = 8. So \(\text{df} = 10 - 8 = 2\). The model is overidentified — it can be tested against the data.

Summary & next tutorial

You can now construct parcels from individual items and fit CFA models on parceled indicators, and you can specify and fit second-order factor models in lavaan. Both techniques appear frequently in published organisational psychology research, particularly in SEM models where the measurement model is not the primary focus. In the next tutorial we move from measurement models to structural models: we fit a path analysis with observed variables — a system of regression equations — and learn why a model with zero degrees of freedom fits perfectly.

Further reading

  • Beaujean (2014) — covers parceling and second-order factors in detail, with lavaan code.
  • Kelloway (2015) — the textbook basis for the original Mplus teaching materials; discusses when to parcel and when not to.
  • Hoyle (2012) — the Handbook of Structural Equation Modeling includes chapters on parceling and hierarchical factor models.

References

Beaujean, A. A. (2014). Latent variable modeling using r: A step-by-step guide. Routledge.
Hoyle, R. H. (2012). Handbook of structural equation modeling.
Kelloway, E. K. (2015). Using mplus for structural equation modeling (2nd ed.). SAGE Publications.