Login

News

Restricted Model Cholesky; how to drop/fix parameters?

Posted on
No user picture. Lux Whitehaven Joined: 02/02/2015
Forums
OpenMx General Help
Hello everyone,

I'm a graduate student currently learning OpenMX, and I am unfortunately still a novice at OpenMX and genetic modeling in general. I've so far successfully run multivariate Cholesky analyses with OpenMX in R studio, and gotten interpretable data.

However, I am now trying to drop parameters on my current Cholesky model, to get a "best fit" reduced model for my data. I am having trouble understanding how I might do this; I understand it would involve fixing certain parameters at zero, but given my code, I am completely at a loss as to how I might go about this.

If anyone has the time, would they be able to help me understand how? Or, alternatively, point me in the direction of a Cholesky model code that could easily fix parameters? My thanks; my code is shown below.

nvar <- 4 # number of variables
tnvar <-8 # number of variables*max family size
nlower <-tnvar*(tnvar+1)/2 # number of free elements in a lower matrix tnvar*tnvar
Vars <- c("v1","v2","v3","v4")
selVars <- paste(Vars,c(rep(1,nvar),rep(2,nvar)),sep="")
# This yields a list of the 6 variables for analysis
#v11 v21 v31 v41 v12 v22 v32 v42
mzData <- as.matrix(subset(twindata, zygosity==1, selVars))
dzData <- as.matrix(subset(twindata, zygosity==2, selVars))
#Print Descriptive Statistics
#-------------------------------------------------------------------------
summary(mzData)
summary(dzData)
### Fit Multivariate ACE Model with RawData and Matrices Input ###
ACE_Cholesky_Model <- mxModel("ACE_Cholesky",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="a" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="c" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="e" ),
# Matrices A, C, and E compute variance components
mxAlgebra( a %*% t(a), name="A" ),
mxAlgebra( c %*% t(c), name="C" ),
mxAlgebra( e %*% t(e), name="E" ),
# Algebra to compute total variances and standard deviations (diagonal only)
mxAlgebra( A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"),
mxAlgebra( solve(sqrt(I*V)), name="iSD"),
# Matrix & Algebra for expected means vector
mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 0, name="Mean" ),
mxAlgebra( cbind(Mean,Mean), name="expMean"),
# Algebra for expected variance/covariance matrix in MZ
mxAlgebra( rbind ( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" ),
# Algebra for expected variance/covariance matrix in DZ
mxAlgebra( 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" ),
mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars )
),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars )
),
mxAlgebra( MZ.objective + DZ.objective, name="modelfit" ),
mxAlgebraObjective("modelfit")
)
ACE_Cholesky_Fit <- mxRun(ACE_Cholesky_Model)
ACE_Cholesky_Summ <- summary(ACE_Cholesky_Fit)
ACE_Cholesky_Summ
### Generate Multivariate Cholesky ACE Output ###
parameterSpecifications(ACE_Cholesky_Fit)
expectedMeansCovariances(ACE_Cholesky_Fit)
tableFitStatistics(ACE_Cholesky_Fit)
ACEpathMatrices <- c("ACE.a","ACE.c","ACE.e","ACE.iSD","ACE.iSD %*% ACE.a","ACE.iSD %*% ACE.c","ACE.iSD %*% ACE.e")
ACEpathLabels <- c("pathEst_a","pathEst_c","pathEst_e","sd","stPathEst_a","stPathEst_c","stPathEst_e")
formatOutputMatrices(ACE_Cholesky_Fit,ACEpathMatrices,ACEpathLabels,Vars,3)
ACEcovMatrices <- c("ACE.A","ACE.C","ACE.E","ACE.V","ACE.A/ACE.V","ACE.C/ACE.V","ACE.E/ACE.V")
ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E")
formatOutputMatrices(ACE_Cholesky_Fit,ACEcovMatrices,ACEcovLabels,Vars,3)
# Genetic and Environmental Correlations
ACEcorMatrices <- c("solve(sqrt(ACE.I*ACE.A)) %*% ACE.A %*% solve(sqrt(ACE.I*ACE.A))",
"solve(sqrt(ACE.I*ACE.C)) %*% ACE.C %*% solve(sqrt(ACE.I*ACE.C))",
"solve(sqrt(ACE.I*ACE.E)) %*% ACE.E %*% solve(sqrt(ACE.I*ACE.E))",
"solve(sqrt(ACE.I*ACE.V)) %*% ACE.V %*% solve(sqrt(ACE.I*ACE.V))")
ACEcorLabels <- c("corA","corC","corE","cor")
formatOutputMatrices(ACE_Cholesky_Fit,ACEcorMatrices,ACEcorLabels, Vars, 4)

Replied on Thu, 03/19/2015 - 12:02
Picture of user. neale Joined: 07/31/2009

If possible, try not to cut down a Cholesky decomposition to find the best fit. The Cholesky is typically not a model per se but just an algorithmic trick to keep the predicted covariance matrix (and indeed the genetic/environmental ones) positive definite. It would be much better to start from a common pathway model with a given number of factors, and reduce parameters from there. Better still would be to have a hypothesized factor structure, and see if that fits as parsimoniously (say by AIC) as the Cholesky.
Log in or register to post comments
Replied on Thu, 03/19/2015 - 15:52
Picture of user. RobK Joined: 04/19/2011

Note that the model in his script fits a tetraphenotype ACE twin model. The only "Cholesky" part is that each 4x4 variance component (A, C, & E) is parametrized in terms of its Cholesky factorization.
Log in or register to post comments
Replied on Thu, 03/19/2015 - 15:43
Picture of user. RobK Joined: 04/19/2011

It sounds like you want to use omxSetParameters(). If you have the OpenMx package loaded into R, you can search for help on this function with

?omxSetParameters

For instance, suppose you want to drop the shared-environmental "C" component. To use omxSetParameters(), you would need to know the labels of the parameters you want to drop. You don't assign any labels to your lower-triangular "c" matrix, but some labels should automatically have been assigned to its free parameters at runtime, so you should see something if you enter ACE_Cholesky_Fit$c@labels at the R prompt. The call to omxSetParameters() would look something like

AE_Cholesky_Model <- omxSetParameters(ACE_Cholesky_Fit, labels= #insert labels here
, free=FALSE, values=0, name="AE_Cholesky")

Here's an example for a similar case, except involving only one phenotype.

BTW, lengthy R scripts posted to the forums are much easier to read if you attach them to your post as a text file, or if you enclose the syntax in the appropriate HTML tags.

Log in or register to post comments
Replied on Mon, 03/23/2015 - 15:55
No user picture. Lux Whitehaven Joined: 02/02/2015

My thanks to the both of you for your insight. The model I'm running is actually a three factor model; I just adapted the Cholesky code for a four factor model that I found on the forums, changing the number of Vars to 3, and total Vars to 6.

The reason I'm using a Cholesky and trying to cut it down to a "best fit" is that I'm trying to replicate the methods of a paper which had previously studied the same thing. This paper is also the only paper at the moment which studies the same construct with genetic analyses.

I will certainly try the ?omxSetParameters. With three factors, however, do you suppose a Common Pathways model would be significantly more fruitful, or indeed, easier, in terms of getting a best fit model?

Thanks so much again.

Log in or register to post comments
Replied on Fri, 03/27/2015 - 21:29
No user picture. Lux Whitehaven Joined: 02/02/2015

I've tried using the omxSetParameters function, but there doesn't appear to be any labels. Closer inspection shows that the "labels" within the ACE_Cholesky_Fit are all "NA", which stops me from setting any parameters in the first place.

>ACE_Cholesky_Fit$c@labels
Error: trying to get slot "labels" from an object of a basic class ("NULL") with no slots

Is there a way to go about labelling each pathway so I can fix an individual parameter? As I see it, for example I'm trying to fix only the "stPathEst_c2" of the Cholesky output, and then run it again in that state. Here's an example of my output:

[1] "Matrix ACE.iSD %*% ACE.c"
stPathEst_c1 stPathEst_c2
ffswanhim 0.418 0.000
ffswanin 0.270 0.034
LogSCT 0.117 0.091

stPathEst_c3
ffswanhim 0.000
ffswanin 0.000
LogSCT 0.000

Is "stPathEst_c2" a label I could just go ahead and fix with a code? Like
ACE_Cholesky_Model <- omxSetParameters(ACE_Cholesky_Fit, labels = stpathEst_c2, free=FALSE, values=0, name"ACE_Cholesky") ?

But then, how would I even fix more specific pathways, like just c2 for "LogSCT"?

Log in or register to post comments
Replied on Mon, 03/30/2015 - 10:44
Picture of user. RobK Joined: 04/19/2011

Sorry, I was mistaken in thinking that OpenMx would assign labels to unlabeled free parameters at runtime. If the umx library has a convenient function to add labels to your matrices, then I second Tim Bates' suggestion above to use umx.

Or, you could do it manually, e.g. define your "c" matrix as

mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.6, name="c",
labels=c("c11","c12","c13","c14","c22","c23","c24","c33","c34","c44"))

I missed the fact that your MxModel "ACE_Cholesky" contains a submodel "ACE." To access the labels of your "c" matrix, I think you would (assuming you're using OpenMx version 1.x) do ACE_Cholesky_Fit$ACE$c@labels.

Everything in your script after ACE_Cholesky_Summ <- summary(ACE_Cholesky_Fit) appears to be output formatting, and is extraneous to OpenMx's actual model-fitting. So, no, "stPathEst_c2" is not a label that OpenMx would understand.

Is "LogSCT" one of your phenotypes?

Log in or register to post comments
Replied on Sat, 03/28/2015 - 04:26
Picture of user. tbates Joined: 07/31/2009

One way to label the cells in your model is to use the umx library of helper functions.

It contains umxLabel() wbich, given a matrix model as input labels each cell with the matrix name, row and column number like "a_r1_c2"

Play with it and its examples a bit to see how it works.

http://tbates.github.io

Log in or register to post comments
  Login

News