1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
|
# Test result of algebras that depend on definition variables
# after mxRun. Contributed by Gregory Carey.
# full text for "Sat.fitFuntion not finite" is
# Error: The job for model 'Sat' exited abnormally with the error message: MxComputeGradientDescent: fitfunction Sat.fitfunction is not finite ()
# In addition: Warning message:
# In model 'Sat' Optimizer returned a non-zero status code 10. Starting values are not feasible. Consider mxTryHard()
#set.seed(101276) # gives "Sat.fitFunction not finite" works w mod 1
#set.seed(50074) # gives "Sat.fitFunction not finite" does not work w mod1
set.seed(120450) # gives "Sat.fitFunction not finite" This works w mod1
#set.seed(270156) # gives "Sat.fitFunction not finite"
h1sq <- .62
h2sq <- .47
rg <- .55
re <- .34
covG <- sqrt(h1sq)*sqrt(h2sq)*rg
rp <- covG + sqrt((1 - h1sq)*(1 - h2sq))*re
pmat <- matrix(c(1, rp, rp, 1), 2, 2)
gmat <- matrix(c(h1sq, covG, covG, h2sq), 2, 2)
SigmaMZ <- rbind(cbind(pmat, gmat),
cbind(gmat, pmat))
SigmaDZ <- rbind(cbind(pmat, .5*gmat),
cbind(.5*gmat, pmat))
thresh <- c(.3, .5, .3, .5)
N <- 250
fcn <- function(N, Sigma, thresh) {
age <- round(25 + 5*rnorm(N))
X <- mvtnorm::rmvnorm(N, rep(0, ncol(Sigma)), Sigma)
for (i in 1:4) {
one <- which(X[,i] > thresh[i])
zero <- which(X[,i] <= thresh[i])
X[one,i] <- 1
X[zero,i] <- 0
}
return(cbind(X, age, age))
}
temp <- fcn(N, SigmaMZ, thresh)
df <- cbind(rep(1, N), temp)
temp <- fcn(N, SigmaMZ, thresh)
temp <- cbind(rep(2, N), temp)
df <- rbind(df, temp)
temp <- fcn(N, SigmaDZ, thresh)
temp <- cbind(rep(3, N), temp)
df <- rbind(df, temp)
temp <- fcn(N, SigmaDZ, thresh)
temp <- cbind(rep(4, N), temp)
df <- rbind(df, temp)
colnames(df) <- c("zyg", "mj1", "nssi1", "mj2",
"nssi2", "age1", "age2")
mj_nssi <- data.frame(df, stringsAsFactors=F)
require(OpenMx)
# --- this from Brooke's post
# http://openmx.ssri.psu.edu/thread/3996
# note that start values have been changed to
# reflect the parameters for the simulated data
#Make Factors
mj_nssi$mj1 <-mxFactor(mj_nssi$mj1, levels=c(0:1), exclude = NA)
mj_nssi$mj2 <-mxFactor(mj_nssi$mj2, levels=c(0:1), exclude = NA)
mj_nssi$nssi1 <-mxFactor(mj_nssi$nssi1, levels=c(0:1), exclude = NA)
mj_nssi$nssi2 <-mxFactor(mj_nssi$nssi2, levels=c(0:1), exclude = NA)
#Selects Vars for Analysis
vars <-c('nssi','mj')
selVars <-c('nssi1', 'mj1', 'nssi2', 'mj2' )
useVars <-c('nssi1', 'mj1', 'nssi2', 'mj2', 'age1', 'age2')
# Select Data for Analysis
mzfData <- subset(mj_nssi, zyg==1, useVars)
mzmData <- subset(mj_nssi, zyg==2, useVars)
dzfData <- subset(mj_nssi, zyg==3, useVars)
dzmData <- subset(mj_nssi, zyg==4, useVars)
nv <- 2 # number of variables per twin
ntv <- nv*2 # number of variables per pair
nth <- 1 # number of max thresholds
# 1) Fits a constrained Polychoric correlation model
# TH same across twins but different across zyg groups
# Age effect is different across variables, but same across thresholds within variables (if c>2) ###?
# There is one overall rPH between var1-2 and the x-trait x-twin correlations are symmetric
# ------------------------------------------------------------------------------------------------------------------------------
# CREATE LABELS & START VALUES as objects(to ease specification)
LabThMZf <-c('Tmz_11f','imz_11f') # THs for var 1 and 2 for a twin individual (mz)
LabThMZm <-c('Tmz_11m','imz_11m')
LabCorMZf <-c('r21f','rMZ1f','MZxtxtf','MZxtxtf','rMZ2f','r21f') #labels for correlation matrix
LabCorMZm <-c('r21m','rMZ1m','MZxtxtm','MZxtxtm','rMZ2m','r21m')
LabThDZf <-c('Tdz_11f','idz_11f')
LabThDZm <-c('Tdz_11m','idz_11m')
LabCorDZf <-c('r21f','rDZ1f','DZxtxtf','DZxtxtf','rDZ2f','r21f')
LabCorDZm <-c('r21m','rDZ1m','DZxtxtm','DZxtxtm','rDZ2m','r21m')
LabCovf <-c('Thnnsif', 'BageThnnsif', 'BageThmjf', 'BageThmjf')
LabCovm <-c('BageThnnsim', 'BageThnnsim', 'BageTmjm', 'BageTmjm')
ThPat <-c(T,T)
# StCorMZf <-c(.1, .6, .2, .2, .6, .1)
# StCorMZm <-c(.1, .6, .2, .2, .6, .1)
# StCorDZf <-c(.1, .2, .1, .1, .2, .1)
# StCorDZm <-c(.1, .2, .1, .1, .2, .1)
StCorMZf <-c(rp, h1sq, covG, covG, h2sq, rp)
StCorMZm <-c(rp, h1sq, covG, covG, h2sq, rp)
StCorDZf <-c(rp, .5*h1sq, .5*covG, .5*covG, .5*h2sq, rp)
StCorDZm <-c(rp, .5*h1sq, .5*covG, .5*covG, .5*h2sq, rp)
#StTH <-c(1.2816, 0.9741)
StTH <- thresh[1:2]
# Define definition variables to hold the Covariates
obsAge1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age1"), name="Age1")
obsAge2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age2"), name="Age2")
# Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations
Mean <-mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" )
betaAf <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCovf, name="BageTHf" )
#betaAm <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=LabCovm, name="BageTHm" )
betaAm <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=T, values=.2, labels=c("Thnnsim", "BageThnnsim"), name="BageTHm" )
inc <-mxMatrix( type="Lower",nrow=nth, ncol=nth, free=F, values=1, name="Low")
Tmzf <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound= -3, ubound=3, labels=LabThMZf, name="ThMZf")
Tmzm <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound= -3, ubound=3, labels=LabThMZm, name="ThMZm") #BMH now give only one lower limit
ThresMZf <-mxAlgebra( expression= cbind(Low%*%ThMZf + BageTHf%x%Age1, Low%*%ThMZf + BageTHf%x%Age2), name="expThresMZf")
ThresMZm <-mxAlgebra( expression= cbind(Low%*%ThMZm + BageTHm%x%Age1, Low%*%ThMZm + BageTHm%x%Age2), name="expThresMZm")
CorMZf <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZf, labels=LabCorMZf, lbound=-.99, ubound=.99, name="expCorMZf")
CorMZm <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorMZm, labels=LabCorMZm, lbound=-.99, ubound=.99, name="expCorMZm")
Tdzf <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound= -3, ubound=3, labels=LabThDZf, name="ThDZf") #BMH now give only one lower limit
Tdzm <-mxMatrix( type="Full", nrow=nth, ncol=nv, free=ThPat, values=StTH, lbound= -3, ubound=3, labels=LabThDZm, name="ThDZm")
ThresDZf <-mxAlgebra( expression= cbind(Low%*%ThDZf + BageTHf%x%Age1, Low%*%ThDZf + BageTHf%x%Age2), name="expThresDZf")
ThresDZm <-mxAlgebra( expression= cbind(Low%*%ThDZm + BageTHm%x%Age1, Low%*%ThDZm + BageTHm%x%Age2), name="expThresDZm")
CorDZf <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZf, labels=LabCorDZf, lbound=-.99, ubound=.99, name="expCorDZf")
CorDZm <-mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=T, values=StCorDZm, labels=LabCorDZm, lbound=-.99, ubound=.99, name="expCorDZm")
# Data objects for Multiple Groups
dataMZf <- mxData( observed=mzfData, type="raw" )
dataMZm <- mxData( observed=mzmData, type="raw" )
dataDZf <- mxData( observed=dzfData, type="raw" )
dataDZm <- mxData( observed=dzmData, type="raw" )
expMZf <- mxExpectationNormal( covariance="expCorMZf", means="M", dimnames=selVars, thresholds="expThresMZf" )
expMZm <- mxExpectationNormal( covariance="expCorMZm", means="M", dimnames=selVars, thresholds="expThresMZm" )
expDZf <- mxExpectationNormal( covariance="expCorDZf", means="M", dimnames=selVars, thresholds="expThresDZf" )
expDZm <- mxExpectationNormal( covariance="expCorDZm", means="M", dimnames=selVars, thresholds="expThresDZm" )
#FitFunctions
fun <- mxFitFunctionML()
fitFunction <- mxFitFunctionMultigroup(c("MZf", "MZm", "DZf", "DZm"))
# Combine Groups
modelMZf <- mxModel( obsAge1, obsAge2, Mean, betaAf, Tmzf, inc, ThresMZf, CorMZf, dataMZf, expMZf, fun, name="MZf" )
modelMZm <- mxModel( obsAge1, obsAge2, Mean, betaAm, Tmzm, inc, ThresMZm, CorMZm, dataMZm, expMZm, fun, name="MZm" )
modelDZf <- mxModel( obsAge1, obsAge2, Mean, betaAf, Tdzf, inc, ThresDZf, CorDZf, dataDZf, expDZf, fun, name="DZf" )
modelDZm <- mxModel( obsAge1, obsAge2, Mean, betaAm, Tdzm, inc, ThresDZm, CorDZm, dataDZm, expDZm, fun, name="DZm" )
SatModel <- mxModel( "Sat", modelMZf, modelMZm, modelDZf, modelDZm, fitFunction)
# my code: where the madness starts
SatModel$MZf$Age1
SatModel$MZf$BageTHf
mxEval(MZf.BageTHf %x% MZf.Age1, SatModel, compute=T)
mxEval(matrix(1, 1, 1) %*% MZf.Age1, SatModel, compute=T)
# -------------------------------------------------------------------------------------------------------------------------------
# 1) RUN Saturated Model
SatFit <- mxRun(SatModel, intervals=F)
summary(SatFit)
# more of my code
mod1 <- mxModel(obsAge1, obsAge2, Mean, betaAf, Tmzf, inc,
ThresMZf, CorMZf, dataMZf, expMZf, fun,
name="MZf", mxFitFunctionML())
mod1Res <- mxRun(mod1)
if (mxOption(NULL, "Default optimizer") == 'NPSOL' &&
max(abs(c(mod1Res$ThMZf$values) - c(0.779, 0.341))) < .01) {
# NPSOL can get stuck in a local minimum
mod1Res <- mxTryHard(mod1Res)
}
# definition variables and dependencies should be set to NA upon return
omxCheckEquals(mod1Res$Age1$values, dataMZf$observed[1,'age1'])
omxCheckEquals(mod1Res$Age2$values, dataMZf$observed[1,'age2'])
omxCheckCloseEnough(c(mod1Res$ThMZf$values), c(0.732, 0.305), .01)
omxCheckCloseEnough(c(mod1Res$BageTHf$values), c(-0.008, 0), .01)
omxCheckEquals(nrow(mod1Res$data$observed), 250)
omxCheckEquals(mod1Res$expThresMZf$result,
mxEval(expThresMZf, mod1Res, compute=T, defvar.row = 1))
|