## You are here

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.

# Interpretation of Threshold Matrix, Ordinal Data

1 post / 0 new
Offline
Joined: 04/04/2018 - 14:52
Interpretation of Threshold Matrix, Ordinal Data
AttachmentSize
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 )