Fri, 10/07/2011 - 18:27

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

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...

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!

still seems to fail under OpenMx 2.0b

This script reads the data from URL, uses OpenMx 2 syntax, has what might be slightly better starts, and is reduced to highlight the problem

It fails because the means and variances are way off. Like 0 instead of 100. This is really pushing likelihoods into effectively zero range, so the error:

occurs. We can tell that this happened from the start because the number of evaluations is only 2:

Hint to my fellow developers: It would be polite of OpenMx to notify the user that bad starting values appear to be the cause of the problem. This could be done whenever the following conditions hold: fit function is not definite, 2 evaluations, number of free parameters > 0.

Miraculously, diag(.6,2,2) for a and e, and diag(.1,2,2) for c are enough to get it going. Better is meanSVnv <- c(100,0). Some of the mxTryHard() iterates fail. I prefer a much smaller scale argument to mxTryHard(), such as scale =.05

Try SVN 3999

Much more helpful!

Perhaps we should edit this thread, as it never was a bug, and is not related to rows of all-missing...

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

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

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.