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(IV)), 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(IV)) %% V %% solve(sqrt(IV)))[1,2],name="Vcorr"),
mxAlgebra((solve(sqrt(IA)) %% A %% solve(sqrt(IA)))[1,2],name="Acorr"),
mxAlgebra((solve(sqrt(IC)) %% C %% solve(sqrt(IC)))[1,2],name="Ccorr"),
mxAlgebra((solve(sqrt(IE)) %% E %% solve(sqrt(IE)))[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