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

5 replies [Last post]
jflournoy's picture
Offline
Joined: 10/07/2011

Hi all, I'm kind of a noob so please forgive me if this is terribly rudimentary.

I get the error message:

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

Traceback yields:

3: .Call("callNPSOL", objective, startVals, constraints, matrices,
parameters, algebras, data, intervalList, communication,
options, state, PACKAGE = "OpenMx")
2: runHelper(model, frontendStart, intervals, silent, suppressWarnings,
unsafe, checkpoint, useSocket, onlyFrontend, useOptimizer)
1: mxRun(bivACEModel, intervals = TRUE)

The code I'm using is:

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")
)

I really appreciate any help!

Cheers,
John

mspiegel's picture
Offline
Joined: 07/31/2009
I found the source of the

I found the source of the bug! I think you have a row of data that contains all missing values for the continuous data, and either missing or non-missing values for the ordinal data (I didn't check). If you have rows that have entirely missing data, you can remove them and the script will work. I'll hand over this bug to Tim, he has more experience with this code. Bug report http://openmx.psyc.virginia.edu/issue/2011/10/handle-row-missing-continu...

tbrick's picture
Offline
Joined: 07/31/2009
Congrats! You found a bug in

Congrats! You found a bug in OpenMx. We'll get right to fixing it, but we'll need a bit more information before we do.

What version of OpenMx are you running? You can just copy and paste the output of the mxVersion() command.

Could you include the lines above what you've listed? Specifically, wherever you set AFac, CFac, EFac, and selVars?

Is this a generated data set, or a real one? If it's generated, could you post it or the script that generates it? If it's real, could you run fakeData() on the dataset you are using and send us the result? The instructions for running fakeData() are here: http://openmx.psyc.virginia.edu/wiki/generating-simulated-data

If we can get those last few details, we'll do our best to fix the problem as quick as we can.

Thanks!

jflournoy's picture
Offline
Joined: 10/07/2011
Files for demonstrating the issue

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

AttachmentSize
fake_mzData.csv 910 bytes
fake_dzData.csv 2.28 KB
bivariate_ACE_w_fakedata.R 6.08 KB
jflournoy's picture
Offline
Joined: 10/07/2011
Thanks for your quick

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

AttachmentSize
fake_dzData.csv 2.28 KB
fake_mzData.csv 906 bytes
mspiegel's picture
Offline
Joined: 07/31/2009
What is "twins_nepsy.csv"?

What is "twins_nepsy.csv"? Also, please include your script as an attachment next time instead of copy/pasting it into the message. This will make it easier to follow the message discussion.