Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings, : BLAS/LAPACK routine 'DPOTRF' gave error code -4

Thanks for your quick reply—pardon my delay. And, thank you for your help.

Version: [1] "1.1.1-1785"

Here's the full script that includes where we set AFac, CFac, EFac, and selVar:

twins <- read.csv("twins_nepsy.csv",header=T)

# Telling R that proband variable is not-numeric
twins[,6] <- as.ordered(twins[,6])

## Creating MZ and DZ data sets ##
# Dividing up twins
twinA <- twins[twins$twin=="A",]
twinB <- twins[twins$twin=="B",]
# Remerging by case to create paired data set
newtwins <- merge(twinA, twinB, by=c("case","Zygosity"),all.x=TRUE, all.y=TRUE,suffixes=c("_A","_B"))
# Making data sets of Just MZ & DZ
MZdata <- as.data.frame(subset(newtwins,Zygosity==1))
DZdata <- as.data.frame(subset(newtwins,Zygosity==2))

## Loading Required/Useful Libraries ##

# Load OpenMX
require(OpenMx) #Loads OpenMx
require(psych) #Loads Psych package
# Load additional GenEpi tools directly from VCU
source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R")

# Defining number of variables as 2 for OpenMx
nv <- 2 #Number of phenotypes
ntv <- nv*2

## Defining Variables / Data for use with OpenMX ##

# Making Sparse Data Frame with only phenotype
mzData <- MZdata[,c(5,6,9,10)]
dzData <- DZdata[,c(5,6,9,10)]
# Define Variable to use in OpenMX
selVars <- c(names(mzData))
Vars <- 'var'

# Fit Bivariate ACE Model with Raw Data Input
# -----------------------------------------------------------------------

# Defining start values, may need to be tweaked if you get MX STATUS RED
meanSVnv <- rep(0,nv)
cholASVnv <- rep(.3,nv*(nv+1)/2)
cholCSVnv <- rep(.1,nv*(nv+1)/2)
cholESVnv <- rep(.3,nv*(nv+1)/2)

# Making labels for OpenMx script, not necessary
rows <- c()
cols <- c()
for (i in 1:nv){
row <- c(i:nv)
col <- rep(i,(nv+1-i))
rows <- c(rows,row)
cols <- c(cols,col)
}

AFac <- paste("A",rows,cols,sep="")
CFac <- paste("C",rows,cols,sep="")
EFac <- paste("E",rows,cols,sep="")
mean <- as.character(c("mean1","mean2"))

bivACEModel <- mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholASVnv, labels=AFac, name="a" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholCSVnv, labels=CFac, name="c" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholESVnv, labels=EFac, name="e" ),
# Matrices A, C, and E compute variance components
mxAlgebra( expression=a %*% t(a), name="A" ),
mxAlgebra( expression=c %*% t(c), name="C" ),
mxAlgebra( expression=e %*% t(e), name="E" ),
# Algebra to compute total variances and standard deviations (diagonal only)
mxAlgebra( expression=A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"),
# Confidence intervals portion for covariance matrices
mxAlgebra( expression=A/V,name="stndVCA"),
mxAlgebra( expression=C/V,name="stndVCC"),
mxAlgebra( expression=E/V,name="stndVCE"),
mxCI(c("stndVCA","stndVCC","stndVCE")),
# Confidence intervals for Phenotypic correlation pieces
mxAlgebra((solve(sqrt(I*V)) %*% V %*% solve(sqrt(I*V)))[1,2],name="Vcorr"),
mxAlgebra((solve(sqrt(I*A)) %*% A %*% solve(sqrt(I*A)))[1,2],name="Acorr"),
mxAlgebra((solve(sqrt(I*C)) %*% C %*% solve(sqrt(I*C)))[1,2],name="Ccorr"),
mxAlgebra((solve(sqrt(I*E)) %*% E %*% solve(sqrt(I*E)))[1,2],name="Ecorr"),
mxCI(c("Vcorr","Acorr","Ccorr","Ecorr")),
## Note that the rest of the mxModel statements do not change for bi/multivariate case
# Matrix & Algebra for expected means vector
mxMatrix( type="Full", nrow=1, ncol=nv, free=c(TRUE,FALSE), values=meanSVnv, labels=mean, name="Mean" ),
mxAlgebra( expression= cbind(Mean,Mean), name="expMean"),
## This is the line that allows OpenMx to deal with the ordinal data, it treats it as a normal variable, centered at 0
## and just finds the threshold for a cutoff where all lower is 0 and all higher is 1
mxMatrix(type="Full", nrow=1, ncol=nv, free=c(TRUE), values=0, name="expThre", labels=c("th"),dimnames=list(c('th'),selVars[c(2,4)])),
# Algebra for expected variance/covariance matrix in MZ
mxAlgebra( expression= rbind ( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" ),
# Algebra for expected variance/covariance matrix in DZ
mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ),
mxModel("MZ",
mxData( observed=mzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThre", threshnames=selVars[c(2,4)] )
),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThre", threshnames=selVars[c(2,4)] )
),
mxAlgebra( expression=MZ.objective + DZ.objective, name="neg2sumll" ),
mxAlgebraObjective("neg2sumll")
)
bivACEFit <- mxRun(bivACEModel,intervals=TRUE)
bivACESumm <- summary(bivACEFit)
bivACESumm

Please ignore my previous reply. I've attached the correct script and the simulated data. They reproduce the error on my system.

Thanks again.

~J