library(OpenMx) # no bounds thisMod <- mxModel(name="testCI") thisMod <- mxModel(thisMod, mxMatrix(name="wt", type="Full", nrow=2, ncol=1, free=FALSE, values=c(147, 147))) thisMod <- mxModel(thisMod, mxMatrix(name="Zobs", type="Full", nrow=2, ncol=1, free=FALSE, values=c(0.5277240, 0.1723423))) thisMod <- mxModel(thisMod, mxMatrix(name="parms", type="Full", nrow=2, ncol=1, free=TRUE, values=c(.4, .05))) thisMod <- mxModel(thisMod, mxMatrix(name="design", type="Full", nrow=2, ncol=2, free=FALSE, values=c(1, .5, 1, 1))) thisMod <- mxModel(thisMod, mxMatrix(name="one", type="Unit", nrow=2, ncol=1)) thisMod <- mxModel(thisMod, mxAlgebra(design %*% parms, name="Rpre")) thisMod <- mxModel(thisMod, mxAlgebra(.5 %x% log((one + Rpre)/(one - Rpre)), name="Zpre")) thisMod <- mxModel(thisMod, mxAlgebra(Zobs - Zpre, name="delta")) thisMod <- mxModel(thisMod, mxAlgebra(sum(wt * delta * delta), name="chisq")) thisMod <- mxModel(thisMod, mxAlgebraObjective("chisq")) thisMod <- mxModel(thisMod, mxCI("parms")) testNoBounds <- mxRun(thisMod, intervals=TRUE) # bounds thatMod <- thisMod thatMod <- mxModel(thatMod, mxMatrix(name="parms", type="Full", nrow=2, ncol=1, free=TRUE, values=c(.4, .05), lbound = c(0, 0))) testBounds <- mxRun(thatMod, intervals=TRUE) # comparison of confidence intervals summary(testNoBounds)$CI summary(testBounds)$CI # examining the condfidence limit of Va over the interval 0:1 Zobs <- thisMod$Zobs@values wt <- thisMod$wt@values fx <- function(x) { Rpre <- c(Va + x, .5*Va + x) if (any(Rpre < -.999) | any(Rpre > .999)) { f <<- 1.1*f return(f) } Zpre <- .5 * log((1 + Rpre)/(1 - Rpre)) delta <- Zobs - Zpre f <- sum(wt * delta * delta) f } un <- matrix(0, 100, 4) cn <- matrix(0, 100, 4) VaSeq <- seq(0, 1, by=.01) for (i in 1:100) { Va <- VaSeq[i] x <- .05 f <- fx(x) + 1 if (Va >= .95) x <- -.1 res <- optim(x, fx, method="BFGS") un[i,] <- c(res$value, Va, res$par, 1 - pchisq(res$value, 1)) if (Va >= .95) x <- .005 f <- fx(x) + 1 res <- optim(x, fx, method="L-BFGS-B", lower=0) cn[i,] <- c(res$value, Va, res$par, 1 - pchisq(res$value, 1)) } # demonstration that the lower CIs are the same for the two solutions cn[1:25,] - un[1:25,]