Hi all,
I am running a simple univariate saturated threshold model to find out the effects of covariates on the variable of interest across three groups - MZ twins, DZ twins and unrelated individuals. The regressions are the following: y ~age + sex + IBS.
As one of the research questions, I would like to specify such a constraint for IBS coefficients: beta_unrel = 2*beta_DZ. This is where I got stuck. I looked at the following thread https://openmx.ssri.psu.edu/thread/446 and wanted to specify an mxAlgebra instead of mxConstraint, but I struggle to modify my original model to test this assumption.
This is the chunk of the code for the regressions from the original model:
meanMZ <- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMeanMZ" ) meanDZ <- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMeanDZ" ) meanUnrel <- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="expMeanUnrel" ) threMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=.8, labels='threMZ', name="threMZ" ) threDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=.8, labels='threDZ', name="threDZ" ) threUnrel <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=.8, labels='threUnrel', name="threUnrel" ) defAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.Zage_s_1", name="Age1" ) defAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.Zage_s_2", name="Age2" ) B_AgeMZ <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaAgeMZ", name="bAgeMZ" ) B_AgeDZ <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaAgeDZ", name="bAgeDZ" ) B_AgeUnrel <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaAgeUnrel", name="bAgeUnrel" ) threAge1MZ <- mxAlgebra( Age1%*%bAgeMZ, name="AgeR1MZ") threAge2MZ <- mxAlgebra( Age2%*%bAgeMZ, name="AgeR2MZ") threAge1DZ <- mxAlgebra( Age1%*%bAgeDZ, name="AgeR1DZ") threAge2DZ <- mxAlgebra( Age2%*%bAgeDZ, name="AgeR2DZ") threAge1Unrel <- mxAlgebra( Age1%*%bAgeUnrel, name="AgeR1Unrel") threAge2Unrel <- mxAlgebra( Age2%*%bAgeUnrel, name="AgeR2Unrel") defSex1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.sex_s_1", name="Sex1" ) defSex2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.sex_s_2", name="Sex2" ) B_SexMZ <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaSexMZ", name="bSexMZ" ) B_SexDZ <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaSexDZ", name="bSexDZ" ) B_SexUnrel <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaSexUnrel", name="bSexUnrel" ) threSex1MZ <- mxAlgebra( Sex1%*%bSexMZ, name="SexR1MZ") threSex2MZ <- mxAlgebra( Sex2%*%bSexMZ, name="SexR2MZ") threSex1DZ <- mxAlgebra( Sex1%*%bSexDZ, name="SexR1DZ") threSex2DZ <- mxAlgebra( Sex2%*%bSexDZ, name="SexR2DZ") threSex1Unrel <- mxAlgebra( Sex1%*%bSexUnrel, name="SexR1Unrel") threSex2Unrel <- mxAlgebra( Sex2%*%bSexUnrel, name="SexR2Unrel") defIBS1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.ibs_gen_s_1", name="IBS1" ) defIBS2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels="data.ibs_gen_s_2", name="IBS2" ) B_IBSMZ <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaIBSMZ", name="bIBSMZ" ) B_IBSDZ <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaIBSDZ", name="bIBSDZ" ) B_IBSUnrel <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaIBSUnrel", name="bIBSUnrel" ) threIBS1MZ <- mxAlgebra( IBS1%*%bIBSMZ, name="IBSR1MZ") threIBS2MZ <- mxAlgebra( IBS2%*%bIBSMZ, name="IBSR2MZ") threIBS1DZ <- mxAlgebra( IBS1%*%bIBSDZ, name="IBSR1DZ") threIBS2DZ <- mxAlgebra( IBS2%*%bIBSDZ, name="IBSR2DZ") threIBS1Unrel <- mxAlgebra( IBS1%*%bIBSUnrel, name="IBSR1Unrel") threIBS2Unrel <- mxAlgebra( IBS2%*%bIBSUnrel, name="IBSR2Unrel") expThreMZ <- mxAlgebra( expression= threMZ + cbind(AgeR1MZ,AgeR2MZ) + cbind(SexR1MZ,SexR2MZ) + cbind(IBSR1MZ,IBSR2MZ), name="expThreMZ" ) expThreDZ <- mxAlgebra( expression= threDZ + cbind(AgeR1DZ,AgeR2DZ) + cbind(SexR1DZ,SexR2DZ) + cbind(IBSR1DZ,IBSR2DZ), name="expThreDZ" ) expThreUnrel <- mxAlgebra( expression= threUnrel + cbind(AgeR1Unrel,AgeR2Unrel) + cbind(SexR1Unrel,SexR2Unrel) + cbind(IBSR1Unrel,IBSR2Unrel), name="expThreUnrel" )
All this, together with other entities, go into three subgroups and into the total model:
modelMZ <- mxModel( name="MZ", meanMZ, covMZ, threMZ, defAge1, defAge2, B_AgeMZ, threAge1MZ, threAge2MZ, defSex1, defSex2, B_SexMZ, threSex1MZ, threSex2MZ, defIBS1, defIBS2, B_IBSMZ, threIBS1MZ, threIBS2MZ, expThreMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", meanDZ, covDZ, threDZ, defAge1, defAge2, B_AgeDZ, threAge1DZ, threAge2DZ, defSex1, defSex2, B_SexDZ, threSex1DZ, threSex2DZ, defIBS1, defIBS2, B_IBSDZ, threIBS1DZ, threIBS2DZ, expThreDZ, dataDZ, expDZ, funML ) modelUnrel <- mxModel( name="Unrel", meanUnrel, covUnrel, threUnrel, defAge1, defAge2, B_AgeUnrel, threAge1Unrel, threAge2Unrel, defSex1, defSex2, B_SexUnrel, threSex1Unrel, threSex2Unrel, defIBS1, defIBS2, B_IBSUnrel, threIBS1Unrel, threIBS2Unrel, expThreUnrel, dataUnrel, expUnrel, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ","Unrel") ) univSatModel <- mxModel("univSat", modelMZ, modelDZ, modelUnrel, multi) univSatFit <- mxRun(univSatModel)
Now I would like to test whether we can impose betaIBSUnrel = 2*betaIBSDZ. I tried the following:
restrModel = mxModel(model1, name='restrM') restrModel = omxSetParameters(restrModel , labels='betaIBSDZ', values=0.01, free=T, newlabels = 'betaIBSDZ') B_IBSDZ_Unrel <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.01, label="betaIBSDZ", name="bIBSDZ" ) B_IBSUnrel_matrix <- mxMatrix( type="Full", nrow=1, ncol=nv, values=.01, label="bIBSUnrel[1,1]", name="bIBSUnrel_matrix" ) B_IBSUnrel_alg = mxAlgebra(2*betaIBSDZ, name='bIBSUnrel') restrModel = mxModel(GenIBSModel, B_IBSDZ_Unrel, B_IBSUnrel_matrix, B_IBSUnrel_alg) restrFit = mxRun(restrModel)
But this doesn't seem to do the right thing. The new 'bIBSUnrel' from B_IBSUnrel_alg is calculated along the old 'bIBSUnrel' in the Unrel subgroup. This is the summary of the model fit:
> summary(restrFit) Summary of restrM free parameters: name matrix row col Estimate Std.Error A 1 betaIBSDZ bIBSDZ 1 1 -0.13358435 0.13235843 2 covT1T2MZ MZ.expCovMZ strain_dich_s_1 strain_dich_s_2 0.42856769 0.11593230 3 threshAll MZ.threMZ 1 1 0.64756220 0.10746847 4 betaAgeMZ MZ.bAgeMZ 1 1 0.07123071 0.07590719 5 betaSexMZ MZ.bSexMZ 1 1 -0.14514120 0.14911503 6 covT1T2DZ DZ.expCovDZ strain_dich_s_1 strain_dich_s_2 0.22069379 0.13846191 7 betaAgeDZ DZ.bAgeDZ 1 1 0.37860942 0.09224387 8 betaSexDZ DZ.bSexDZ 1 1 0.07751815 0.14411199 9 threshUnrel Unrel.threUnrel 1 1 0.96161404 0.03937058 10 betaAgeUnrel Unrel.bAgeUnrel 1 1 0.10529660 0.02691893 11 betaSexUnrel Unrel.bSexUnrel 1 1 -0.13111884 0.05280346 12 betaIBSUnrel Unrel.bIBSUnrel 1 1 -0.49658313 0.09702530 > mxEval(c(betaIBSDZ, bIBSUnrel), restrFit) [,1] [1,] -0.1335843 [2,] -0.2671687
You can see that there is still a matrix bIBSUnrel in Unrel group which element was estimated freely (-0.497), whereas the new bIBSUnrel matrix that I created for the submodel is equal to two times betaIBSDZ.
Could you please help me to code this across group constraint and to run a restricted model with it?
Thank you in advance!
Julia
Your syntax doesn't work because
B_IBSUnrel_matrix
is going into the container model, not the submodelmodelUnrel
, so the submodel's expectation never "sees" it. Plus, your syntax never fixes the parameter labeled "bIBSUnrel", which is still there in that submodel. I think the easiest thing is probably to just define the regression coefficient for IBS among unrelateds as an algebra,, or perhaps,
, then redefine
modelUnrel
to use it instead ofB_IBSUnrel
,and then rebuild the multigroup model using the redefined
modelUnrel
.Hi Rob. Thanks a bunch!!! This does seem to work! Thank you again!
Glad to hear it!