You are here

Home › Forums › OpenSEM Forums › Behavioral Genetics Models

The OpenMx website will be down for maintenance from 9 AM EDT on Tuesday, September 17th, and is expected to return by the end of the day on Wednesday, September 18th. During this period, the backend will be updated and the website will get a refreshed look.

Sex-Limitation model in two different cohorts

1 post / 0 new
Log in or register to post comments
Wed, 03/30/2022 - 11:21
#1
JuanJMV's picture
JuanJMV
Offline
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 (https://pubmed.ncbi.nlm.nih.gov/15355151/) 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
 
# PREPARE MODEL
# 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), 
                     values=c(0,1),
                     labels=c("t1","inc"), 
                     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',
               'VAfTWO','VCfTWO','VEfTWO','SAfTWO','SCfTWO','SEfTWO','VAmTWO','VCmTWO','VEmTWO','SAmTWO','SCmTWO','SEmTWO','rgTWO','rcTWO'),each=nv)
estUS <- mxAlgebra( expression=cbind(VAfONE,VCfONE,VEfONE,VAfONE/VfONE,VCfONE/VfONE,VEfONE/VfONE,VAmONE+VAmsONE,VCmONE+VCmsONE,VEmONE,(VAmONE+VAmsONE)/VmONE,(VCmONE+VCmsONE)/VmONE,VEmONE/VmONE,rgONE,rcONE,
                                     VAfTWO,VCfTWO,VEfTWO,VAfTWO/VfTWO,VCfTWO/VfTWO,VEfTWO/VfTWO,VAmTWO+VAmsTWO,VCmTWO+VCmsTWO,VEmTWO,(VAmTWO+VAmsTWO)/VmTWO,(VCmTWO+VCmsTWO)/VmTWO,VEmTWO/VmTWO,rgTWO,rcTWO), name="US",
                    dimnames=list(rowUS,colUS))
# 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

Top

Copyright © 2007-2024 The OpenMx Project