You are here

Sex-Limitation model in two different cohorts

1 post / 0 new
JuanJMV's picture
Joined: 07/20/2016 - 13:13
Sex-Limitation model in two different cohorts

Dear all,

I am trying to fit a sex-limitation model (5 groups) but in two different cohorts (1944-1980) and (1981-2000). So, I have 10 groups (MZM, MZF, DZM, DZF and DZO x 2 cohorts).

My variable is ordinal with 3 categories (2 thresholds). I know that for model identification, I have to: (1) fix the mean (or regression intercept) and variance, and free all thresholds; (2) free the mean (or regression intercept) and fix the variance and one of the thresholds; or (3) free the mean (or regression intercept) and the variance, and fix two of the thresholds.

I would like to compare raw variances (VA, VC, VE and VT) among cohorts, so if I fix the variance to one in the four groups, I cannot do it. I have then fixed the two thresholds and free mean and variance. However, I am not 100% sure this is a valid approach in a model with two different cohorts. I would really appreciate if you could guide me with this.

I have also seen this paper ( that says that “having more than 2 categories, threshold information can be used to estimate the mean and variance of the liability rather than fixing them” but I am not sure how to do it.

What would be the most suitable approach?

Here the script that I have adapted as reference:

nth <- 2 # number of thresholds
# Select Variables for Analysis
vars <- 'St' # list of variables names
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
# Select Data for Analysis
mzfDataONE <- subset(ordData, Zyg==3 & Cohort==0, c(selVars))
dzfDataONE <- subset(ordData, Zyg==4 & Cohort==0, c(selVars))
mzmDataONE <- subset(ordData, Zyg==1 & Cohort==0, c(selVars))
dzmDataONE <- subset(ordData, Zyg==2 & Cohort==0, c(selVars))
dzoDataONE <- subset(ordData, Zyg==5 & Cohort==0, c(selVars))
mzfDataTWO <- subset(ordData, Zyg==3 & Cohort==1, c(selVars))
dzfDataTWO <- subset(ordData, Zyg==4 & Cohort==1, c(selVars))
mzmDataTWO <- subset(ordData, Zyg==1 & Cohort==1, c(selVars))
dzmDataTWO <- subset(ordData, Zyg==2 & Cohort==1, c(selVars))
dzoDataTWO <- subset(ordData, Zyg==5 & Cohort==1, c(selVars))
mzfDataONE$St1 <- mxFactor(mzfDataONE$St1, levels =c(0,1,2))
mzfDataONE$St2 <- mxFactor(mzfDataONE$St2, levels =c(0,1,2))
dzfDataONE$St1 <- mxFactor(dzfDataONE$St1, levels =c(0,1,2))
dzfDataONE$St2 <- mxFactor(dzfDataONE$St2, levels =c(0,1,2))
mzmDataONE$St1 <- mxFactor(mzmDataONE$St1, levels =c(0,1,2))
mzmDataONE$St2 <- mxFactor(mzmDataONE$St2, levels =c(0,1,2))
dzmDataONE$St1 <- mxFactor(dzmDataONE$St1, levels =c(0,1,2))
dzmDataONE$St2 <- mxFactor(dzmDataONE$St2, levels =c(0,1,2))
dzoDataONE$St1 <- mxFactor(dzoDataONE$St31, levels =c(0,1,2))
dzoDataONE$St2 <- mxFactor(dzoDataONE$St2, levels =c(0,1,2))
mzfDataTWO$St1 <- mxFactor(mzfDataTWO$St1, levels =c(0,1,2))
mzfDataTWO$St2 <- mxFactor(mzfDataTWO$St2, levels =c(0,1,2))
dzfDataTWO$St1 <- mxFactor(dzfDataTWO$St1, levels =c(0,1,2))
dzfDataTWO$St2 <- mxFactor(dzfDataTWO$St2, levels =c(0,1,2))
mzmDataTWO$St1 <- mxFactor(mzmDataTWO$St1, levels =c(0,1,2))
mzmDataTWO$St2 <- mxFactor(mzmDataTWO$St2, levels =c(0,1,2))
dzmDataTWO$St1 <- mxFactor(dzmDataTWO$St1, levels =c(0,1,2))
dzmDataTWO$St2 <- mxFactor(dzmDataTWO$St2, levels =c(0,1,2))
dzoDataTWO$St1 <- mxFactor(dzoDataTWO$St1, levels =c(0,1,2))
dzoDataTWO$St2 <- mxFactor(dzoDataTWO$St2, levels =c(0,1,2))
# Set Starting Values
svBe <- 0.01 # start value for regressions
svLTh <- -1.5 # start value for first threshold
svITh <- 1 # start value for increments
svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds
lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds
svPa <- .4 # start value for path coefficient
svPe <- .8 # start value for path coefficient for e
# ACE Model
# Create Algebra for expected Mean Matrices
meanfONE    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1,label="mfONE", name="MeanfONE" )
meanmONE    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1,label="mmONE", name="MeanmONE" )
meanoONE    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1,label="mmONE", name="MeanoONE" )
meanfTWO    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1,label="mfTWO", name="MeanfTWO" )
meanmTWO    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1,label="mmTWO", name="MeanmTWO" )
meanoTWO    <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1,label="moTWO", name="MeanoTWO" )
expMeanZfONE <- mxAlgebra( expression= MeanfONE , name="expMeanZfONE" )
expMeanZmONE <- mxAlgebra( expression= MeanmONE , name="expMeanZmONE" )
expMeanZoONE <- mxAlgebra( expression= MeanoONE , name="expMeanZoONE" )
expMeanZfTWO <- mxAlgebra( expression= MeanfTWO , name="expMeanZfTWO" )
expMeanZmTWO <- mxAlgebra( expression= MeanmTWO , name="expMeanZmTWO" )
expMeanZoTWO <- mxAlgebra( expression= MeanoTWO , name="expMeanZoTWO" )
Inc      <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="Inc" )
Thre    <- mxMatrix( type="Full", nrow=nth, ncol=nv, free=c(F,F), 
                     name="Thre" )
threshold    <- mxAlgebra( expression= cbind(Inc %*% Thre,Inc %*% Thre), name="threshold" )
covAfONE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VAf11ONE", name="VAfONE" )
covCfONE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VCf11ONE", name="VCfONE" )
covEfONE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="VEf11ONE", name="VEfONE" )
covAmONE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VAm11ONE", name="VAmONE" )
covCmONE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VCm11ONE", name="VCmONE" )
covEmONE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="VEm11ONE", name="VEmONE" )
covAmsONE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=0, label="VAms11ONE", name="VAmsONE" )
covCmsONE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=FALSE, values=0, label="VCms11ONE", name="VCmsONE" )
covAfTWO <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VAf11TWO", name="VAfTWO" )
covCfTWO <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VCf11TWO", name="VCfTWO" )
covEfTWO <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="VEf11TWO", name="VEfTWO" )
covAmTWO <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VAm11TWO", name="VAmTWO" )
covCmTWO <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VCm11TWO", name="VCmTWO" )
covEmTWO <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="VEm11TWO", name="VEmTWO" )
covAmsTWO <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=0, label="VAms11TWO", name="VAmsTWO" )
covCmsTWO <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=FALSE, values=0, label="VCms11TWO", name="VCmsTWO" )
# Produce a vector which is = to 1 if variance is positive and -1 if variance is negative
signAONE <- mxAlgebra( ((-1)^omxLessThan(VAfONE,0))*((-1)^omxLessThan(VAmONE,0)), name="signAONE")
signCONE <- mxAlgebra( ((-1)^omxLessThan(VCfONE,0))*((-1)^omxLessThan(VCmONE,0)), name="signCONE")
signATWO <- mxAlgebra( ((-1)^omxLessThan(VAfTWO,0))*((-1)^omxLessThan(VAmTWO,0)), name="signATWO")
signCTWO <- mxAlgebra( ((-1)^omxLessThan(VCfTWO,0))*((-1)^omxLessThan(VCmTWO,0)), name="signCTWO")
# Calculate absolute covariation between Males and Females and then un-absoluting the product
lower <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=FALSE, values=1, name="lower" )
covAosONE <- mxAlgebra( signAONE*(sqrt(abs(VAfONE%*%lower))*t(sqrt(abs(VAmONE%*%lower)))), name="VAosONE")
covCosONE <- mxAlgebra( signCONE*(sqrt(abs(VCfONE%*%lower))*t(sqrt(abs(VCmONE%*%lower)))), name="VCosONE")
covAosTWO <- mxAlgebra( signATWO*(sqrt(abs(VAfTWO%*%lower))*t(sqrt(abs(VAmTWO%*%lower)))), name="VAosTWO")
covCosTWO <- mxAlgebra( signCTWO*(sqrt(abs(VCfTWO%*%lower))*t(sqrt(abs(VCmTWO%*%lower)))), name="VCosTWO")
# Calculate rg/rc from reparameterized model
pathRgONE <- mxAlgebra( signAONE*(sqrt(abs(VAfONE%*%lower))*t(sqrt(abs(VAmONE%*%lower))))/sqrt(VAfONE*(VAmONE+VAmsONE)), name="rgONE")
pathRcONE <- mxAlgebra( signCONE*(sqrt(abs(VCfONE%*%lower))*t(sqrt(abs(VCmONE%*%lower))))/sqrt(VCfONE*(VCmONE+VCmsONE)), name="rcONE")
pathRgTWO <- mxAlgebra( signATWO*(sqrt(abs(VAfTWO%*%lower))*t(sqrt(abs(VAmTWO%*%lower))))/sqrt(VAfTWO*(VAmTWO+VAmsTWO)), name="rgTWO")
pathRcTWO <- mxAlgebra( signCTWO*(sqrt(abs(VCfTWO%*%lower))*t(sqrt(abs(VCmTWO%*%lower))))/sqrt(VCfTWO*(VCmTWO+VCmsTWO)), name="rcTWO")
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covPfONE <- mxAlgebra( expression= VAfONE+VCfONE+VEfONE, name="VfONE" )
covPmONE <- mxAlgebra( expression= VAmONE+VCmONE+VEmONE+VAmsONE+VCmsONE, name="VmONE" )
covMZfONE <- mxAlgebra( expression= VAfONE+VCfONE, name="cMZfONE" )
covDZfONE <- mxAlgebra( expression= 0.5%x%VAfONE+ VCfONE, name="cDZfONE" )
covMZmONE <- mxAlgebra( expression= VAmONE+VCmONE+VAmsONE+VCmsONE, name="cMZmONE" )
covDZmONE <- mxAlgebra( expression= 0.5%x%VAmONE+ VCmONE+ 0.5%x%VAmsONE+VCmsONE, name="cDZmONE" )
covDZoONE <- mxAlgebra( expression= 0.5%x%VAosONE+VCosONE, name="cDZoONE" )
expCovMZfONE <- mxAlgebra( expression= rbind( cbind(VfONE, cMZfONE), cbind(t(cMZfONE), VfONE)), name="expCovMZfONE" )
expCovDZfONE <- mxAlgebra( expression= rbind( cbind(VfONE, cDZfONE), cbind(t(cDZfONE), VfONE)), name="expCovDZfONE" )
expCovMZmONE <- mxAlgebra( expression= rbind( cbind(VmONE, cMZmONE), cbind(t(cMZmONE), VmONE)), name="expCovMZmONE" )
expCovDZmONE <- mxAlgebra( expression= rbind( cbind(VmONE, cDZmONE), cbind(t(cDZmONE), VmONE)), name="expCovDZmONE" )
expCovDZoONE <- mxAlgebra( expression= rbind( cbind(VfONE, cDZoONE), cbind(t(cDZoONE), VmONE)), name="expCovDZoONE" )
covPfTWO <- mxAlgebra( expression= VAfTWO+VCfTWO+VEfTWO, name="VfTWO" )
covPmTWO <- mxAlgebra( expression= VAmTWO+VCmTWO+VEmTWO+VAmsTWO+VCmsTWO, name="VmTWO" )
covMZfTWO <- mxAlgebra( expression= VAfTWO+VCfTWO, name="cMZfTWO" )
covDZfTWO <- mxAlgebra( expression= 0.5%x%VAfTWO+ VCfTWO, name="cDZfTWO" )
covMZmTWO <- mxAlgebra( expression= VAmTWO+VCmTWO+VAmsTWO+VCmsTWO, name="cMZmTWO" )
covDZmTWO <- mxAlgebra( expression= 0.5%x%VAmTWO+ VCmTWO+ 0.5%x%VAmsTWO+VCmsTWO, name="cDZmTWO" )
covDZoTWO <- mxAlgebra( expression= 0.5%x%VAosTWO+VCosTWO, name="cDZoTWO" )
expCovMZfTWO <- mxAlgebra( expression= rbind( cbind(VfTWO, cMZfTWO), cbind(t(cMZfTWO), VfTWO)), name="expCovMZfTWO" )
expCovDZfTWO <- mxAlgebra( expression= rbind( cbind(VfTWO, cDZfTWO), cbind(t(cDZfTWO), VfTWO)), name="expCovDZfTWO" )
expCovMZmTWO <- mxAlgebra( expression= rbind( cbind(VmTWO, cMZmTWO), cbind(t(cMZmTWO), VmTWO)), name="expCovMZmTWO" )
expCovDZmTWO <- mxAlgebra( expression= rbind( cbind(VmTWO, cDZmTWO), cbind(t(cDZmTWO), VmTWO)), name="expCovDZmTWO" )
expCovDZoTWO <- mxAlgebra( expression= rbind( cbind(VfTWO, cDZoTWO), cbind(t(cDZoTWO), VmTWO)), name="expCovDZoTWO" )
# Create Data Objects for Multiple Groups
dataMZfONE <- mxData( observed=mzfDataONE, type="raw" )
dataDZfONE <- mxData( observed=dzfDataONE, type="raw" )
dataMZmONE <- mxData( observed=mzmDataONE, type="raw" )
dataDZmONE <- mxData( observed=dzmDataONE, type="raw" )
dataDZoONE <- mxData( observed=dzoDataONE, type="raw" )
dataMZfTWO <- mxData( observed=mzfDataTWO, type="raw" )
dataDZfTWO <- mxData( observed=dzfDataTWO, type="raw" )
dataMZmTWO <- mxData( observed=mzmDataTWO, type="raw" )
dataDZmTWO <- mxData( observed=dzmDataTWO, type="raw" )
dataDZoTWO <- mxData( observed=dzoDataTWO, type="raw" )
# Create Expectation Objects for Multiple Groups
expMZfONE <- mxExpectationNormal( covariance="expCovMZfONE", means="expMeanZfONE", dimnames=selVars, thresholds="threshold" )
expDZfONE <- mxExpectationNormal( covariance="expCovDZfONE", means="expMeanZfONE", dimnames=selVars, thresholds="threshold" )
expMZmONE <- mxExpectationNormal( covariance="expCovMZmONE", means="expMeanZmONE", dimnames=selVars, thresholds="threshold" )
expDZmONE <- mxExpectationNormal( covariance="expCovDZmONE", means="expMeanZmONE", dimnames=selVars, thresholds="threshold" )
expDZoONE <- mxExpectationNormal( covariance="expCovDZoONE", means="expMeanZoONE", dimnames=selVars, thresholds="threshold" )
expMZfTWO <- mxExpectationNormal( covariance="expCovMZfTWO", means="expMeanZfTWO", dimnames=selVars, thresholds="threshold" )
expDZfTWO <- mxExpectationNormal( covariance="expCovDZfTWO", means="expMeanZfTWO", dimnames=selVars, thresholds="threshold" )
expMZmTWO <- mxExpectationNormal( covariance="expCovMZmTWO", means="expMeanZmTWO", dimnames=selVars, thresholds="threshold" )
expDZmTWO <- mxExpectationNormal( covariance="expCovDZmTWO", means="expMeanZmTWO", dimnames=selVars, thresholds="threshold" )
expDZoTWO <- mxExpectationNormal( covariance="expCovDZoTWO", means="expMeanZoTWO", dimnames=selVars, thresholds="threshold" )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
parsZfONE <- list(  meanfONE, Inc,threshold,Thre, covAfONE, covCfONE, covEfONE, covPfONE )
parsZmONE <- list(  meanmONE, Inc,threshold,Thre, covAmONE, covCmONE, covEmONE, covPmONE, covAmsONE, covCmsONE )
parsZoONE <- list( parsZmONE, parsZfONE, meanoONE, Inc,threshold,Thre, signAONE, signCONE, lower, covAosONE, covCosONE, pathRgONE, pathRcONE )
modelMZfONE <- mxModel( parsZfONE, expMeanZfONE, covMZfONE, expCovMZfONE, dataMZfONE, expMZfONE, funML, name="MZfONE" )
modelDZfONE <- mxModel( parsZfONE, expMeanZfONE, covDZfONE, expCovDZfONE, dataDZfONE, expDZfONE, funML, name="DZfONE" )
modelMZmONE <- mxModel( parsZmONE, expMeanZmONE, covMZmONE, expCovMZmONE, dataMZmONE, expMZmONE, funML, name="MZmONE" )
modelDZmONE <- mxModel( parsZmONE, expMeanZmONE, covDZmONE, expCovDZmONE, dataDZmONE, expDZmONE, funML, name="DZmONE" )
modelDZoONE <- mxModel( parsZoONE, expMeanZoONE, covDZoONE, expCovDZoONE, dataDZoONE, expDZoONE, funML, name="DZoONE" )
parsZfTWO <- list(  meanfTWO, Inc,Thre,threshold, covAfTWO, covCfTWO, covEfTWO, covPfTWO )
parsZmTWO <- list(  meanmTWO, Inc,Thre,threshold, covAmTWO, covCmTWO, covEmTWO, covPmTWO, covAmsTWO, covCmsTWO )
parsZoTWO <- list( parsZmTWO, parsZfTWO, meanoTWO, Inc,Thre,threshold, signATWO, signCTWO, lower, covAosTWO, covCosTWO, pathRgTWO, pathRcTWO )
modelMZfTWO <- mxModel( parsZfTWO, expMeanZfTWO, covMZfTWO, expCovMZfTWO, dataMZfTWO, expMZfTWO, funML, name="MZfTWO" )
modelDZfTWO <- mxModel( parsZfTWO, expMeanZfTWO, covDZfTWO, expCovDZfTWO, dataDZfTWO, expDZfTWO, funML, name="DZfTWO" )
modelMZmTWO <- mxModel( parsZmTWO, expMeanZmTWO, covMZmTWO, expCovMZmTWO, dataMZmTWO, expMZmTWO, funML, name="MZmTWO" )
modelDZmTWO <- mxModel( parsZmTWO, expMeanZmTWO, covDZmTWO, expCovDZmTWO, dataDZmTWO, expDZmTWO, funML, name="DZmTWO" )
modelDZoTWO <- mxModel( parsZoTWO, expMeanZoTWO, covDZoTWO, expCovDZoTWO, dataDZoTWO, expDZoTWO, funML, name="DZoTWO" )
multi <- mxFitFunctionMultigroup( c("MZfONE","DZfONE","MZmONE","DZmONE","DZoONE","MZfTWO","DZfTWO","MZmTWO","DZmTWO","DZoTWO") )
# Create Algebra for Variance Components
rowUS <- rep('US',nv)
colUS <- rep(c('VAfONE','VCfONE','VEfONE','SAfONE','SCfONE','SEfONE','VAmONE','VCmONE','VEmONE','SAmONE','SCmONE','SEmONE','rgONE','rcONE',
# Create Confidence Interval Objects
ciACE <- mxCI( c("US") )
# Build Model with Confidence Intervals
modelACEra <- mxModel( "oneACEra5voa", parsZfONE, parsZmONE, parsZoONE, modelMZfONE, modelDZfONE, modelMZmONE, modelDZmONE, modelDZoONE, multi, estUS, ciACE, parsZfTWO, parsZmTWO, parsZoTWO, modelMZfTWO, modelDZfTWO, modelMZmTWO, modelDZmTWO, modelDZoTWO )

Thank you so much in advance