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
Printer-friendly version