Using a MML estimator in a CFA

[Jeffrey Henry]

Dear all,

I would like to know relevant syntax specifications and/or R packages I should download to produce a marginal maximum likelihood (MML) estimator in a confirmatory factor analysis (CFA). My CFA syntax (with a MLR estimator) is ready. Could you please point me to some relevant references, if any? I looked up online with no luck.

Thanks in advance,

Jeffrey Henry

[yrosseel]

On 03/31/2014 01:17 PM, Jeffrey Henry wrote:
> Dear all,
>
> I would like to know relevant syntax specifications and/or R packages I
> should download to produce a marginal maximum likelihood (MML) estimator
> in a confirmatory factor analysis (CFA).

I'm not sure what you are asking for. Syntax? Nothing special, just
describe your model as usual. To use MML, simply add estimator="MML" to
the cfa/sem/lavaan call.

For example:

# create binary data
HS9 <-
HolzingerSwineford1939[,c("x1","x2","x3","x4","x5","x6","x7","x8","x9")]
HSbinary <- as.data.frame( lapply(HS9, cut, 2, labels=FALSE) )

# describe model (here a one-dimensional cfa)
model <- ' trait =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 '

# fit model
fit <- cfa(model, data=HSbinary, ordered=names(HSbinary), std.lv=TRUE,
estimator="MML")

# look at results
summary(fit)


Other R packages? Look at the 'mirt' package. Which is currently much,
much faster (using C++), but MML only.

References? Try anything with 'IRT'.

Some books:

van der Linden, W. J., & Hambleton, R. K. (Eds.). (1997). Handbook of
modern item response theory. New York, NY: Springer.

Hambleton, R.K., and Swaminathan, H. Item Response Theory: Principles
and Applications. Hingham, MA: Kluwer, Nijhoff, 1984.

Skrondal, A., & Rabe-Hesketh, S. (2004). Generalized latent variable
modeling: Multilevel, longitudinal, and structural equation models. CRC
Press.

Bartholomew, D. J., Knott, M., & Moustaki, I. (2011). Latent variable
models and factor analysis: a unified approach (Vol. 899). John Wiley &
Sons.

Yves.

[Jeffrey Henry]

Dear Mr Rosseel,

Thanks for this quick response. I had already tried using the command estimator="MML" in my model. However, I get the following error:

Error in fml_deriv1(Sigma.hat = Sigma.hat[[g]], TH = TH[[g]], th.idx = th.idx[[g]],  : 
  not implemented

N.B.: Everything works fine when I use a ML or a MLR estimator.

Maybe the kind of model I am testing requires further adjustments. I am testing a three-bifactor model where all items (24) load on a general factor, but also on three, uncorrelated subfactors (i.e., each subfactor includes loadings for a specific subset of items). I am testing this model in a twin sample. Therefore, I included all relevant specifications for interchangeable dyads (see the Olsen and Kenny paper attached) + MZ/DZ zygosity groups. It is possible that further adjustments should be made considering the specific nature of this model? Or maybe the MML estimator is not suitable for this kind of analysis?

Thanks again for your help,

Jeffrey Henry

[yrosseel]

On 04/01/2014 01:46 PM, Jeffrey Henry wrote:
> Dear Mr Rosseel,
>
> Thanks for this quick response. I had already tried using the command
> estimator="MML" in my model. However, I get the following error:
>
> Error in fml_deriv1(Sigma.hat = Sigma.hat[[g]], TH = TH[[g]], th.idx = th.idx[[g]], :
> not implemented

You need to update lavaan (or maybe even R) so that you have the current
CRAN version: 0.5-16

Yves.

[Jeffrey Henry]

Thanks, I just did (switched from R 3.0.2 to R 3.0.3, and to lavaan 05.-16). However, I get a slightly different error message now:

Error in lav_partable_flat(FLAT, meanstructure = meanstructure, int.ov.free = int.ov.free,  : 
  lavaan ERROR: parameterization mml not supported

It seems there is a problem with me freeing the intercepts of the observed variables, but I get the same error message when I try to fix them to 0 or 1. Would you like me to send you my syntax? It may be easier this way.

Thanks again,

Jeffrey

[yrosseel]

On 04/01/2014 02:56 PM, Jeffrey Henry wrote:
> It seems there is a problem with me freeing the intercepts of the
> observed variables, but I get the same error message when I try to fix
> them to 0 or 1. Would you like me to send you my syntax? It may be
> easier this way.

Yes.

Yves.

[Jeffrey Henry]

Thanks for taking the time. Here it is. It is a big syntax, but quite simple really. All items load on a general factor + their respective, uncorrelated subfactors. Each factor is represented twice. For example, CU1 = twin 1, and CU2 = twin 2. The same goes for the observed variables, assessed for both twins 1 and 2. There are also other minor specifications (e.g., equality constraints on observed variables' intercepts, equality constraints in factor loadings from family (twin 1, twin 2) and zygosity group (MZ, DZ)); see the Olsen and Kelly paper attached in my other message). Grouping is made according to zygosity (Group 1 = MZ; Group 2 = DZ).

biftwinzyg <- '
#Factor loadings
  CU1 =~ c(b46,b46)*ppbhicut021 + c(b1,b1)*ppbhicut041 + c(b2,b2)*ppbhicut071 + c(b3,b3)*ppbhicut08r1 + c(b4,b4)*ppbhicut09r1 + c(b5,b5)*ppbhicut101 + c(b6,b6)*ppbhicut111 + c(b7,b7)*ppbhicut121 + c(b8,b8)*ppbhicut181 + c(b9,b9)*ppbhicut201 + c(b10,b10)*ppbhicut211 + c(b11,b11)*ppbhicut03r1 + c(b12,b12)*ppbhicut05r1 + c(b13,b13)*ppbhicut13r1 + c(b14,b14)*ppbhicut15r1 + c(b15,b15)*ppbhicut16r1 + c(b16,b16)*ppbhicut17r1 + c(b17,b17)*ppbhicut23r1 + c(b18,b18)*ppbhicut24r1 + c(b19,b19)*ppbhicut01r1 + c(b20,b20)*ppbhicut061 + c(b21,b21)*ppbhicut071 + c(b22,b22)*ppbhicut14r1 + c(b23,b23)*ppbhicut19r1 + c(b24,b24)*ppbhicut221
  CU2 =~ c(b46,b46)*ppbhicut022 + c(b1,b1)*ppbhicut042 + c(b2,b2)*ppbhicut072 + c(b3,b3)*ppbhicut08r2 + c(b4,b4)*ppbhicut09r2 + c(b5,b5)*ppbhicut102 + c(b6,b6)*ppbhicut112 + c(b7,b7)*ppbhicut122 + c(b8,b8)*ppbhicut182 + c(b9,b9)*ppbhicut202 + c(b10,b10)*ppbhicut212 + c(b11,b11)*ppbhicut03r2 + c(b12,b12)*ppbhicut05r2 + c(b13,b13)*ppbhicut13r2 + c(b14,b14)*ppbhicut15r2 + c(b15,b15)*ppbhicut16r2 + c(b16,b16)*ppbhicut17r2 + c(b17,b17)*ppbhicut23r2 + c(b18,b18)*ppbhicut24r2 + c(b19,b19)*ppbhicut01r2 + c(b20,b20)*ppbhicut062 + c(b21,b21)*ppbhicut072 + c(b22,b22)*ppbhicut14r2 + c(b23,b23)*ppbhicut19r2 + c(b24,b24)*ppbhicut222
  callousness1 =~ c(b47,b47)*ppbhicut021 + c(b25,b25)*ppbhicut041 + c(b26,b26)*ppbhicut071 + c(b27,b27)*ppbhicut08r1 + c(b28,b28)*ppbhicut09r1 + c(b29,b29)*ppbhicut101 + c(b30,b30)*ppbhicut111 + c(b31,b31)*ppbhicut121 + c(b32,b32)*ppbhicut181 + c(b33,b33)*ppbhicut201 + c(b34,b34)*ppbhicut211
  callousness2 =~ c(b47,b47)*ppbhicut022 + c(b25,b25)*ppbhicut042 + c(b26,b26)*ppbhicut072 + c(b27,b27)*ppbhicut08r2 + c(b28,b28)*ppbhicut09r2 + c(b29,b29)*ppbhicut102 + c(b30,b30)*ppbhicut112 + c(b31,b31)*ppbhicut122 + c(b32,b32)*ppbhicut182 + c(b33,b33)*ppbhicut202 + c(b34,b34)*ppbhicut212
  uncaring1 =~ c(b48,b48)*ppbhicut03r1 + c(b35,b35)*ppbhicut05r1 + c(b36,b36)*ppbhicut13r1 + c(b37,b37)*ppbhicut15r1 + c(b38,b38)*ppbhicut16r1 + c(b39,b39)*ppbhicut17r1 + c(b40,b40)*ppbhicut23r1 + c(b41,b41)*ppbhicut24r1
  uncaring2 =~ c(b48,b48)*ppbhicut03r2 + c(b35,b35)*ppbhicut05r2 + c(b36,b36)*ppbhicut13r2 + c(b37,b37)*ppbhicut15r2 + c(b38,b38)*ppbhicut16r2 + c(b39,b39)*ppbhicut17r2 + c(b40,b40)*ppbhicut23r2 + c(b41,b41)*ppbhicut24r2
  unemotional1 =~ c(b49,b49)*ppbhicut01r1 + c(b42,b42)*ppbhicut061 + c(b43,b43)*ppbhicut14r1 + c(b44,b44)*ppbhicut19r1 + c(b45,b45)*ppbhicut221
  unemotional2 =~ c(b49,b49)*ppbhicut01r2 + c(b42,b42)*ppbhicut062 + c(b43,b43)*ppbhicut14r2 + c(b44,b44)*ppbhicut19r2 + c(b45,b45)*ppbhicut222
#Residual variances observed variables
  ppbhicut021 ~~ c(v5,v5)*ppbhicut021
  ppbhicut022 ~~ c(v5,v5)*ppbhicut022
  ppbhicut041 ~~ c(v6,v6)*ppbhicut041
  ppbhicut042 ~~ c(v6,v6)*ppbhicut042
  ppbhicut071 ~~ c(v7,v7)*ppbhicut071
  ppbhicut072 ~~ c(v7,v7)*ppbhicut072
  ppbhicut08r1 ~~ c(v8,v8)*ppbhicut08r1
  ppbhicut08r2 ~~ c(v8,v8)*ppbhicut08r2
  ppbhicut09r1 ~~ c(v9,v9)*ppbhicut09r1
  ppbhicut09r2 ~~ c(v9,v9)*ppbhicut09r2
  ppbhicut101 ~~ c(v10,v10)*ppbhicut101
  ppbhicut102 ~~ c(v10,v10)*ppbhicut102
  ppbhicut111 ~~ c(v11,v11)*ppbhicut111
  ppbhicut112 ~~ c(v11,v11)*ppbhicut112
  ppbhicut121 ~~ c(v12,v12)*ppbhicut121
  ppbhicut122 ~~ c(v12,v12)*ppbhicut122
  ppbhicut181 ~~ c(v13,v13)*ppbhicut181
  ppbhicut182 ~~ c(v13,v13)*ppbhicut182
  ppbhicut201 ~~ c(v14,v14)*ppbhicut201
  ppbhicut202 ~~ c(v14,v14)*ppbhicut202
  ppbhicut211 ~~ c(v15,v15)*ppbhicut211
  ppbhicut212 ~~ c(v15,v15)*ppbhicut212
  ppbhicut03r1 ~~ c(v16,v16)*ppbhicut03r1
  ppbhicut03r2 ~~ c(v16,v16)*ppbhicut03r2
  ppbhicut05r1 ~~ c(v17,v17)*ppbhicut05r1
  ppbhicut05r2 ~~ c(v17,v17)*ppbhicut05r2
  ppbhicut13r1 ~~ c(v18,v18)*ppbhicut13r1
  ppbhicut13r2 ~~ c(v18,v18)*ppbhicut13r2
  ppbhicut15r1 ~~ c(v19,v19)*ppbhicut15r1
  ppbhicut15r2 ~~ c(v19,v19)*ppbhicut15r2
  ppbhicut16r1 ~~ c(v20,v20)*ppbhicut16r1
  ppbhicut16r2 ~~ c(v20,v20)*ppbhicut16r2
  ppbhicut17r1 ~~ c(v21,v21)*ppbhicut17r1
  ppbhicut17r2 ~~ c(v21,v21)*ppbhicut17r2
  ppbhicut23r1 ~~ c(v22,v22)*ppbhicut23r1
  ppbhicut23r2 ~~ c(v22,v22)*ppbhicut23r2
  ppbhicut24r1 ~~ c(v23,v23)*ppbhicut24r1
  ppbhicut24r2 ~~ c(v23,v23)*ppbhicut24r2
  ppbhicut01r1 ~~ c(v24,v24)*ppbhicut01r1
  ppbhicut01r2 ~~ c(v24,v24)*ppbhicut01r2
  ppbhicut061 ~~ c(v25,v25)*ppbhicut061
  ppbhicut062 ~~ c(v25,v25)*ppbhicut062
  ppbhicut14r1 ~~ c(v26,v26)*ppbhicut14r1
  ppbhicut14r2 ~~ c(v26,v26)*ppbhicut14r2
  ppbhicut19r1 ~~ c(v27,v27)*ppbhicut19r1
  ppbhicut19r2 ~~ c(v27,v27)*ppbhicut19r2
  ppbhicut221 ~~ c(v28,v28)*ppbhicut221
  ppbhicut222 ~~ c(v28,v28)*ppbhicut222
#Factor variance(s)
  CU1 ~~ 1*CU1
  callousness1 ~~ 1*callousness1
  uncaring1 ~~ 1*uncaring1
  unemotional1 ~~ 1*unemotional1
  CU2 ~~ 1*CU2
  callousness2 ~~ 1*callousness2
  uncaring2 ~~ 1*uncaring2
  unemotional2 ~~ 1*unemotional2
#Factor covariance(s)
  CU1 ~~ CU2
  callousness1 ~~ callousness2
  uncaring1 ~~ uncaring2
  unemotional1 ~~ unemotional2
#Residual covariances observed variables
  ppbhicut021 ~~ ppbhicut022
  ppbhicut041 ~~ ppbhicut042
  ppbhicut071 ~~ ppbhicut072
  ppbhicut08r1 ~~ ppbhicut08r2
  ppbhicut09r1 ~~ ppbhicut09r2
  ppbhicut101 ~~ ppbhicut102
  ppbhicut111 ~~ ppbhicut112
  ppbhicut121 ~~ ppbhicut122
  ppbhicut181 ~~ ppbhicut182
  ppbhicut201 ~~ ppbhicut202
  ppbhicut211 ~~ ppbhicut212
  ppbhicut01r1 ~~ ppbhicut01r2
  ppbhicut061 ~~ ppbhicut062
  ppbhicut14r1 ~~ ppbhicut14r2
  ppbhicut19r1 ~~ ppbhicut19r2
  ppbhicut221 ~~ ppbhicut222
  ppbhicut03r1 ~~ ppbhicut03r2
  ppbhicut05r1 ~~ ppbhicut05r2
  ppbhicut13r1 ~~ ppbhicut13r2
  ppbhicut15r1 ~~ ppbhicut15r2
  ppbhicut16r1 ~~ ppbhicut16r2
  ppbhicut17r1 ~~ ppbhicut17r2
  ppbhicut23r1 ~~ ppbhicut23r2
  ppbhicut24r1 ~~ ppbhicut24r2
#Intercepts obvserved variables
ppbhicut021 ~ c(i1,i1)*1
ppbhicut022 ~ c(i1,i1)*1
ppbhicut041 ~ c(i2,i2)*1
ppbhicut042 ~ c(i2,i2)*1
ppbhicut071 ~ c(i3,i3)*1
ppbhicut072 ~ c(i3,i3)*1
ppbhicut08r1 ~ c(i4,i4)*1
ppbhicut08r2 ~ c(i4,i4)*1
ppbhicut09r1 ~ c(i5,i5)*1
ppbhicut09r2 ~ c(i5,i5)*1
ppbhicut101 ~ c(i6,i6)*1
ppbhicut102 ~ c(i6,i6)*1
ppbhicut111 ~ c(i7,i7)*1
ppbhicut112 ~ c(i7,i7)*1
ppbhicut121 ~ c(i8,i8)*1
ppbhicut122 ~ c(i8,i8)*1
ppbhicut181 ~ c(i9,i9)*1
ppbhicut182 ~ c(i9,i9)*1
ppbhicut201 ~ c(i10,i10)*1
ppbhicut202 ~ c(i10,i10)*1
ppbhicut211 ~ c(i11,i11)*1
ppbhicut212 ~ c(i11,i11)*1
ppbhicut03r1 ~ c(i12,i12)*1
ppbhicut03r2 ~ c(i12,i12)*1
ppbhicut05r1 ~ c(i13,i13)*1
ppbhicut05r2 ~ c(i13,i13)*1
ppbhicut13r1 ~ c(i14,i14)*1
ppbhicut13r2 ~ c(i14,i14)*1
ppbhicut15r1 ~ c(i15,i15)*1
ppbhicut15r2 ~ c(i15,i15)*1
ppbhicut16r1 ~ c(i16,i16)*1
ppbhicut16r2 ~ c(i16,i16)*1
ppbhicut17r1 ~ c(i17,i17)*1
ppbhicut17r2 ~ c(i17,i17)*1
ppbhicut23r1 ~ c(i18,i18)*1
ppbhicut23r2 ~ c(i18,i18)*1
ppbhicut24r1 ~ c(i19,i19)*1
ppbhicut24r2 ~ c(i19,i19)*1
ppbhicut01r1 ~ c(i20,i20)*1
ppbhicut01r2 ~ c(i20,i20)*1
ppbhicut061 ~ c(i21,i21)*1
ppbhicut062 ~ c(i21,i21)*1
ppbhicut14r1 ~ c(i22,i22)*1
ppbhicut14r2 ~ c(i22,i22)*1
ppbhicut19r1 ~ c(i23,i23)*1
ppbhicut19r2 ~ c(i23,i23)*1
ppbhicut221 ~ c(i24,i24)*1
ppbhicut222 ~ c(i24,i24)*1'
biftwinzyg.fit <- lavaan(biftwinzyg, data = TEDSr, group="pzygos",std.lv=TRUE,estimator="MML")


Jeffrey

[Jeffrey Henry]

Mr. Rosseel,

Did you get my syntax alright (its in my last post, from last tuesday)?

Thanks in advance,

Jeffrey Henry

[yrosseel]

On 04/01/2014 03:33 PM, Jeffrey Henry wrote:
> Thanks for taking the time. Here it is. It is a big syntax, but quite
> simple really. All items load on a general factor + their respective,
> uncorrelated subfactors. Each factor is represented twice. For example,
> CU1 = twin 1, and CU2 = twin 2. The same goes for the observed
> variables, assessed for both twins 1 and 2. There are also other minor
> specifications (e.g., equality constraints on observed variables'
> intercepts, equality constraints in factor loadings from family (twin 1,
> twin 2) and zygosity group (MZ, DZ)); see the Olsen and Kelly paper
> attached in my other message). Grouping is made according to zygosity
> (Group 1 = MZ; Group 2 = DZ).

Are the indicators categorical (binary/ordered)?

If yes, so need to specify the thresholds in the syntax (or use the
sem() function instead of the lavaan() function if you are unsure which
parameters are needed; or use the auto.th=TRUE flag in the lavaan()
function). Also, categorical indicators do not have 'free' residual
variances (at least not under the default 'delta' parameterization), so
you can not impose equality constraints on them. You might consider
parameterization="theta" instead. Intercepts are also not needed, since
you have thresholds (which you may constrain to be equal across groups).

If no, and all your data is continuous, please use estimator="ML".
Estimator MML currently (0.5-16) only works for categorical data. And
for continuous data, MML == ML.

Yves.