Transcript slides
OpenMx
Frühling Rijsdijk
OpenMx Models
Steven Boker1 Michael Neale2 Hermine Maes2 Paras Mehta3 Michael Wilde4 Timothy Brick1 Jerey Spies1 Michael Spiegel1 Ryne Estabrook1 Sarah Kenny4 John Fox5 Timothy Bates6 • 1University of Virginia; • 2Virginia Commonwealth University; • 3University of Houston; • 4University of Chicago, Argonne National Labs; • 5McMasters University; 6University of Edinburgh
Open Mx is
OpenMx
is free open source software for fitting Structural Equation Models (SEM) to observed data.
OpenMx
offers the features you would expect in an SEM software package, but
OpenMx
works in ways that will make your modeling jobs easier and will allow you to do things that other SEM packages don't.
Integration with R OpenMx
works as an integral part of the R statistical software system. You have available the full power of the R statistical software system for data manipulation, graphics, simulation, and report generation http://openmx.psyc.virginia.edu
Open Mx is
1. A free, full-featured, open source SEM package.
2. Runs on Windows, Mac OS-X, and Linux.
3. Runs inside the R statistical programming environment.
OpenMx features: 1. A new approach to model specification.
2. Allows both path-style and matrix-style scripting.
3. Web-based forums, tutorials, and a wiki.
http://openmx.psyc.virginia.edu
Install OpenMx
What computers run OpenMx?
• You can run OpenMx on computers using Windows XP, Windows Vista, Mac (Intel or PPC) OS-X 10.5 or later, and varieties of Linux (32 and 64 bit).
What do I need to do first?
• In order to install OpenMx, you will need R version 2.9.x or 2.10.1.
How do I install OpenMx?
• • Open up an R session and copy the following line into the R command line: – source('http://openmx.psyc.virginia.edu/getOpenMx.R') execute line • A few lines of R output will scroll by and your OpenMx (and snow and Matrix) library will be installed
Biometrical Genetic Theory Observed Data
model building
Twin Model Summary Statistics System of Linear Equations Covariance Algebra Path Diagrams Predicted Var/Cov of the Model Path Tracing Rules
SEM
Observed Var/Cov of the Data
Path Diagrams for the Classical ACE Twin Model
1
Model for MZ or DZ Pairs Reared Together
E 1 C
e c
1 A
a
1 1 /.5
1 A 1 C
a c
1 E
e
PTwin 1 PTwin 2
Predicted Var-Cov Matrices
Cov MZ
Tw1 Tw2
Tw1
a
2
a
2
c
2
c
2
e
2
Tw2
a
2
a
2
c
2
c
2
e
2
Cov DZ
Tw1 Tw2
a
2
1 2
a
Tw1
2
c
2
c
2
e
2
a
2
1 2
a
2
Tw2
c
2
c
2
e
2
Saturated Twin Model
To get Twin correlations (MZ and DZ) Testing e.g. equality in means and variances across Twin 1 and 2 and MZ and DZ groups
L1
r
Variance / Covariance
L2
S 1 S 2
V1tw1 V1tw2 Mean1 Mean2
Stand 2x2
1
r r
1
Diag 2x2
S
1 0 0
S
2
Variance / Covariance
L1
r
L2
S 1 S 2
V1tw1 V1tw2
S
1 0 0
S
2 * 1
r r
1 *
S
1 0 0
S
2
V1T1 V1T2 Mean1 Mean2 V1T1
S 1 2 S 1 * r * S 2
V1T2
S 1 * r * S 2 S 2 2
Means and Variances Multinormal Probability Density Function: |2πΣ |-n/2 e -.5((xi - μ) Σ-1 (xi - μ)’) make use of all available data get unbiased estimates if missing data are missing at random Use Maximum Likelihood to estimate free Parameters: 2 means, 2 variances, 1 covariance
ACEuniv.R
require(OpenMx) source("GenEpiHelperFunctions.R") # ---------------------------------------------------------------------- # Prepare Data # ---------------------------------------------------------------------- NCdata <- read.table ('N-CortisolNA.csv', header=T, sep=',') names (NCdata) str(NCdata) nv <- 1 ntv <- nv*2 Vars <-('ncomp') selVars <- c('ncomp1','ncomp2') summary(NCdata) mzData <- subset(NCdata, zyg==1, selVars) dzData <- subset(NCdata, zyg==2, selVars) summary(mzData) summary(dzData) colMeans(mzData,na.rm=TRUE) cov(mzData,use="complete") colMeans(dzData,na.rm=TRUE) cov(dzData,use="complete")
Specify and Run Saturated Model with RawData and Matrix-style Input twinSatModel <- mxModel("twinSat", mxModel("MZ", mxMatrix(type="Full", nrow=1, ncol=ntv, free=T, values=c(30,30), name="expMeanMZ"), mxMatrix(type="Diag", nrow=ntv, ncol=ntv, free=TRUE, values=8, lbound=.001, name= "expSDMZ" ), mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=.8, lbound=-.99, ubound=.99, name= "expCorMZ" ), mxAlgebra( expression= expSDMZ %*% expCorMZ %*% expSDMZ, name= "expCovMZ" ), mxData(mzData, type="raw"), mxFIMLObjective("expCovMZ", "expMeanMZ", selVars)), mxModel("DZ", mxMatrix("Full", 1, 2, T, 30, name="expMeanDZ"), mxMatrix(type="Diag", nrow=ntv, ncol=ntv, free=TRUE, values=8, lbound=.001, name="expSDDZ" ), mxMatrix(type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=.6, lbound=-.99, ubound=.99, name="expCorDZ" ), mxAlgebra( expression= expSDDZ %*% expCorDZ %*% expSDDZ, name="expCovDZ" ),
L1
r mxData(dzData, type="raw"), mxFIMLObjective("expCovDZ", "expMeanDZ", selVars)), mxAlgebra(MZ.objective + DZ.objective, name="-2sumLL"), mxAlgebraObjective("-2sumLL")) twinSatFit <- mxRun(twinSatModel) S 1
V1tw1 L2
S 2
V1tw2 Mean1 Mean2
TESTS ASSUMPTION OF EQUAL VARIANCES # ----------------------------------------------------------------------------------------------------------------- # Specify and Run Saturated SubModel 1 equating Variances ACROSS TWINS # --------------------------------------------------------------------------------------------------------------- twinSatModelSub1 <- twinSatModel twinSatModelSub1$MZ$expSDMZ <- mxMatrix("Diag", 2, 2, T, 8, "sdMZ" , name="expSDMZ") twinSatModelSub1$DZ$expSDDZ <- mxMatrix("Diag", 2, 2, T, 8, "sdDZ" , name="expSDDZ") twinSatFitSub1 <- mxRun(twinSatModelSub1) #Check!!
ExpSDMZ <- mxEval(MZ.expSDMZ, twinSatFitSub1) ExpSDMZ ExpSDDZ <- mxEval(DZ.expSDDZ, twinSatFitSub1) ExpSDDZ
TESTS ASSUMPTION OF EQUAL VARIANCES # ----------------------------------------------------------------------------------------------------------------- # Specify and Run Saturated SubModel 1 equating Variances ACROSS TWINS and ZYG # --------------------------------------------------------------------------------------------------------------- twinSatModelSub1 <- twinSatModel twinSatModelSub1$MZ$expSDMZ <- mxMatrix("Diag", 2, 2, T, 8, "sd" , name="expSDMZ") twinSatModelSub1$DZ$expSDDZ <- mxMatrix("Diag", 2, 2, T, 8, "sd" , name="expSDDZ") twinSatFitSub1 <- mxRun(twinSatModelSub1) #Check!!
ExpSDMZ <- mxEval(MZ.expSDMZ, twinSatFitSub1) ExpSDMZ ExpSDDZ <- mxEval(DZ.expSDDZ, twinSatFitSub1) ExpSDDZ
SPECIFY AND RUN ACE MODEL WITH RAWDATA AND MATRIX STYLE INPUT univACEModel <- mxModel ( "univACE", mxModel( "ACE" , # Matrix for expected means vector for MZ and DZ twins mxMatrix("Full", 1, 2, T, 20, "mean", name="expMean"), # Matrices a11, c11, and e11 to store the a, c, and e path coefficients mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=3, label="a11", name="a"), mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=3, label="c11", name="c"), mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=3, label="e11", name="e"), # Matrices A, C, and E to 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 SD
E 1 C
mxAlgebra( expression=A+C+E, name="V" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="sd"), # Algebra to model expected variance/covariance matrix in MZ mxAlgebra(rbind (cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name= "expCovMZ" ), # Algebra to model expected variance/covariance matrix in DZ
e
mxAlgebra(rbind (cbind(A+C+E , .5%x%A+C), cbind(.5%x%A+C , A+C+E)), name= "expCovDZ" ) ),
c PTwin 1 1 a A
mxModel( "MZ" , mxData(mzData, type="raw"), mxFIMLObjective("ACE.expCovMZ", "ACE.expMean",selVars) ), mxModel( "DZ" , mxData(dzData, type="raw"), mxFIMLObjective("ACE.expCovDZ", "ACE.expMean",selVars) ), mxAlgebra(MZ.objective + DZ.objective, name="-2sumLL"), mxAlgebraObjective("-2sumLL") ) # end of model univACE univACEFit <- mxRun(univACEModel) # ---------------------------------------------------------------------- # Generate ACE Model Output # ---------------------------------------------------------------------- LL_ACE <- mxEval(objective, univACEFit) LL_ACE
Extract Information
univTwinSatFit all information in MxModel fitted object univTwinSatFit@algebras univTwinSatFit@algebras$’-2sumll’ univTwinSatFit@submodels univTwinSatFit@submodels$MZ =univTwinSatFit$MZ list of algebras of container model specific algebra list of all items in all child models list of all items in specific child model univTwinSatFit$MZ@matrices univTwinSatFit$MZ@matrices$CholMZ =univTwinSatFit$MZ$CholMZ list of all matrices in specific child model specific matrix in specific child model univTwinSatFit$MZ@algebras univTwinSatFit$MZ@algebras$expCovMZ =univTwinSatFit$MZ$expCovMZ list of all algebras in specific child model specific algebra in specific child model univTwinSatFit$MZ$objective univTwinSatFit$MZ$data objective of specific child model data of specific child model
R Matrix Operators
http://openmx.psyc.virginia.edu/wiki/matrix-operators-and-functions For example: Inverse Transpose Row x col Mx ‘ ~ * Dot Prod .
Kronecker Prod @ Vert adhesion Hor adhesion | _ OpenMx solve() t() %*% * %x% rbind() cbind()
GenEpiHelperFunctions.R
• Function "parameterSpecifations()" prints labels of a MxMatrix with square brackets surrounding free parameters; returns a matrix of strings • Function "expectedMeansCovariances()" prints expected means and expected covariance matrices for all submodels • Function "formatOutputMatrices()" prints matrix with specified labels and number of decimals • Function "formatMatrix()" returns a matrix with specified dimnames and of decimal places • Function "tableFitStatistics()" prints fit statistics with labels for Full Model and list of Nested Models