Nothing
# ==================================================================================
# = This is a 1-factor model of 2 continuous variables, and 3 factors =
# = Two definition variables are created (coded with a "c" suffix ?) =
# = In RAM models, you add a definition variable as a fake latent with no variance =
# = with its mean set to "data.defnVar" =
# ==================================================================================
# TODO: Change the file name: "exoPredWLS.R" is not discoverable
library(OpenMx)
library(testthat)
# Make the results repeatable (unclear if this is needed?)
set.seed(1)
data("jointdata", package ="OpenMx", verbose= TRUE)
# ==========================================
# = Here's what our input data looks like: =
# ==========================================
str(jointdata)
# Make z1c = z1 + some noise;
jointdata$z1c <- with(jointdata, z1 * .1 + rnorm(nrow(jointdata)))
# and z2c which is just a column of noise with a mean = the mean factor level of z2 (no clue why as yet...)
jointdata$z2c <- with(jointdata, rnorm(nrow(jointdata), mean=unclass(z2)*.2))
# =============================================
# = Build a function to create WLS RAM models =
# =============================================
# note: This makes the script more complex, but will allow some programatic generation below
buildModel <- function(manifests= paste0('z', 1:5), type = 'WLS', slopes= TRUE) {
jr <- mxModel("JointRAM", type= "RAM",
manifestVars = manifests,
# TODO: is this legit (to use a var name found in the data as a latent name?)
latentVars = c('G','z1c','z2c'),
mxPath('one', c('z1', 'z3'), free= TRUE, labels= c('z1','z3')),
mxPath(paste0('z', c(1,3)), arrows= 2, free= TRUE, values= .5),
# ordinal variables have fixed variance and no mean
# TODO: why var = .5?
mxPath(paste0('z', c(2,4,5)), arrows= 2, free= FALSE, values= .5),
# latent scaled to var of 1 (mean is fixed at zero by default)
mxPath('G', arrows= 2, values= 1, free= FALSE),
mxPath('G', to = manifests, free= TRUE, values= 1, lbound= 0),
mxThreshold('z2', 1, free=TRUE, labels="z2_thresh1"),
mxThreshold('z4', 3, free=TRUE, labels=paste0("z4_thresh",1:3)),
mxThreshold('z5', 2, free=TRUE, labels=paste0("z5_thresh",1:2)),
# Note: Data are raw, despite our using WLS
mxData(jointdata, type = "raw", verbose= 0L),
# Note: this call to mxFitFunctionWLS is all that's
# required to make a model into WLS!
mxFitFunctionWLS(type)
)
# TODO: "slopes" is means?
if(slopes){
jr <- mxModel(jr,
mxPath('one', to = 'z1c', free= FALSE, labels= "data.z1c"),
mxPath('one', to = 'z2c', free= FALSE, labels= "data.z2c"),
mxPath('z1c', to = 'z1', labels= "r1"),
mxPath('z2c', to = 'z2', labels= "r2"))
}
# TODO: Shouldn't we have a call to mxExpectationRAM or LISREL??? here???
# mxExpectationRAM(M = NA, dimnames = NA, thresholds = "T", threshnames = ???)
return(jr)
}
jointRAM1 <- buildModel()
jointRAM1 <- mxRun(jointRAM1)
summary(jointRAM1)
expect_equal(length(mxGetExpected(jointRAM1, 'standvector')),
summary(jointRAM1)$observedStatistics)
# plot(jointRAM1) # (if using umx)
# Where do these come from?
# cat(deparse(round(unname(coef(jointRAM1)), 4)))
param <- c(0.5809, 0.5898, 0.6528, 0.6061, 0.1744, 0.0871, 0.0504, 0.5475,
0.4638, 7.9037, 2.0759, 0.0624, -0.3663, 0.127, 0.7916, -0.6475,
-0.296)
omxCheckCloseEnough(coef(jointRAM1), param, 1e-3)
param.se <- c(0.0613, 0.1056, 0.0684, 0.0941, 0.0665, 0.0541, 0.0648,
0.0559, 0.0645, 0.0736, 0.0593, 0.0745, 0.0776, 0.0726, 0.0921,
0.0655, 0.0585)
omxCheckCloseEnough(c(jointRAM1$output$standardErrors), param.se, 1e-3)
# ===============================================================
# = Switch jointRAM1 from MxFitFunctionWLS to an ML fitfunction =
# = to compare the estimates of these two approaches =
# ===============================================================
jointRAM2 <- mxModel(jointRAM1, mxFitFunctionML())
jointRAM2 <- mxRun(jointRAM2)
summary(jointRAM2)
estCmp <- cbind(coef(jointRAM2), coef(jointRAM1))
omxCheckCloseEnough(cor(estCmp)[2,1], 1, 1e-4)
seCmp <- cbind(jointRAM2$output$standardErrors, jointRAM1$output$standardErrors)
omxCheckCloseEnough(cor(seCmp)[2,1], 1, .04)
# ===============================================================
# = Iterate over model types allowed by the buildModel function =
# ===============================================================
mani = paste0('z', 5:1)
for (slopes in c(TRUE,FALSE)) {
for (type in c('WLS','DWLS','ULS')) {
jm3 <- buildModel(type=type, slopes=slopes)
#jm3$data$verbose <- 1L
jm3 <- mxRun(jm3)
jm4 <- mxModel(buildModel(manifests = mani, type=type, slopes=slopes), jm3$data)
jm4 <- mxRun(jm4)
}
}
# =========================
# = Make a LISREL version =
# =========================
jointLISREL <- mxModel("JointLISREL", type="LISREL",
manifestVars = list(endogenous= paste0('z',1:5)),
latentVars = list(endogenous= c('G','z1c','z2c')),
mxPath('one', paste0('z', c(1,3)), free= TRUE, labels= c('z1','z3')),
mxPath(paste0('z', c(1,3)), arrows= 2, free=TRUE, values= .5),
mxPath(paste0('z', c(2,4,5)), arrows= 2, free=FALSE, values= .5),
mxPath('G', arrows= 2, values= 1, free= FALSE),
mxPath('G', paste0('z', 1:5), free= TRUE, values= 1, lbound= 0, ubound= 4),
mxPath('one', 'z1c', free= FALSE, labels= "data.z1c"),
mxPath('one', 'z2c', free= FALSE, labels= "data.z2c"),
mxPath('z1c', 'z1', labels= "r1"),
mxPath('z2c', 'z2', labels= "r2"),
mxThreshold('z2', 1, free=TRUE, labels="z2_thresh1"),
mxThreshold('z4', 3, free=TRUE, labels=paste0("z4_thresh",1:3)),
mxThreshold('z5', 2, free=TRUE, labels=paste0("z5_thresh",1:2)),
mxData(jointdata, "raw", verbose=0L),
mxFitFunctionWLS()
# Shouldn't we have a call to mxExpectationRAM or LISREL??? here???
# mxExpectationRAM(M = NA, dimnames = NA, thresholds = "T", threshnames = ???)
)
# =================================================================
# = TODO: How would a user ever discover this?: What does it do?? =
# =================================================================
# jointLISREL$expectation$verbose <- 1L
jointLISREL <- mxRun(jointLISREL)
expect_equal(length(mxGetExpected(jointLISREL, 'standvector')),
summary(jointLISREL)$observedStatistics)
# Compare parameter estimates from the RAM and LISREL models
print(max(abs(coef(jointLISREL) - coef(jointRAM1))))
omxCheckCloseEnough(coef(jointLISREL), coef(jointRAM1), 2e-5)
# tmp <- c(jointRAM1$output$standardErrors)
# names(tmp) <- c()
# cat(deparse(round(tmp,4)))
numPeople <- 100
personData <- data.frame(
snp=rnorm(numPeople),
isMale=as.numeric(rbinom(numPeople,1,.5)),
phenotype=rnorm(numPeople),
snpsex=0)
m1 <- mxModel(
"gwsem", type="RAM",
manifestVars = c('phenotype'),
latentVars = c('snp', 'sex', 'snpsex'),
# residual variances
mxPath(c('phenotype'), arrows=2, values=1),
mxPath('one', 'sex', free=FALSE, labels="data.isMale"),
mxPath('one', 'snp', free=FALSE, labels="data.snp"),
mxAlgebra(data.snp * data.isMale, name="snpsexAlg",
dimnames=list(NULL, 'snpsex')),
mxPath('one', 'snpsex', free=FALSE, labels="data.snpsex"),
mxPath('one', 'phenotype'),
mxPath(c('snp','sex', 'snpsex'), 'phenotype'),
mxData(personData, type="raw", algebra='snpsexAlg'),
mxFitFunctionWLS(allContinuousMethod = "marginals"))
m1 <- mxRun(m1)
personData$snpsex <- with(personData, snp * isMale)
c1 <- coef(lm(phenotype ~ isMale + snp + snpsex, personData))
omxCheckCloseEnough(c1[-1], coef(m1)[c('gwsem.A[1,3]', 'gwsem.A[1,2]', 'gwsem.A[1,4]')], 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.