Nothing
#
# Copyright 2007-2019 by the individuals mentioned in the source code history
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# -----------------------------------------------------------------------
# Program: UnivariateTwinAnalysis_MatrixRaw.R
# Author: Hermine Maes
# Date: 2009.08.01
#
# ModelType: ACE
# DataType: Twin
# Field: Human Behavior Genetics
#
# Purpose:
# Univariate Twin Analysis model to estimate causes of variation
# Matrix style model input - Raw data input
#
# RevisionHistory:
# Hermine Maes -- 2009.10.08 updated & reformatted
# Ross Gore -- 2011.06.06 added Model, Data & Field metadata
# Hermine Maes -- 2014.11.04 piecewise specification
# -----------------------------------------------------------------------------
require(OpenMx)
# Load Library
# -----------------------------------------------------------------------------
# Load Data
data(twinData)
# Select Variables for Analysis
Vars <- 'bmi'
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") #c('bmi1','bmi2')
# Select Data for Analysis
mzData <- subset(twinData, zyg==1, selVars)
dzData <- subset(twinData, zyg==3, selVars)
# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
# Prepare Data
# -----------------------------------------------------------------------------
require(OpenMx)
# Set Starting Values
svMe <- 20 # start value for means
svPa <- .5 # start value for path coefficients (sqrt(variance/#ofpaths))
# ACE Model
# Matrices declared to store a, d, and e Path Coefficients
pathA <- mxMatrix( type="Full", nrow=nv, ncol=nv,
free=TRUE, values=svPa, label="a11", name="a" )
pathC <- mxMatrix( type="Full", nrow=nv, ncol=nv,
free=TRUE, values=svPa, label="c11", name="c" )
pathE <- mxMatrix( type="Full", nrow=nv, ncol=nv,
free=TRUE, values=svPa, label="e11", name="e" )
# Matrices generated to hold A, C, and E computed Variance Components
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )
# Algebra to compute total variances
covP <- mxAlgebra( expression=A+C+E, name="V" )
# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv,
free=TRUE, values=svMe, label="mean", name="expMean" )
covMZ <- mxAlgebra( expression=rbind( cbind(V, A+C),
cbind(A+C, V)), name="expCovMZ" )
covDZ <- mxAlgebra( expression=rbind( cbind(V, 0.5%x%A+C),
cbind(0.5%x%A+C , V)), name="expCovDZ" )
# Data objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
# Objective objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean",
dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean",
dimnames=selVars )
funML <- mxFitFunctionML()
# Combine Groups
pars <- list( pathA, pathC, pathE, covA, covC, covE, covP )
modelMZ <- mxModel( pars, meanG, covMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, meanG, covDZ, dataDZ, expDZ, funML, name="DZ" )
fitML <- mxFitFunctionMultigroup(c("MZ", "DZ") )
twinACEModel <- mxModel( "ACE", pars, modelMZ, modelDZ, fitML )
# Run Model
twinACEFit <- mxRun(twinACEModel, intervals=T)
twinACESum <- summary(twinACEFit)
twinACESum
# Fit ACE Model with RawData and Matrices Input
# -----------------------------------------------------------------------------
twinACEFit <- mxRun(twinACEModel)
# Run ACE Model
# -----------------------------------------------------------------------------
# Generate ACE Model Output
estMean <- mxEval(mean, twinACEFit) # expected mean
estCovMZ <- mxEval(expCovMZ, twinACEFit$MZ) # expected covariance matrix for MZ's
estCovDZ <- mxEval(expCovDZ, twinACEFit$DZ) # expected covariance matrix for DZ's
estVA <- mxEval(a*a, twinACEFit) # additive genetic variance, a^2
estVC <- mxEval(c*c, twinACEFit) # dominance variance, d^2
estVE <- mxEval(e*e, twinACEFit) # unique environmental variance, e^2
estVP <- (estVA+estVC+estVE) # total variance
estPropVA <- estVA/estVP # standardized additive genetic variance
estPropVC <- estVC/estVP # standardized dominance variance
estPropVE <- estVE/estVP # standardized unique environmental variance
estACE <- rbind(cbind(estVA,estVC,estVE), # table of estimates
cbind(estPropVA,estPropVC,estPropVE))
LL_ACE <- mxEval(objective, twinACEFit) # likelihood of ADE model
# Get Model Output
# -----------------------------------------------------------------------------
Mx.A <- 0.6173023
Mx.C <- 5.595822e-14
Mx.E <- 0.1730462
Mx.M <- 21.39293
Mx.LL_ACE <- 4067.663
# Mx answers hard-coded
# -----------------------------------------------------------------------------
omxCheckCloseEnough(LL_ACE,Mx.LL_ACE,.001)
omxCheckCloseEnough(estVA,Mx.A,.001)
omxCheckCloseEnough(estVC,Mx.C,.001)
omxCheckCloseEnough(estVE,Mx.E,.001)
omxCheckCloseEnough(estMean,Mx.M,.001)
# Compare OpenMx results to Mx results
# -----------------------------------------------------------------------------
# Change Model
twinAEModel <- mxModel( twinACEFit, name="AE" )
twinAEModel <- omxSetParameters( twinAEModel, labels="c11", free=FALSE, values=0 )
twinAEFit <- mxRun(twinAEModel)
# Run AE Model
# -----------------------------------------------------------------------------
# Generate AE Model Output
estVA <- mxEval(a*a, twinAEFit) # additive genetic variance, a^2
estVE <- mxEval(e*e, twinAEFit) # unique environmental variance, e^2
estVP <- (estVA+estVE) # total variance
estPropVA <- estVA/estVP # standardized additive genetic variance
estPropVE <- estVE/estVP # standardized unique environmental variance
estAE <- rbind(cbind(estVA,estVE), # table of estimates
cbind(estPropVA,estPropVE))
LL_AE <- mxEval(objective, twinAEFit) # likelihood of AE model
LRT_ACE_AE <- LL_AE - LL_ACE
# Get Model Output
# -----------------------------------------------------------------------------
Mx.A <- 0.6173023
Mx.C <- 0
Mx.E <- 0.1730462
Mx.M <- 21.39293
Mx.LL_AE <- 4067.663
# Hard-code Mx results
# -----------------------------------------------------------------------------
omxCheckCloseEnough(LL_AE,Mx.LL_AE,.001)
omxCheckCloseEnough(estVA,Mx.A,.001)
omxCheckCloseEnough(estVC,Mx.C,.001)
omxCheckCloseEnough(estVE,Mx.E,.001)
omxCheckCloseEnough(estMean,Mx.M,.001)
# Compare OpenMx results to Mx results:
# -----------------------------------------------------------------------------
estACE
estAE
LRT_ACE_AE
# Print relevant output
# -----------------------------------------------------------------------------
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.