# Interpretation of Threshold Matrix, Ordinal Data

1 post / 0 new Offline
Joined: 04/04/2018 - 14:52
Interpretation of Threshold Matrix, Ordinal Data
AttachmentSize drawing of distribution that I think the output describes1.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])  0.38561 > pnorm(fitACE$output$algebras$MZ.threG[2,2])
 0.6121303
> pnorm(fitACE$output$algebras$MZ.threG[2,2])-pnorm(fitACE$output$algebras$MZ.threG[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 )