#
# 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.
# ===========
# = history =
# ===========
# ~/bin/OpenMx/inst/models/failing/SimpleConfidenceIntervals.R
# 2017-04-14 04:35PM TBATES: update for mxFitFunctionMultigroup(c("MZ", "DZ"))
# 2019-03-24 05:06PM TBATES: Still FAILING at line > 167)
require(OpenMx)
# j repeatedly gets set to a large value at subnp.cpp line 893.
#
# This is the same algorithmic problem with RSOLNP and CSOLNP we talked
# about some time ago. So, CSOLNP needs boundaries for some problems to
# get to the optimum. If you remove the boundary, it will get stuck in a
# loop.
if (mxOption(NULL,"Default optimizer") != 'CSOLNP') stop("SKIP")
# mxOption(NULL, "Default optimizer", "NPSOL")
#Prepare Data
# -----------------------------------------------------------------------
data(twinData)
summary(twinData)
selVars <- c('bmi1','bmi2')
mzfData <- as.matrix(subset(twinData, zyg==1, c(bmi1,bmi2)))
dzfData <- as.matrix(subset(twinData, zyg==3, c(bmi1,bmi2)))
colMeans(mzfData,na.rm=TRUE)
colMeans(dzfData,na.rm=TRUE)
cov(mzfData,use="complete")
cov(dzfData,use="complete")
#Fit ACE Model with RawData and Matrices Input
# -----------------------------------------------------------------------
twinACEModel <- mxModel("twinACE",
mxModel("common",
# Matrices X, Y, and Z to store a, c, and e path coefficients
mxMatrix("Full", nrow=1, ncol=1, free=T, values=.6, lbound=.1, label="a", name="X"),
mxMatrix("Full", nrow=1, ncol=1, free=T, values=sqrt(.6), label="c", name="Y"),
mxMatrix("Full", nrow=1, ncol=1, free=T, values=.6, lbound=.1, label="e", name="Z"),
# Matrices A, C, and E compute variance components
mxAlgebra(expression=X %*% t(X), name="A"),
mxAlgebra(expression=Y %*% t(Y), name="C"),
mxAlgebra(expression=Z %*% t(Z), name="E"),
mxMatrix("Full", nrow=1, ncol=2, free=T, values= 20, lbound=0, ubound=30, label="mean", name="expMean"),
# Algebra for expected variance/covariance matrix in MZ
mxAlgebra(name="expCovMZ",
rbind(cbind(A+C+E, A+C),
cbind(A+C , A+C+E))
),
# Algebra for expected variance/covariance matrix in DZ
# note use of 0.5, converted to 1*1 matrix
mxAlgebra(name="expCovDZ", expression=
rbind(cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E))
)
),
mxModel("MZ",
mxData(observed=mzfData, type="raw"),
mxFitFunctionML(),
mxExpectationNormal("common.expCovMZ", means="common.expMean", dimnames=selVars)
),
mxModel("DZ",
mxData(observed=dzfData, type="raw"),
mxFitFunctionML(),
mxExpectationNormal("common.expCovDZ", means="common.expMean", dimnames=selVars)
),
mxFitFunctionMultigroup(c("MZ", "DZ")),
mxCI(c("common.A", "common.C[,]", "common.E[1,1]"))
)
#Run ACE model
# -----------------------------------------------------------------------
twinACENoIntervals <- mxRun(twinACEModel, suppressWarnings = TRUE)
twinACENoIntervals <- mxRun(twinACEModel)
twinACEFit <- mxRun(twinACEModel, intervals=TRUE)
# summary(twinACEFit)
detail <- twinACEFit$compute$steps[['CI']]$output[['detail']]
print(detail)
omxCheckTrue(is.factor(detail[['side']]))
omxCheckEquals(levels(detail[['side']]), c('upper', 'lower'))
ci <- twinACEFit$output$confidenceIntervals
# cat(deparse(round(ci[,'ubound'],4)))
omxCheckCloseEnough(ci[-2,'lbound'], c(0.556, 0.1537), .01)
omxCheckCloseEnough(ci[,'ubound'], c(0.683, 0.052, 0.1956), .005)
iterateMxRun <- function(model, maxIterations) {
model <- mxOption(model, "Optimality tolerance", 1e-6)
return(iterateMxRunHelper(mxRun(mxModel(model, mxComputeGradientDescent(maxMajorIter=150L))),
maxIterations, 1))
}
iterateMxRunHelper <- function(model, maxIterations, iteration) {
if (length(model$output) > 0 && model$output$status[[1]] == 0) {
return(model)
} else if (iteration < maxIterations) {
return(iterateMxRunHelper(mxRun(mxModel(model, mxComputeGradientDescent(maxMajorIter=150L))),
maxIterations, iteration + 1))
} else {
return(model)
}
}
twinACEIntervals <- twinACEFit
# twinACEIntervals$output <- list()
maxMisfit <- mxEval(objective, twinACEFit)[1,1] + 3.84
CIaupper <- mxModel(twinACEIntervals, name = 'A_CIupper',
mxConstraint(MZ.objective + DZ.objective < maxMisfit),
mxAlgebra(- common.A, name="upperCIa"),
mxFitFunctionAlgebra("upperCIa")
)
CIalower <- mxModel(twinACEIntervals, name = 'A_CIlower',
mxConstraint(MZ.objective + DZ.objective < maxMisfit),
mxAlgebra(common.A, name="lowerCIa"),
mxFitFunctionAlgebra("lowerCIa")
)
runCIalower <- suppressWarnings(iterateMxRun(CIalower, 3))
runCIaupper <- suppressWarnings(iterateMxRun(CIaupper, 3))
CIcupper <- mxModel(twinACEIntervals, name = 'C_CIupper',
mxConstraint(MZ.objective + DZ.objective < maxMisfit),
mxAlgebra(- common.C,name="upperCIc"),
mxFitFunctionAlgebra("upperCIc")
)
runCIcupper <- suppressWarnings(iterateMxRun(CIcupper, 3))
CIeupper <- mxModel(twinACEIntervals, name = 'E_CIupper',
mxConstraint(MZ.objective + DZ.objective < maxMisfit),
mxAlgebra(- common.E,name="upperCIe"),
mxFitFunctionAlgebra("upperCIe")
)
CIelower <- mxModel(twinACEIntervals, name = 'E_CIlower',
mxConstraint(MZ.objective + DZ.objective < maxMisfit),
mxAlgebra(common.E,name="lowerCIe"),
mxFitFunctionAlgebra("lowerCIe")
)
runCIelower <- suppressWarnings(iterateMxRun(CIelower, 3))
runCIeupper <- suppressWarnings(iterateMxRun(CIeupper, 3))
omxCheckCloseEnough(twinACEFit$output$confidenceIntervals[1, 'lbound'], mxEval(common.A, runCIalower), .01)
omxCheckCloseEnough(twinACEFit$output$confidenceIntervals[1, 'ubound'], mxEval(common.A, runCIaupper), .01)
# Can go either way
# omxCheckTrue(is.na(twinACEFit$output$confidenceIntervals[2, 'lbound']))
# Next check still failing under 2.7.9
# Still failing as of OpenMx version: 2.12.2.233 [GIT v2.12.2-233-ga7a310a]
omxCheckCloseEnough(twinACEFit$output$confidenceIntervals[2, 'ubound'], mxEval(common.C, runCIcupper), .001)
omxCheckCloseEnough(twinACEFit$output$confidenceIntervals[3, 'lbound'], mxEval(common.E, runCIelower), .005)
omxCheckCloseEnough(twinACEFit$output$confidenceIntervals[3, 'ubound'], mxEval(common.E, runCIeupper), .005)
twinACEParallel <- omxParallelCI(twinACENoIntervals)
if (0) {
twinACEFit$output$confidenceIntervals - twinACEParallel$output$confidenceIntervals
}
mask <- !is.na(twinACEFit$output$confidenceIntervals)
omxCheckCloseEnough(twinACEFit$output$confidenceIntervals[mask],
twinACEParallel$output$confidenceIntervals[mask], .001)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.