Nothing
#------------------------------------------------------------------------------
# Author: Michael D. Hunter
# Date: 2020-03-05 16:53:37
# Filename: RampartDimnames.R
# Purpose: Define a behavior genetics ACE model as a
# Relational SEM (rSEM). This is done for one phenotype (height).
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
library(OpenMx)
library(testthat)
#------------------------------------------------------------------------------
# Prepare Data
# Take wide format data, put it into tall format for multilevel model.
# Also, add coefficient of relatedness "rel" and its inverse "relu"
# Wide to tall
data("twinData", package="OpenMx")
selVars <- c('ht1', 'ht2', 'zyg')
wideData <- subset(twinData, zyg %in% c(1, 3), selVars)
wideData$rel <- c(1, NA, .5)[wideData$zyg]
tallData <- reshape(wideData, varying=list(c('ht1', 'ht2')), v.names=c('ht'), timevar='twin', times=1:2, idvar='fam', direction='long')
# Add a couple new variables, sort, and drop vars that are not needed.
tallData$personID <- 1:nrow(tallData)
tallData$relu <- 1-tallData$rel
tallData <- tallData[order(tallData$fam, tallData$twin), c('fam', 'personID', 'twin', 'rel', 'relu', 'ht')]
head(tallData)
# Split data into "between" and "within" data sets.
wData <- tallData
bData <- tallData[!duplicated(tallData$fam), c('fam', 'rel', 'relu')]
dim(wData)
dim(bData)
dzData <- wideData[wideData$zyg == 3,]
mzData <- wideData[wideData$zyg == 1,]
#------------------------------------------------------------------------------
# Between Model
#Version 1 that writes everything out
bModel1 <- mxModel('between', type="RAM",
mxData(type="raw", observed=bData, primaryKey="fam"),
latentVars = c("C1", "AC1"),
mxPath(c("C1", "AC1"), arrows=2, values=1,
free=FALSE, labels=c(NA, "data.rel")))
bModel1Matrix <- mxModel('mbetween',
mxData(type="raw", observed=bData, primaryKey="fam"),
mxMatrix(name='A', type='Zero', nrow=2, ncol=2),
mxMatrix(name='S', type='Symm', nrow=2, ncol=2, labels=c(NA, NA, 'data.rel'), values=c(1, 0, 1)),
mxMatrix(name='F', type='Zero', nrow=0, ncol=2),
mxExpectationRAM(A='A', S='S', F='F', dimnames=c("C1", "AC1")),
mxFitFunctionML()
)
#------------------------------------------------------------------------------
# Within Model
#Version 1 that writes everything out
wModel1 <- mxModel('within', type="RAM", bModel1,
mxData(type="raw", observed=wData, sort=FALSE),
manifestVars = c('ht'),
latentVars = c("AU1", "E1"),
mxPath(from="one", to=c('ht'), arrows=1, free=TRUE, values=1.5, labels=c("meanHt")),
mxPath(c("AU1", "E1"), arrows=2, values=1,
free=FALSE, labels=c("data.relu", NA)),
mxPath('AU1', c('ht'), values=.8, labels=c('a1Ht')),
mxPath('E1', c('ht'), values=.8, labels=c('e1Ht')),
mxPath('between.C1', c('ht'), values=.8,
labels=c('c1Ht'), joinKey="fam"),
mxPath('between.AC1', c('ht'), values=.8, arrows=1,
labels=c('a1Ht'), joinKey="fam")
)
wModel1Matrix <- mxModel('mwithin', bModel1Matrix,
mxData(type="raw", observed=wData, sort=FALSE),
mxMatrix(name="A", type="Full", nrow=3, ncol=3,
values=c(0, .8, .8,
0, 0, 0,
0, 0, 0),
free=c(FALSE, TRUE, TRUE,
FALSE, FALSE, FALSE,
FALSE, FALSE, FALSE),
labels=c(NA, "a1Ht", "e1Ht", NA, NA, NA, NA, NA, NA),
byrow=TRUE),
mxMatrix(name="S", type="Diag", nrow=3, ncol=3, label=c(NA, "data.relu", NA), values=c(0, 1, 1)),
mxMatrix(name="F", type="Full", nrow=1, ncol=3, values=c(1, 0, 0)),
mxMatrix(name="M", type="Full", nrow=1, ncol=3, values=c(1.5, 0, 0), labels=c("meanHt", NA, NA), free=c(TRUE, FALSE, FALSE)),
mxMatrix(name="cross_model", type="Full", nrow=3, ncol=2,
values=c(.8, .8,
0, 0,
0, 0),
free=c(TRUE, TRUE,
FALSE, FALSE,
FALSE, FALSE),
labels=c("c1Ht", "a1Ht", NA, NA, NA, NA),
byrow=TRUE, dimnames=list(c("ht", "AU1", "E1"), c("C1", "AC1")),
joinKey="fam", joinModel="mbetween"),
mxExpectationRAM(A='A', S='S', F='F', M='M', between="cross_model", dimnames=c("ht", "AU1", "E1")),
mxFitFunctionML()
)
#------------------------------------------------------------------------------
# Run 'em
wRun1 <- mxRun(wModel1)
wRun1Matrix <- mxRun(wModel1Matrix)
coef(wRun1)
coef(wRun1Matrix)
expect_equal(wRun1$output$fit, wRun1Matrix$output$fit, 1e-6)
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.