Interaction effect between categorical variable and latent variable

[Kss2k]

👋 Welcome!

We’re using Discussions as a place to connect with other members of our community. We hope that you:

• Ask questions you’re wondering about.
• Share ideas.
• Engage with other community members.
• Welcome others and are open-minded. Remember that this is a community we
  build together 💪.

To get started, comment below with an introduction of yourself and tell us about what you do with this community.

[claboeck]

Hello! I am new to R, trying to use it to estimate SEM with two latent IVs (both from ordinal indicators), a binary observed IV for a moderator, and two observed ordinal DVs.

The latents are structural and sociocultural. I am trying to model (1) the path of structural (latent) --> acad_satisfaction, (2) sociocultural (latent) --> belonging, (3) Pell as a moderator for paths (1) and (2).

I see that in Lavaan the : is only for interactions between 2 observed variables, and it was not clear how to interact an observed and latent with indProd.

Here is what I tried with modsem:
try1 <- '
structural=~ strcl_1 + strcl_2 + strcl_3 + strcl_4
sociocultural=~soccl_1 + soccl_2 + soccl_4 + soccl_5 + soccl_6 + soccl_7 + soccl_8

acad_s ~ structural + pell + structural:pell
belonging ~ sociocultural + pell + sociocultural:pell
'
est1 <- modsem(try1, oneInt, method = "pind", std.lv=TRUE)
summary(est1, standardized=TRUE)

I receive the error message "Unable to find observed variables in data: [1] "strcl_1" "strcl_2" "strcl_3" "strcl_4" "soccl_1" "soccl_2" "soccl_4"
[8] "soccl_5" "soccl_6" "soccl_7" "soccl_8"

This seems like an easy syntax fix...what am I doing wrong?

I also want to know how (or if) I need to specify that the indicators for the latents (strcl_1-4, soccl_1-8) are categorical like with Lavaan, and if/how to specify WLMSV.

Thank you!

[Kss2k]

Hey!

In this case you need to replace oneInt with your own dataset. oneInt is a dataset included in the modsem package.
However, in this case, I wouldn't use modsem at all. Since "pell" is a binary variable you should probably use a multigroup
model instead. From experience the product indicator approaches don't work that well with nominal variables. Here is an example using a generated dataset, which you can replace with your own. Using a multigroup model, "pell" is of course
treated as a categorical variable. Wheter or not you should treat your indicators as categorical is up to you. Depending on the number of categories, and their distribution (and other assumptions), it is often okay to treat them as continuous, and using ordered variables with lavaan can be a bit tricky if not all categories are present in both groups. Regardless I have included two examples, one treating the indicators as continuous, and the other as categegorical.

library(mvtnorm)
library(lavaan)

# function for ordinalizing a continuous variable (return as integer)
ordinalize <- function(x) {
  ifelse(abs(x) < 2.5, round(x, 0), sign(x) * 2) # really simple
}

ordinalize_data <- function(x) apply(x, MARGIN = 2, ordinalize)

# function for generating indicators
gen_indicators <- function(lV, k, prefix, sigma = sqrt(2), 
                           lambda = 0.9, rand = 0.1) {
  theta <- diag(runif(k, min = sigma - rand, max = sigma + rand))
  factorloadings <- c(1, runif(k - 1, min = lambda - rand, 
                               max = lambda - rand))
  residuals <- as.matrix(rmvnorm(N, mean = rep(0, k), sigma = theta))
  indicators <- lV %*% t(factorloadings) + residuals
  colnames(indicators) <- paste0(prefix, seq_len(k))
  indicators
}


# generating the dataset
set.seed(122)
N <- 1000
phi <- diag(1, 2) # covariance matrix of exogenous LVs
phi[2, 1] <- phi[1, 2] <- 0.3
exogenous <- as.data.frame(rmvnorm(N, mean = c(0, 0), sigma = phi))
structural <- exogenous[[1]]
sociocultural <- exogenous[[2]]
pell <- sample(c(0, 1), N, replace = TRUE)
# ordinal indicators (note is integer, and treated as continuous in the model)
strcl <- ordinalize_data(gen_indicators(structural, k = 4, "strcl_"))
soccl <- ordinalize_data(gen_indicators(sociocultural, k = 8, "soccl_"))

acad_s <- 0.2 * structural + # simple main effect
  0.3 * pell + # simple main effect
  0.1 * pell * structural + # interaction
  rnorm(N, sd = 0.5) # residual
belonging <- 0.1* sociocultural + # simple main effect
  0.6 * pell + # simple main effect
  0.3 * pell * sociocultural + # interaction
  rnorm(N, sd = 0.4) # residual

my_dataset <- cbind(data.frame(pell, acad_s, belonging), strcl, soccl)


# estimating the model in lavaan
model <- '
# messurement/outer model
  structural =~ strcl_1 + strcl_2 + strcl_3 + strcl_4
  sociocultural =~ soccl_1 + soccl_2 + soccl_3 + soccl_4 + 
    soccl_5 + soccl_6 + soccl_7 + soccl_8

# structural/inner model
  acad_s ~ c(ia_1, ia_2) * 1 + c(st_1, st_2) * structural
  belonging ~ c(ib_1, ib_2) * 1 + c(sc_1, sc_2) * sociocultural

# interaction effects and simple main effects for pell
  structural_pell := st_2 - st_1 # acad_s ~ structural:pell
  sociocultural_pell := sc_2 - sc_1 # belonging  ~ sociocultural:pell 
  main_pell_acad_s := ia_2 - ia_1 # acad_s ~ pell
  main_pell_belonging := ib_2 - ib_1 # belonging ~ pell
'

# treating indicators as continuous
est_cont <- sem(model, data = my_dataset, group = "pell")
summary(est_cont) 

# treating indicators as ordinal 
ordinal_vars <- colnames(my_dataset)[grepl("strcl_|soccl_", colnames(my_dataset))]
est_ord <- sem(model, data = my_dataset, group = "pell", ordered = ordinal_vars)
summary(est_ord)

[claboeck]

Thank you so much for your help! May I please ask some follow-up questions?

1. when I asked for fit statistics (summary(est_ord, fit.measures = TRUE, rsquare = TRUE, standardized = TRUE), I do not get fit statistics that makes sense. I ran a SEM without the interaction effect, and I would expect the RMSEA, CFI, TLI to be similar. I am not sure what exactly is going on with the code you shared, but is there something in there that would change how a RMSEA is calculated?

2. Is the code comparable to how a person would run the model without the interaction term with the WLSMV estimator? I ask because the factor loadings for both groups look different (and lower!) than the factor loadings for my SEM without the interaction term. I specified "WLSMV" in the model command but it did not seem to change anything (est_ord <- sem(model, data = sem_pell_data, group = "pell", ordered = ordinal_vars, std.lv=TRUE, estimator = "WLSMV"))

3. what do the 0.1, 0.3, and 0.6 numbers multiplying Pell represent in these lines:

acad_s <- 0.2 * structural + # simple main effect
0.3 * pell + # simple main effect
0.1 * pell * structural + # interaction
rnorm(N, sd = 0.5) # residual
belonging <- 0.1* sociocultural + # simple main effect
0.6 * pell + # simple main effect
0.3 * pell * sociocultural + # interaction
rnorm(N, sd = 0.4) # residual

3. at the bottom of the output, I see sociocltrl_pll, structural_pll, main_pell_cd_s, and man_pll_belonging. are the first two the interaction effects (sociocultrl_pll = the interaction of sociocultural*pell--> belonging?)

I attached my output if that is helpful. Sorry to bother you - there is a lot to learn! I am used to Stata and Mplus.
SEM_Clssism_Pell output.docx

[Kss2k]

1. For the fit measures, it would depend on how you ran the model without the interaction? But in general the interaction effects for the simulated data are quite large, and you would thus expect a relatively big difference between fit measures for H0 and H1. I also see that TLI and CFI are slightly above 1 (which should be interpreted as 1), and that the RMSEA is < .001. The fact that the fit-measures are so good makes sense, given that the data is generated from the model (i.e., we have the true model for explaining our observed data).
2. Im not quite sure what you mean here, but you should be able to run a model without an interaction term yes. The way I would do that would be something like this.

# null model 
null_model <- '
# messurement/outer model
  structural =~ strcl_1 + strcl_2 + strcl_3 + strcl_4
  sociocultural =~ soccl_1 + soccl_2 + soccl_3 + soccl_4 + 
    soccl_5 + soccl_6 + soccl_7 + soccl_8

# structural/inner model
  acad_s ~ c(ia_1, ia_2) * 1 + c(st_1, st_2) * structural
  belonging ~ c(ib_1, ib_2) * 1 + c(sc_1, sc_2) * sociocultural

# simple main effects for pell
  main_pell_acad_s := ia_2 - ia_1 # acad_s ~ pell
  main_pell_belonging := ib_2 - ib_1 # belonging ~ pell
# interaction effect equal accross groups
  sc_1 == sc_2 
  st_1 == st_2
'

est_cont_null <- sem(null_model, data = my_dataset, group = "pell")
summary(est_cont_null, fit.measures = TRUE)

3. These numbers represent the population coefficients. e.g,. acad_s <- 0.2 * structural + # simple main effect is the simple main effect of structural upon acad_s, where the coefficient is 0.2.
4. The first two are the interaction effects yes, whilst the last two the main effects for pell (hint: look at the comments behind #). NOTE: The first two lines were named incorrectly, i.e., the first line should have been named structural_pell not sociocultural_pell, and the same for the second line. This has been corrected now.
5. I realize that estimating the interaction effects using a multigroup model, might be a 'bit' confusing, if you are new to both R and lavaan. You could also consider using the QML approach, it does not support categorical variables per se, but it is quite robust to normality violations. You can see that it performs quite well for the simulated data (i.e., the the coefficients are close to the population parameters used to simulate the data). NOTE if you use the QML approach it will correspond to a regression model where both of the variables in the interaction term have been centered before hand. It also only works with the latest GitHub version.

model_qml <- '
# messurement/outer model
  STRUCTURAL =~ strcl_1 + strcl_2 + 
    strcl_3 + strcl_4
  SOCIOCULTURAL =~ soccl_1 + soccl_2 + soccl_4 + 
    soccl_5 + soccl_6 + soccl_7 + soccl_8
  ACAD_S =~ acad_s 
  BELONGING =~ belonging
  PELL =~ pell

# structural/inner model
  ACAD_S ~ STRUCTURAL + PELL + STRUCTURAL:PELL
  BELONGING ~ SOCIOCULTURAL + PELL + SOCIOCULTURAL:PELL 
'

est_qml <- modsem(model_qml, my_dataset, method = "qml")
summary(est_qml)