| Attachment | Size |
|---|---|
 drawing of distribution that I think the output describes | 1.3 MB |
Hello,
I'm working with binary/ordinal models after many years of only working with continuous data in OpenMx. I have working scripts for liability threshold ACE models for both a binary outcome, and an ordinal outcome with three levels. I've referenced some of the lectures from prior IBG workshops but want to make sure I am interpreting my threshold output properly.
So a threshold can be interpreted as a Z-score, so if my threshold in a binary model is -1.76, that would correspond to a percentile of 3.9, meaning that 3.9% of my sample falls under the threshold and 96.1% are above the threshold on the underlying liability dimension?
And on the other hand, for an ordinal variable, if I have 3 categories, I am estimating two parameters the first threshold, and then the distance between threshold 1 and 2 (assuming I've gone the route of fixing the liability distribution mean/variance to identify the model). Code for my ordinal script is below.
Is it correct that thinG is freely estimating the first threshold and the increment between threshold 1 and 2, and then threG is the an algebra that constrains the thresholds to be increasing (inc %*% thinG) , and then the resulting threG is the first and second threshold, both of which I can interpret as Z scores?
Here's the ordinal output for thinG and threG
> fitACE$output$matrices$MZ.thinG [,1] [,2] [1,] -0.2907795 -0.2907795 [2,] 0.5756552 0.5756552 > fitACE$output$algebras$MZ.threG [,1] [,2] [1,] -0.2907795 -0.2907795 [2,] 0.2848757 0.2848757 > pnorm(fitACE$output$algebras$MZ.threG[1,1]) [1] 0.38561 > pnorm(fitACE$output$algebras$MZ.threG[2,2]) [1] 0.6121303 > pnorm(fitACE$output$algebras$MZ.threG[2,2])-pnorm(fitACE$output$algebras$MZ.threG[1,1]) [1] 0.2265203
So for this output, threshold 1 is -0.29 and threshold 2 is 0.28, that would mean ~39% of the distribution falls before threshold 1 and ~61% are before threshold 2? And the increment between the z scores for threshold 1 and 2 is 0.57, corresponding to 22% falling between thresholds 1 and 2? I've also attached a drawing of the distribution with the output mapped on, mainly for my own benefit, because this is throwing it back to stats 101 and its been a while.
# Select variables #note phenotype is an ordered factor with 3 levels Vars <- c("Phenotype") nv <- length(Vars) ntv <- nv*2 nth <- 2 selVars <- paste(Vars, c(rep(1,nv), rep(2,nv)), sep="") DefVars <- c("sex", "agetwinpair") # Select Data for Analysis MzData <- subset(TUKUni_widedf, zyg == 1, select = c(selVars, DefVars)) DzData <- subset(TUKUni_widedf, zyg == 2, select = c(selVars, DefVars)) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, values=0, labels="mean", name="meanG" ) regCoefSex <- mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.2, labels="Sexbeta", name="regCoefSex") regCoefAge <- mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.2, labels="Agebeta", name="regCoefAge") SexCov <- mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, labels="data.sex", name="SexCov") AgeCov <- mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, labels="data.agetwinpair", name="AgeCov") means2 <- mxAlgebra(expression=(meanG+(regCoefSex%*%SexCov)+(regCoefAge%*%AgeCov)), name="means2") expMean <- mxAlgebra(t(rbind(means2, means2)), name="expMean") thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntv, free=TRUE, values=.4, labels=c("threshold", "increment", "threshold", "increment"), name="thinG", byrow = F ) inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ) threG <- mxAlgebra( expression= inc %*% thinG, name="threG" ) # Create Matrices for Path Coefficients A <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=.5, labels="a", name="A" ) C <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=.5, labels="c", name="C" ) E <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=.7, labels="e", name="E" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression= A+C+E, name="V" ) covMZ <- mxAlgebra( expression= A+C, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" ) #constrain variance of binary variables var1 <- mxConstraint(expression=diag2vec(V)==1, name="var1") # Putting MZ/DZ data into OpenMx objects DataMZ <- mxData( observed=MzData, type="raw" ) DataDZ <- mxData( observed=DzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="threG" ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG" ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list( meanG, means2, expMean, threG, A, C, E, covP, var1, thinG, inc, regCoefSex, regCoefAge, SexCov, AgeCov ) modelMZ <- mxModel( pars, covMZ, expCovMZ, DataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, covDZ, expCovDZ, DataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Build Model with Confidence Intervals modelACE <- mxModel( "oneACEc", modelMZ, modelDZ, multi ) # Run ACE Model fitACE <- mxRun( modelACE ) sumACE <- summary( fitACE )
Printer-friendly version