R forced closed due to non pd cov matrix

1 reply [Last post]
wuhao_osu's picture
Offline
Joined: 09/07/2010

Code as below. I think it is due to a non positive definite covariance input. When this is fixed, all is fine.

A pop -up says R for windows GUI front end has stopped working.

require(OpenMx)

varnames<-c("BIC1","SSC1","BIC2","SSC2");
dataMZ<-dataDZ<-matrix(diag(1,4),nrow=4,ncol=4,dimnames=list(varnames,varnames));
dataDZ[1,2]<-dataDZ[2,1]<-1.2;

p<-2;
pstar<-p*(p+1)/2;
diagind<-(diag(p)==1);

nMZ<-84;
nDZ<-33;
varnames<-c("BIC1","SSC1","BIC2","SSC2");

LA<-mxMatrix(type="Lower",nrow=p,ncol=p,free=T,values=0,
label=paste("a",1:pstar,sep=""),name="LA");
LC<-mxMatrix(type="Lower",nrow=p,ncol=p,free=T,values=0,
label=paste("c",1:pstar,sep=""),name="LC");
LE<-mxMatrix(type="Lower",nrow=p,ncol=p,free=T,values=c(1,.5,sqrt(3)/2),
label=paste("e",1:pstar,sep=""),name="LE");

LA@lbound[diagind]<-0;
LC@lbound[diagind]<-0;
LE@lbound[diagind]<-0;

A<-mxAlgebra(expression=LA%*%t(LA),name="A");
C<-mxAlgebra(expression=LC%*%t(LC),name="C");
E<-mxAlgebra(expression=LE%*%t(LE),name="E");

ACESigmaMZ<-mxAlgebra(expression=rbind(cbind(A+C+E,A+C),cbind(A+C,A+C+E)),
name="ACESigmaMZ");
ACESigmaDZ<-mxAlgebra(expression=rbind(cbind(A+C+E,.5%x%A+C),cbind(.5%x%A+C,A+C+E)),
name="ACESigmaDZ");

ACEModelMZ<-mxModel("ACEModelMZ",
mxMLObjective(covariance="ACEModelTwin.ACESigmaMZ",dimnames=varnames),
mxData(observed=dataMZ,type="cov",numObs=nMZ));
ACEModelDZ<-mxModel("ACEModelDZ",
mxMLObjective(covariance="ACEModelTwin.ACESigmaDZ",dimnames=varnames),
mxData(observed=dataDZ,type="cov",numObs=nDZ));
twin<-mxAlgebra(expression=ACEModelMZ.objective+ACEModelDZ.objective,name="twin");
ACEModelTwin<-mxModel("ACEModelTwin",LA,LC,LE,A,C,E,
ACESigmaMZ,ACESigmaDZ,ACEModelMZ,ACEModelDZ,twin,
mxAlgebraObjective("twin"));

mxRun(ACEModelTwin)

mspiegel's picture
Offline
Joined: 07/31/2009
This is a bug in the backend.

This is a bug in the backend. Your intuition is correct, it is caused by a non-positive-definite observed covariance matrix. However we definitely shouldn't be throwing a segmentation fault in these cases. We're working on a fix.