RI-CLPM with Auxiliary Variables in lavaan: Theta Warning and Saturated Correlates Approach

[Ronja Steinhauser]

Hello everyone, 

I am currently estimating a random-intercept cross-lagged panel model (RI-CLPM) in lavaan and would like to incorporate auxiliary variables (dummy-coded disciplinary background variables) to strengthen the plausibility of the Missing at Random (MAR) assumption.

Using semTools::cfa.auxiliary(), the model converges without any problems. However, I receive the following warning:

lavaan WARNING: the covariance matrix of the residuals of the observed variables (theta) is not positive definite.

The smallest eigenvalue of the theta matrix was −0.0326.  

Despite this warning, the model appears to be functioning. The overall model fit is good, there are no Heywood cases or negative residual variances, and the parameter estimates are virtually identical to those obtained from the model without auxiliary variables.

While searching for an explanation, I found a previous discussion (https://groups.google.com/g/lavaan/c/kr0JWacZLT8/m/rfGquwzSCgAJ). Based on the suggestions in that discussion, I attempted to implement the auxiliary variables manually using the saturated correlates approach, specifying them directly in the lavaan syntax ( specifying covariances among each of the observed covariates, covariances between the observed covariates and all within person factors, and variances for the observed covariates)

However, this approach raises a new set of issues:

• The model's degrees of freedom change, which (according to the previous discussion:  https://groups.google.com/g/lavaan/c/kr0JWacZLT8/m/rfGquwzSCgAJ ) suggests that I have not specified the saturated correlates structure correctly.
• I also receive a different warning indicating that the parameter variance-covariance matrix (vcov) is not positive definite, with the smallest eigenvalue being 6.23 × 10⁻¹⁵.

At this point, I am unsure whether the problem lies in my manual specification or whether the warning produced by cfa.auxiliary() might be expected in this type of model.

I would therefore greatly appreciate any advice on the following questions:

1. Is semTools::cfa.auxiliary() intended to be used with RI-CLPMs, or is it preferable to manually implement the saturated correlates approach manually?
2. If cfa.auxiliary() is appropriate, should the Theta not positive definite warning be considered problematic in this context, even though the model fit and parameter estimates appear reasonable?
3. If manual implementation is the recommended approach, is there anything obviously missing or incorrectly specified in my syntax that could explain the changed degrees of freedom?

I have attached my model syntax (with cfa.auxiliary and manual) and would be very grateful for any suggestions or comments. 

Thank you very much for your time and help!

Approach using cfa.auxiliary:

RICLPM_lifesat <- '

  #####################
  # MEASUREMENT MODEL #
  #####################

# Factor models for perceived success at 4 waves (constrained)
  psuc_t1 =~ 1*psuc1_t1 + a*psuc2_t1 + b*psuc3_t1
  psuc_t2 =~ 1*psuc1_t2 + a*psuc2_t2 + b*psuc3_t2
  psuc_t3 =~ 1*psuc1_t3 + a*psuc2_t3 + b*psuc3_t3
  psuc_t4 =~ 1*psuc1_t4 + a*psuc2_t4 + b*psuc3_t4
 
 
# Constrained intercepts over time (this is necessary for strong factorial
# invariance; without these contraints we have weak factorial invariance)
  psuc1_t1 + psuc1_t2 + psuc1_t3 + psuc1_t4 ~ i1*1
  psuc2_t1 + psuc2_t2 + psuc2_t3 + psuc2_t4 ~ i2*1
  psuc3_t1 + psuc3_t2 + psuc3_t3 + psuc3_t4 ~ i3*1

# Free latent means from t = 2 onward (only do this in combination with
# constraints on intercepts; without these, this would not be identified)
  psuc_t2 + psuc_t3 + psuc_t4 ~ 1

# Correlated residuals of identical indicators across waves
  psuc1_t1 ~~ psuc1_t2 + psuc1_t3 + psuc1_t4
  psuc1_t2 ~~ psuc1_t3 + psuc1_t4
  psuc1_t3 ~~ psuc1_t4

  psuc2_t1 ~~ psuc2_t2 + psuc2_t3 + psuc2_t4
  psuc2_t2 ~~ psuc2_t3 + psuc2_t4
  psuc2_t3 ~~ psuc2_t4

  psuc3_t1 ~~ psuc3_t2 + psuc3_t3 + psuc3_t4
  psuc3_t2 ~~ psuc3_t3 + psuc3_t4
  psuc3_t3 ~~ psuc3_t4


  ################
  # BETWEEN PART #
  ################

# Create between factors (random intercepts)
  RI_psuc =~ 1*psuc_t1 + 1*psuc_t2 + 1*psuc_t3 + 1*psuc_t4
  RI_lsat =~ 1*lsat_t1 + 1*lsat_t2 + 1*lsat_t3 + 1*lsat_t4


# Set residual variances of repeated measures (RI_psuc, RI_lsat) to zero
  psuc_t1 ~~ 0*psuc_t1
  psuc_t2 ~~ 0*psuc_t2
  psuc_t3 ~~ 0*psuc_t3
  psuc_t4 ~~ 0*psuc_t4

  lsat_t1 ~~ 0*lsat_t1
  lsat_t2 ~~ 0*lsat_t2
  lsat_t3 ~~ 0*lsat_t3
  lsat_t4 ~~ 0*lsat_t4


  ###############
  # WITHIN PART #
  ###############

# Create within-components
  W_psuc1 =~ 1*psuc_t1
  W_psuc2 =~ 1*psuc_t2
  W_psuc3 =~ 1*psuc_t3
  W_psuc4 =~ 1*psuc_t4

  W_lsat1 =~ 1*lsat_t1
  W_lsat2 =~ 1*lsat_t2
  W_lsat3 =~ 1*lsat_t3
  W_lsat4 =~ 1*lsat_t4

# Specify lagged effects between within-person centered latent variables
  W_psuc2 + W_lsat2 ~ W_psuc1 + W_lsat1
  W_psuc3 + W_lsat3 ~ W_psuc2 + W_lsat2
  W_psuc4 + W_lsat4 ~ W_psuc3 + W_lsat3

# Estimate correlations within same wave
  W_psuc1 ~~ W_lsat1
  W_psuc2 ~~ W_lsat2
  W_psuc3 ~~ W_lsat3
  W_psuc4 ~~ W_lsat4


  ##########################
  # ADDITIONAL CONSTRAINTS #
  ##########################

  # Set correlations between between-factors (random intercepts) and within-
  # factors at wave 1 at 0
 
  RI_psuc + RI_lsat ~~ 0*W_psuc1 + 0*W_lsat1
 
'

fit_RICLPM_lifesat <- cfa.auxiliary(
  RICLPM_lifesat,
  data = data,
  aux = c("fg_Andere", "fg_Ingenieur", "fg_MatNat", "fg_Medizin", "fg_Sozial"), #Reference group Humanities
  missing = "FIML",
  estimator = "MLR",
  fixed.x = FALSE
)


summary(
  fit_RICLPM_lifesat,
  standardized = TRUE,
  fit.measures = TRUE,
  rsquare = TRUE
)

Manual approach:

RICLPM_lifesat <- '

  #####################
  # MEASUREMENT MODEL #
  #####################

# Factor models for perceived success at 4 waves (constrained)
  psuc_t1 =~ 1*psuc1_t1 + a*psuc2_t1 + b*psuc3_t1
  psuc_t2 =~ 1*psuc1_t2 + a*psuc2_t2 + b*psuc3_t2
  psuc_t3 =~ 1*psuc1_t3 + a*psuc2_t3 + b*psuc3_t3
  psuc_t4 =~ 1*psuc1_t4 + a*psuc2_t4 + b*psuc3_t4
 
 
# Constrained intercepts over time (this is necessary for strong factorial
# invariance; without these contraints we have weak factorial invariance)
  psuc1_t1 + psuc1_t2 + psuc1_t3 + psuc1_t4 ~ i1*1
  psuc2_t1 + psuc2_t2 + psuc2_t3 + psuc2_t4 ~ i2*1
  psuc3_t1 + psuc3_t2 + psuc3_t3 + psuc3_t4 ~ i3*1

# Free latent means from t = 2 onward (only do this in combination with
# constraints on intercepts; without these, this would not be identified)
  psuc_t2 + psuc_t3 + psuc_t4 ~ 1

# Correlated residuals of identical indicators across waves
  psuc1_t1 ~~ psuc1_t2 + psuc1_t3 + psuc1_t4
  psuc1_t2 ~~ psuc1_t3 + psuc1_t4
  psuc1_t3 ~~ psuc1_t4

  psuc2_t1 ~~ psuc2_t2 + psuc2_t3 + psuc2_t4
  psuc2_t2 ~~ psuc2_t3 + psuc2_t4
  psuc2_t3 ~~ psuc2_t4

  psuc3_t1 ~~ psuc3_t2 + psuc3_t3 + psuc3_t4
  psuc3_t2 ~~ psuc3_t3 + psuc3_t4
  psuc3_t3 ~~ psuc3_t4


  ################
  # BETWEEN PART #
  ################

# Create between factors (random intercepts)
  RI_psuc =~ 1*psuc_t1 + 1*psuc_t2 + 1*psuc_t3 + 1*psuc_t4
  RI_lsat =~ 1*lsat_t1 + 1*lsat_t2 + 1*lsat_t3 + 1*lsat_t4


# Set residual variances of repeated measures (RI_psuc, RI_lsat) to zero
  psuc_t1 ~~ 0*psuc_t1
  psuc_t2 ~~ 0*psuc_t2
  psuc_t3 ~~ 0*psuc_t3
  psuc_t4 ~~ 0*psuc_t4

  lsat_t1 ~~ 0*lsat_t1
  lsat_t2 ~~ 0*lsat_t2
  lsat_t3 ~~ 0*lsat_t3
  lsat_t4 ~~ 0*lsat_t4


  ###############
  # WITHIN PART #
  ###############

# Create within-components
  W_psuc1 =~ 1*psuc_t1
  W_psuc2 =~ 1*psuc_t2
  W_psuc3 =~ 1*psuc_t3
  W_psuc4 =~ 1*psuc_t4

  W_lsat1 =~ 1*lsat_t1
  W_lsat2 =~ 1*lsat_t2
  W_lsat3 =~ 1*lsat_t3
  W_lsat4 =~ 1*lsat_t4

# Specify lagged effects between within-person centered latent variables
  W_psuc2 + W_lsat2 ~ W_psuc1 + W_lsat1
  W_psuc3 + W_lsat3 ~ W_psuc2 + W_lsat2
  W_psuc4 + W_lsat4 ~ W_psuc3 + W_lsat3

# Estimate correlations within same wave
  W_psuc1 ~~ W_lsat1
  W_psuc2 ~~ W_lsat2
  W_psuc3 ~~ W_lsat3
  W_psuc4 ~~ W_lsat4


  ##########################
  # ADDITIONAL CONSTRAINTS #
  ##########################

  # Set correlations between between-factors (random intercepts) and within-
  # factors at wave 1 at 0
 
  RI_psuc + RI_lsat ~~ 0*W_psuc1 + 0*W_lsat1
 
 
 
   #######################
  # AUXILIARY VARIABLES #
  #######################

  # variances of auxiliary variables
  #humanities is reference category
  fg_Andere    ~~ fg_Andere
  fg_Ingenieur ~~ fg_Ingenieur
  fg_MatNat    ~~ fg_MatNat
  fg_Medizin   ~~ fg_Medizin
  fg_Sozial    ~~ fg_Sozial

  # covariances among auxiliary variables
  fg_Andere ~~ fg_Ingenieur + fg_MatNat + fg_Medizin + fg_Sozial
  fg_Ingenieur ~~ fg_MatNat + fg_Medizin + fg_Sozial
  fg_MatNat ~~ fg_Medizin + fg_Sozial
  fg_Medizin ~~ fg_Sozial

  # saturated correlations with latent within-person factors
  W_psuc1 + W_psuc2 + W_psuc3 + W_psuc4 ~~
    fg_Andere + fg_Ingenieur + fg_MatNat + fg_Medizin + fg_Sozial

  W_lsat1 + W_lsat2 + W_lsat3 + W_lsat4 ~~
    fg_Andere + fg_Ingenieur + fg_MatNat + fg_Medizin + fg_Sozial
'

fit_RICLPM_lifesat <- cfa(
  RICLPM_lifesat,
  data = data,
  missing = "FIML",
  estimator = "MLR",
  fixed.x = FALSE
)


summary(
  fit_RICLPM_lifesat,
  standardized = TRUE,
  fit.measures = TRUE,
  rsquare = TRUE
)

[Terrence Jorgensen]

Dear Ronja,

It looks like your manual script is correlating the auxiliaries with common factors instead of their indicators. If you pass each model to lavInspect(), the default output will show you which elements of each parameter matrix are being estimated, so you can compare to see how they differ. 

FYI, the latest development version of lavaan now implements this approach with its own aux= argument. Could you try installing that to compare it with semTools?

remotes::install_github("yrosseel/lavaan")

Best,

Terrence

[Ronja]

Dear Terrence,

Thank you very much for your helpful response.

First, I tested the new aux= argument in the latest development version of lavaan.

 However, lavaan removed all of my auxiliary variables because they are binary dummy variables and issued the following warning:

“auxiliary (aux) variable(s) are binary or ordered/categorical; this is not supported and they are removed.”

Therefore, I was unfortunately unable to compare the new lavaan implementation directly with semTools::cfa.auxiliary().

I also followed your suggestion and compared the two specifications using lavInspect().

This confirmed your suspicion: with cfa.auxiliary(), the auxiliary variables are correlated with the observed variables—that is, with the indicators of perceived success (psuc1–psuc3 at each wave) and with the observed life-satisfaction (lsat) variables. In my manual specification, by contrast, I had correlated the auxiliary variables with the latent within-person factors.

I had based this manual specification on examples of RI-CLPMs with manifest repeated measures, in which the auxiliary variables were correlated with the corresponding single-indicator within-person factors. My model, however, contains one construct measured by multiple indicators and one manifest construct. 

Is there a preferred approach in this situation? If cfa.auxiliary() is the recommended approach, should the Theta warning be interpreted as indicating a remaining issue?

Thank you again for your time,

Ronja

[Terrence Jorgensen]

I had based this manual specification on examples of RI-CLPMs with manifest repeated measures, in which the auxiliary variables were correlated with the corresponding single-indicator within-person factors. My model, however, contains one construct measured by multiple indicators and one manifest construct. 

Is there a preferred approach in this situation?

 

It is not the saturated correlates approach unless your df are unaffected, which is what requires saturating all of your observed variables with the auxiliaries. 

If cfa.auxiliary() is the recommended approach, should the Theta warning be interpreted as indicating a remaining issue?

I think that is unrelated.  Dummy codes follow a Bernoulli distribution, which is binomial with n=1 trial.  In that case, the mean and variance are based on the same parameter p (the probability of "success"): https://en.wikipedia.org/wiki/Bernoulli_distribution

That means their sampling distributions are perfectly correlated, which leads to the error about the asymptotic sample-covariance matrix (i.e., vcov() output) being NPD.  You could prevent the warning by labeling both parameters and specifying the constraint explicitly in the syntax (for each dummy code):

dummyVar ~ intercept*1

dummyVar ~~ variance*dummyVar

variance == intercept * (1 - intercept)

But the goal of the saturated correlates approach is to leave your df unaffected, so I think you can just ignore it.  But you could check that there are no other causes of an NPD matrix by implementing the constraints to check there is no further warning.

Best,

Terrence D. Jorgensen    (he, him, his)
Assistant Professor, Methods and Statistics
Research Institute for Child Development and Education, the University of Amsterdam
http://www.uva.nl/profile/t.d.jorgensen