Attachment | Size |
---|---|
Postscript.txt [6] | 4.57 KB |
I still consider myself a BG novice, and just when I think I've gone as deep as I'll ever go, something else comes along! Awhile ago, the folks on this board graciously helped me do a bivariate Cholesky decomposition of two variables (caster and ogo). The two variables were assessed at different time points, so one variable needed to be adjusted for age1 and sex, the other variable for age2 and sex. We came up with a successful script (attached).
I am now interested in adapting this work so that the "ogo" variable is an ordinal variable, while caster is still continuous. I've read some about thresholds and mxFactor and such, and I've found a number of example scripts that seem to be doing what I'm trying to do (minus the covariates)....still, I know that something is clearly wrong in what I've put together, as I keep getting the error message "Unknown thresholds name 'expThre' detected in the expectation function of model 'MZ'."
mzData <- subset(mydata, ZYGOSITY=="0") dzData <- subset(mydata, ZYGOSITY=="1") # Select Variables for Analysis Vars <- c('caster','ogo') nv <- 2 # number of variables ntv <- nv*2 # number of total variables selVars <- c("caster1", "ogo1", "caster2", "ogo2") bivACEModel <- mxModel("bivACE", mxModel("ACE", # Matrices a, c, and e to store a, c, and e path coefficients mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=c("a11","a21","a22"), name="a" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=c("c11","c21","c22"),name="c" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=c("e11","e21","e22"),name="e" ), # Matrices a, c, and e to store a, c, and e path coefficients mxAlgebra( expression=a %*% t(a), name="A" ), mxAlgebra( expression=c %*% t(c), name="C" ), mxAlgebra( expression=e %*% t(e), name="E" ), # Declare a matrix for the definition variable regression parameters, called beta mxMatrix( type="Full", nrow=1, ncol=2, free=c(TRUE,FALSE), values= c(.4,0), label=c("bage_caster", "bage_ogo"), name="betaage"), mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .1, label=c("bsex_caster", "bsex_ogo"), name="betasex"), mxMatrix( type="Full", nrow=1, ncol=2, free=c(FALSE,TRUE), values= c(0,.1), label=c("bage2_caster", "bage2_ogo"), name="betaage2"), # Matrices A, C, and E compute variance components mxAlgebra( expression=A+C+E, name="V" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"), mxConstraint(V[1,2]==1, name="Vconst"), mxAlgebra(iSD %*% a, name="sta"), mxAlgebra(iSD %*% c, name="stc"), mxAlgebra(iSD %*% e, name="ste"), mxAlgebra(A/V, name="StandardizedA"), mxAlgebra(C/V, name="StandardizedC"), mxAlgebra(E/V, name="StandardizedE"), mxAlgebra(solve(sqrt(I*A)) %&% A, name="CorA"), mxAlgebra(solve(sqrt(I*C)) %&% C, name="CorC"), mxAlgebra(solve(sqrt(I*E)) %&% E, name="CorE"), mxAlgebra(solve(sqrt(I*V)) %&% V, name="CorP"), # Matrix & Algebra for expected means vector mxMatrix( type="Full", nrow=1, ncol=ntv, free=c(T,F,T,F), values=c(.6,0,.6,0), labels=c("caster1","ogo1","caster2","ogo2"), name="expMean"), mxMatrix(type="Full", nrow=1, ncol=2, free=T, values=0.5, label="thre", name="expThre"), # Algebra for expected variance/covariance matrix in MZ mxAlgebra( expression= rbind ( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ), # Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) ), mxModel("MZ", mxData( observed=mzData, type="raw" ), # Algebra for making the means a function of the definition variables age and sex mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.age2T1", "data.age2T2"), name="obs_age2" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ), mxAlgebra( expression= ACE.expMean + obs_sex%x%ACE.betasex + obs_age%x%ACE.betaage + obs_age2%x%ACE.betaage2 , name="expMean"), mxFIMLObjective( covariance="ACE.expCovMZ", means="expMean", dimnames=selVars, thresholds="expThre", threshnames=c('ogo1','ogo2') ) ), mxModel("DZ", mxData( observed=dzData, type="raw" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.age2T1", "data.age2T2"), name="obs_age2" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ), mxAlgebra( expression=ACE.expMean + obs_sex%x%ACE.betasex + obs_age%x%ACE.betaage + obs_age2%x%ACE.betaage2 , name="expMean"), mxFIMLObjective( covariance="ACE.expCovDZ", means="expMean", dimnames=selVars, thresholds="expThre", threshnames=c('ogo1','ogo2') ) ), mxAlgebra( expression=MZ.objective + DZ.objective, name="minus2sumloglikelihood" ), mxAlgebraObjective("minus2sumloglikelihood"), mxCI(c('ACE.sta', 'ACE.stc', 'ACE.ste')), mxCI(c('ACE.StandardizedA', 'ACE.StandardizedC', 'ACE.StandardizedE')), mxCI(c('ACE.CorA', 'ACE.CorC', 'ACE.CorE', 'ACE.CorP')) )