Define constant to use as multiplier in equality constraint (then estimate it)

[JB]

Dear all,
 
I am currently working in quantitative genetics with lavaan and everything has run great so far! In this field there are theoretical constraints, e.g. genes are equal in identical twins but, in average, only half the genes are the same in non-identical twins. We thus use 1* and 0.5* multipliers a lot. In one model however, I need to estimate this constant instead of fixing it. Below, I use the HS.model as a mock example to try and do sth similar but I am unsure of my syntax.


The first step: use a integer (I chose 2 because it was roughly the ratio between the two loadings I labeled l1, l2).

HS.model1 <- '
visual =~ x1 + l1*x2 + x3
textual =~ x4 + l2*x5 + x6
speed =~ x7 + x8 + x9
l2==2*l1
'
fit1 <- cfa(HS.model1, data=HolzingerSwineford1939)
summary(fit1)

The second step: define a multiplier

HS.model2 <- '
visual =~ x1 + l1*x2 + x3
textual =~ x4 + l2*x5 + x6
speed =~ x7 + x8 + x9
multi:=2
l2==multi*l1
'
fit2 <- cfa(HS.model2, data=HolzingerSwineford1939)
summary(fit2)

This one is strictly equivalent to the previous one. And I have no clue on how to "free" the label "multi" to estimate the multiplier instead of fixing it to 2.

Third step: I tried with a phantom variable approach that seem to work. 

HS.model3 <- '
visual =~ x1 + l1*x2 + x3
textual =~ x4 + l2*x5 + x6
speed =~ x7 + x8 + x9
phantom =~ 0   # create a phantom latent variable
phantom ~~ 0*phantom # fix the residual variance of f to zero
phantom~multi*1+NA*1#fit the interecept to 1
l2==multi*l1
'
fit3 <- cfa(HS.model3, data=HolzingerSwineford1939)
summary(fit3) #you get an estimate of multi, slightly superior to 2

The question is thus two-fold: is the last model correct and is there a more elegant way to do that.

Thanks for lavaan and for your help!!!

JB

[Alex Schoemann]

So it sounds like you want to estimate the ratio between two parameter estimates. To get an estimate of multi from the equation l2==multi*l1, just solve for multi so the new syntax should read:

HS.model2 <- '

visual =~ x1 + l1*x2 + x3

textual =~ x4 + l2*x5 + x6

speed =~ x7 + x8 + x9 

#multi:=2

multi:=l2/l1

'

This will create the new parameter multi and estimate it (with the same value as from the phantom variable model.

Alex

[JB]

Hi Alex,
Thanks for your quick answer!!
My real model is slightly more complicated that the one I posted (sorry for that). In fact, the multiplier is part of an equality constraint where parameters in one group are defined by parameters in another group. Normally, the multiplier is fixed to 0.5 and in one model it has to be estimated. So here is an example closer to my real model.

In the first one, I use the phantom approach. As you will see, two parameters in the second group (b3 and b5) are defined by a path using two parameters in the first group and the multiplier. I get a multiplier of 2.061. When I check the equalities, this multiplier seems correctly estimated.

HS.model <- 'visual =~ c(a1,b1)*x1 + c(a2,b2)*x2 + c(a3,b3)*x3
textual =~ c(a4,b4)*x4 + c(a5,b5)*x5 + c(a6,b6)*x6
speed =~ c(a7,b7)*x7 + c(a8,b8)*x8 + c(a9,b9)*x9

phantom =~ 0   # create a phantom latent variable
phantom ~~ 0*phantom # fix the residual variance of f to zero

phantom~multi*1+NA*1#label the multiplier and free it
b3==a3*a5*multi
b5==a2*a6*multi
'
fit <- cfa(HS.model, data = HolzingerSwineford1939, group = "school")
summary(fit)
 
Now, I use the define operator and define the multiplier based on other parameters and this time the multiplier is 2.873 and, in my example, is way off the other estimation.

HS.model1 <- 'visual =~ c(a1,b1)*x1 + c(a2,b2)*x2 + c(a3,b3)*x3
textual =~ c(a4,b4)*x4 + c(a5,b5)*x5 + c(a6,b6)*x6
speed =~ c(a7,b7)*x7 + c(a8,b8)*x8 + c(a9,b9)*x9
multi:=b3/(a3*a5)
multi:=b5/(a2*a6)
'
fit1 <- cfa(HS.model1, data = HolzingerSwineford1939, group = "school")
summary(fit1)

To my understanding the two models are not the same. In the first, I don't ask the model to estimate b3 and b5 because they are estimated based on parameters in other groups (which is what I need to do). In the second approach, b3 and b5 without constraint and the multi is simply calculated from the estimated parameters. Two questions: Any insight on why the models are different? Is there a more simple way to fit the first model (thus including the inequality constraint).
Many Thanks!!!
JB

[Alex Schoemann]

Interesting. I would have expected those two models to be identical. I did a bit of playing around and it looks like if you change the last 3 lines of fit1 to this:

multi:=b3/(a3*a5)

multi1:=b5/(a2*a6)

multi == multi1

Then the models are identical (same estimates, same fit). I have no idea why your second model (fit1) doesn't give the same information (it has the same number of estimated parameters). My best guess is that it involves the estimator being used but that's above my pay grade. 

Hope this helps,

Alex

[JB]

Hi Alex,
Thanks!! It does yield the same results in the example. When I adapt your multi, multi1 one syntaxe to my data, I get an error message:

  In nlminb(start = x.par, objective = auglag, gradient = fgrad, control = control, : NA/NaN function evaluation

So it may be that the phantom and the define approaches come to the same results in the Hs.model by different routes? And one of them does not work in my model because it is quite complex... Maybe the optimizer is not the same as the constraints are different. But it is also above my pay grade :) I think I'll try to estimate this model in OpenMx and see if I can get results close to the ones I get with the phantom method, the only one that works for my data. If you have new insights, please post them, your comments were very helpful!!
JB

[Alex Schoemann]

I'm guessing you've run into a starting value issue with the constrained estimator. I found this from a discussion with Yves in October (the full discussion is here: https://groups.google.com/forum/?hl=en#!searchin/lavaan/chime/lavaan/Zm8_LNC91Co/DbNRlSWp7rcJ):

" (for the current version of lavaan, 0.5-14): as 
soon as "==" or "<" or ">" statements are added to the model syntax, 
lavaan 'switches' to the constrained optimizer (nlminb.constr), which is 
very, very sensitive to starting values. "

The good news is that there's a quick and easy way to get good starting values. Just estimate the model without the equality constrain in place (e.g. in the above model take out the mult == mult1 line). Then use the results of that model as starting values for the constrained model (using the start option, which can take a lavaan object as an argument).

Alex