BASIC Question about running ACE model

5 replies [Last post]
mdiaspar's picture
Offline
Joined: 05/26/2011

Hello!,
I have prepared a script to compute a factor model (path diagram for Phillips and Fulker model ) of phenotypes of parents and twin children given by Neale.

I get an error message indicating the following sentence:
+ mxAlgebra( expression= (A%*%(S%*%t(C)))+(C%*%(t(S)%*%t(A))), name="J" ),

# MZ Twin
+ mxAlgebra( expression= (A%&%G%)+(C%&%R%)+J+(N%*%t(N)), name="U1" ),

Error: unexpected '&' in:
"mxAlgebra( expression= (A%*%(S%*%t(C)))+(C%*%(t(S)%*%t(A))), name="J" ),
mxAlgebra( expression= (A%&%G%)+(C%&"

> mxAlgebra( expression= H%x%(A%&%(G+H%x%(T%&%(t(D)+D)))+(C%&%R)+J+((H%x%H)%x%(N%*%t(N))), name="V" ),# DZ Twin
Error: unexpected ',' in "mxAlgebra( expression= H%x%(A%&%(G+H%x%(T%&%(t(D)+D)))+(C%&%R)+J+((H%x%H)%x%(N%*%t(N))),"
>
> # Algebra - Constraints
> #===========================================================
># Phenotypic Variance constraint
> mxAlgebra( expression= A&G+C&R+E%*%t(E)+A%*%S%*%t(C)+C%*%t(S)%*%t(A)+N%*%t(N), name="Z1" )
Error: unexpected ',' in "mxAlgebra( expression= A&G+C&R+E%*%t(E)+A%*%S%*%t(C)+C%*%t(S)%*%t(A)+N%*%t(N), name="Z1" ),"
> mxConstraint(as.mxMatrix(P)==as.mxMatrix(Z1)),
Error: unexpected ',' in "mxConstraint(as.mxMatrix(P)==as.mxMatrix(Z1)),"

========================================================= FIN
I think that the following operations are supported in mxAlgebra:
t( ), %*%, %x% and %&%.

Can the openMx team give me some advice ? Thank you in advance! :-)
Best Regards,
Maikol

mspiegel's picture
Offline
Joined: 07/31/2009
There are extra '%'

There are extra '%' characters following the 'G' and 'R' in (A%&%G%)+(C%&%R%)

mdiaspar's picture
Offline
Joined: 05/26/2011
basic question

Thank very much!

mspiegel's picture
Offline
Joined: 07/31/2009
R is case-sensitive. T(M) is

R is case-sensitive. T(M) is not the same as t(M).

mdiaspar's picture
Offline
Joined: 05/26/2011
Hi, Thank very much! I have

Hi,
Thank very much!

I have another basic question.

The rawData is generated "every run" using the multivariate normal function. However,
I get an error message indicating the following sentence:
Running ACEN
Error: The job for model 'ACEN' exited abnormally with the error message: Expected covariance matrix is not positive-definite in data row 169 at iteration 0.
> summary(twinACENFit)
Error in summary(twinACENFit) :
error in evaluating the argument 'object' in selecting a method for function 'summary': Error: object 'twinACENFit' not found

Any help will be appreciated!
Anyone could tell me what happened?

###********************** The RawData will be generated using the multivariate normal function
require(OpenMx)
require(MASS)
set.seed(200)
rs=.7
n1<-500
n2<-n1/2
n3<-n2+1
xy <- mvrnorm (n1, c(0,0), matrix(c(1,rs,rs,1),2,2))
wy <- mvrnorm (n1, c(0,0), matrix(c(1,rs,rs,1),2,2))
data<-cbind(xy,wy)
Vars<- c('lbmip1','lbmim1', 'lbmi8T1', 'lbmi8T2')
dimnames(data) <- list(NULL, Vars)
mzData <-data[1:n2,]
dzData <-data[n3:n1,]
dimnames(mzData) <- list(NULL, Vars)
dimnames(dzData) <- list(NULL, Vars)

#=====================================================================
#
#Fit ACEN Model with RawData input and Matrices style (multivariate)
#
# -----------------------------------------------------------------------
twinACENModel <- mxModel("ACEN",
#===========================================================
# Matrices & Model Parameters
#===========================================================
mxMatrix( type="Full", nrow=2, ncol=2, free=TRUE, values=c(.41,.42,.39,.38), name="D" ), # assortive mating delta paths
mxMatrix( type="Symm", nrow=2, ncol=2, free=TRUE, values=.5, name="P" ), # within person covariance (Rp)
mxMatrix( type="Lower", nrow=2, ncol=2, free=TRUE, values=c(.41,.42,.39), name="A" ), # additive genetic paths
mxMatrix( type="Lower", nrow=2, ncol=2, free=TRUE, values=c(.41,.42,.39), name="C" ), # common environment paths
mxMatrix( type="Full", nrow=2, ncol=2, free=TRUE, values=c(.41,.42,.39,.38), name="F" ), # paternal cultural transmission
mxMatrix( type="Full", nrow=2, ncol=2, free=TRUE, values=c(.41,.42,.39,.38), name="M" ), # maternal cultural transmission
mxMatrix( type="Symm", nrow=2, ncol=2, free=TRUE, values=.5, name="G" ), # additive genetic covariance (Ra)
mxMatrix( type="Full", nrow=2, ncol=2, free=TRUE, values=c(.41,.42,.39,.38), name="S" ), # A-C covariance
mxMatrix( type="Symm", nrow=2, ncol=2, free=TRUE, values=.45, name="R" ), # common environmental covariance (Rc)
mxMatrix( type="Lower", nrow=2, ncol=2, free=TRUE, values=c(.41,.42,.39), name="E" ), # specific enviromental paths
mxMatrix( type="Lower", nrow=2, ncol=2, free=F, values=0, name="N" ), # non-additive paths # drop at 0
mxMatrix( type="Full", nrow=1, ncol=1, free=F, values=.5, name="H" ), # scalar
mxMatrix( type="Iden", nrow=2, ncol=2, free=F, name="I" ), # Identity matrix
mxMatrix( type="Iden", nrow=2, ncol=2, free=F, name="B" ), # common env residual variance - Identity matrix
mxMatrix( type="Full", nrow=1, ncol=8, free=T, values=1,
label="mean", name="expMean" ), # Exp Means

# Algebra - Covariance Matrices
#===========================================================
mxAlgebra( expression= P%*%D%*%t(P), name="W" ), # Mother-Father covariance

mxAlgebra( expression= G%*%t(A), name="GA" ),
mxAlgebra( expression= S%*%t(C), name="SC" ),
mxAlgebra( expression= GA+SC, name="T" ), # Mother-Child Covariance

mxAlgebra( expression= (P%*%t(F)+ W%*%t(M))%*%t(C), name="O1" ),
mxAlgebra( expression= (I+P%*%t(D))%*%t(T)%*%(H%x%t(A)), name="O2" ),
mxAlgebra( expression= O1+ O2, name="O" ), # Father-Child Covariance

mxAlgebra( expression= (P%*%t(M)+W%*%t(F))%*%t(C), name="Q1" ),
mxAlgebra( expression= Q1+O2, name="Q" ), # Mother-Child Covariance

mxAlgebra( expression= (A%*%(S%*%t(C)))+(C%*%(t(S)%*%t(A))), name="J" ),
mxAlgebra( expression= (A%&%G), name="AG" ),
mxAlgebra( expression= C%&%R, name="CR" ),
mxAlgebra( expression= N%*%t(N), name="NN" ),
mxAlgebra( expression= AG+CR+J+NN, name="U1" ), # MZ Twin

mxAlgebra( expression= t(D)+D, name="DD" ),
mxAlgebra( expression= (T%*%DD)%*%t(T), name="TDD" ),
mxAlgebra( expression= H%x%(A%&%(G+H%x%TDD))+CR+J+((H%x%H)%x%NN), name="V1" ),
mxAlgebra( expression= ((H%x%H)%x%NN), name="V2" ),
mxAlgebra( expression= V1+CR+J+V2, name="V" ), # DZ Twin

# Algebra - Constraints
#===========================================================
mxAlgebra( expression=U1+E%*%t(E), name="Z1" ), # Phenotypic Variance constraint
mxConstraint(P==Z1),

mxAlgebra( expression=TDD%*%t(T), name="Z3" ),
mxAlgebra( expression=H%x%(G+(H%x%Z3)+I) , name="Z2" ), # Genetic constraint
mxConstraint(G==Z2),

mxAlgebra( expression=H%x%T%*%(t(M)+t(F)+D%*%P%*%t(M)+t(D)%*%P%*%t(F)) , name="Z3" ), # A-C constraint
mxConstraint(S==Z3),

mxAlgebra( expression=M%*%P%*%t(M), name="MM" ),
mxAlgebra( expression=F%*%P%*%t(F), name="FF" ),
mxAlgebra( expression=(M%*%W%*%t(F))+(F%*%t(W)%*%t(M))+B , name="RES" ),
mxAlgebra( expression=MM+FF+RES , name="Z4" ), # Common Environment Constraint
mxConstraint(R==Z4),

# MZ & DZ TWINS AND PARENTS
#=================================================================
mxAlgebra( expression= rbind( cbind(P, t(W), O, O),
cbind(W, P, Q, Q),
cbind(t(O), t(Q), P, U1),
cbind(t(O), t(Q), t(U1), P)), name="expCOVMZP"),
mxAlgebra( expression= rbind( cbind(P, t(W), O, O),
cbind(W, P, Q, Q),
cbind(t(O), t(Q), P, V),
cbind(t(O), t(Q), t(V), P)), name="expCOVDZP"),

# Objective function
#===============================================================
mxModel("MZ",
mxData( observed=mzData, type="raw" ),
mxFIMLObjective( covariance="ACEN.expCOVMZP", means="ACEN.expMean", dimnames=c('lbmip1', 'lbmip1', 'lbmim1', 'lbmim1', 'lbmi8T1', 'lbmi8T1', 'lbmi8T2','lbmi8T2') ) ),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACEN.expCOVDZP",means="ACEN.expMean", dimnames=c('lbmip1', 'lbmip1', 'lbmim1', 'lbmim1', 'lbmi8T1', 'lbmi8T1', 'lbmi8T2','lbmi8T2') ) ),
mxAlgebra( expression=MZ.objective + DZ.objective, name="twin" ),
mxAlgebraObjective("twin")
)
#===========================================================
# run THE ACEN MODEL
#==========================================================

twinACENFit <- mxRun(twinACENModel)
summary(twinACENFit)

mspiegel's picture
Offline
Joined: 07/31/2009
See