require(OpenMx, polycor) getThresh <- function(df) { tbl <- table(df[,1], df[,2]) p <- (2*tbl[1,1] + tbl[1,2] + tbl[2,1])/(2*sum(tbl)) qnorm(p) } dz <- read.csv('dz.csv') dz[,1] <- mxFactor(dz[,1],c(1,2)) dz[,2] <- mxFactor(dz[,2],c(1,2)) dnames <- colnames(dz)[1:2] dzAlphaStart <- .9*polychor(table(dz[,1:2])) dzBetaStart <- .02 dztStart <- getThresh(dz) # these values should give an error because the # predicted matrix is not positive definite dztStart <- 0.1 dzAlphaStart <- -1.1 dzBetaStart <- 0.6 dzModel <- mxModel("dzModel", mxData(dz, type="raw"), mxMatrix(name="one", type="Iden", 1), mxMatrix(name="means", type="Full", 1, 2, dimnames=list("mean", dnames)), mxMatrix(name="dzt", type="Full", 1, 1, free=TRUE, value=dztStart), mxMatrix(name="alphaDZ", type="Full", 1, 1, free=TRUE, values=dzAlphaStart), mxMatrix(name="betaDZ", type="Full", 1, 1, free=TRUE, values=dzBetaStart), mxAlgebra(name="rPreDZ", alphaDZ + betaDZ*data.SF), mxAlgebra(rbind(cbind(one, rPreDZ), cbind(rPreDZ, one)), name="preDZ", dimnames=list(dnames, dnames)), mxAlgebra(cbind(dzt, dzt), name="dzThreshold", dimnames=list("th", dnames)), mxFIMLObjective(covariance="preDZ", means="means", thresholds="dzThreshold")) dzModel <- mxRun(dzModel) # note that when SF = 1, the predicted covariance # is negative definite preCov <- dzModel$preDZ@result eigen(preCov) det(preCov) # note that pmvnorm chokes on the bivariate normal integral require(mvtnorm) t <- dzModel$dzt@values[1,1] pmvnorm(lower=c(-Inf, -Inf), upper=c(t, t), mean=c(0, 0), sigma=preCov) omxMnor(preCov, c(0, 0), c(-Inf, -Inf), c(t, t)) # this is ok omxMnor(matrix(c(1,0,0,1),2,2), c(0, 0), c(-Inf, -Inf), c(0,0)) # so is this omxMnor(matrix(c(1,0,0,1),2,2), c(0, 0), c(-Inf, -Inf), c(1.282,1.282)) # looks like it should complain about this: omxMnor(matrix(c(1,90,90,1),2,2), c(0, 0), c(-Inf, -Inf), c(1.282,1.282)) # and about this omxMnor(matrix(c(1,90,1,90),2,2), c(0, 0), c(-Inf, -Inf), c(1.282,1.282))