# Pscyhometric model, code by Michel Nivard
#

# code by Michel Nivard
#

# specified as described in Neale & Cardon (2004)
#

# questions mailto: m.g.nivard@vu.nl
#


dataDZ <- read.table("dataDZ.txt")
dataMZ <- read.table("dataMZ.txt")

colnames(dataMZ) <- c("fam","zyg","m1","f1","m2","f2")
colnames(dataDZ) <- c("fam","zyg","m1","f1","m2","f2")


SelVarsMZ=colnames(dataMZ[3:6])

SelVarsDZ=colnames(dataDZ[3:6])

#
#    Below you will find an example of the Psychometric multi rater model as described in Neale & Cardon (2003)
#    its path diagram and matrix specification is found on page 317 of the book.
#
require(OpenMx)


MxdataMZ <- mxData( observed=dataMZ, type="raw")
MxdataDZ <- mxData( observed=dataDZ, type="raw")


### Here we specify the means of the observed variables ###
MeansMZ <- mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=.0, 
                     label=c("m1","m1","m1", "m1"),name="MMZ")
MeansDZ <- mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=.0, 
                     label=c("m1","m1","m1", "m1"),name="MDZ")


### Here we specify the a,c and e path for the rater agreement ###
a <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="a", name="X")
c <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="c", name="Y")
e <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="e", name="Z")

### Here we square the a,c and e path for the rater agreement to obtain the varience of A,C and E  ###
VA <- mxAlgebra(X * t(X), name="A")
VC <- mxAlgebra(Y * t(Y), name="C")
VE <- mxAlgebra(Z * t(Z), name="E")  


### Here we specify the loadings of the latent "agreement" varaible on the observed variables ###
Lambda <- mxMatrix(type="Full", nrow=4, ncol=2, free=F,c(1,0
                                                        ,1,0
                                                        ,0,1
                                                        ,0,1), byrow=T, name="Ly")

### Here we create the latent covariance matrix for MZ twins ####
PsyMZ <-  mxAlgebra(rbind(cbind(A+C+E   , A+C),
                          cbind(A+C     , A+C+E)), name="PsyMZ")

### Here we create the latent covariance matrix for DZ twins ####
PsyDZ <-  mxAlgebra(rbind(cbind( A+C+E , (.5%x%A)+C),
                          cbind((.5%x%A)+C        , A+C+E)), name="PsyDZ")




### Here we specify the a,c and e paths for the rater DISagreement for the mother ratings  ###
am <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="am", name="Xm")
cm <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="cm", name="Ym")
em <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="em", name="Zm")

### Here we square the a,c and e path for the rater DISagreement to obtain the varience of A,C and E for the mother ratings  ###
VAm <- mxAlgebra(Xm * t(Xm), name="Am")
VCm <- mxAlgebra(Ym * t(Ym), name="Cm")
VEm <- mxAlgebra(Zm * t(Zm), name="Em")



### Here we specify the a,c and e paths for the rater DISagreement for the father ratings  ###
af <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="af", name="Xf")
cf <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="cf", name="Yf")
ef <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.6, label="ef", name="Zf")

### Here we square the a,c and e path for the rater DISagreement to obtain the varience of A,C and E for the father ratings  ###
VAf <- mxAlgebra(Xf * t(Xf), name="Af")
VCf <- mxAlgebra(Yf * t(Yf), name="Cf")
VEf <- mxAlgebra(Zf * t(Zf), name="Ef")  


### Here we specify the residual matrix for MZ twins ###
ThetaMZ <- mxAlgebra(rbind(cbind(Am+Cm+Em,0       ,Am+Cm    ,0),
                           cbind(0       ,Af+Cf+Ef,0        ,Af+Cf),
                           cbind(Am+Cm   ,0       ,Am+Cm+Em ,0),
                           cbind(0       ,Af+Cf   ,0        ,Af+Cf+Ef)), name="ThetaMZ")
  
### Here we specify the residual matrix for DZ twins ###  
ThetaDZ <- mxAlgebra(rbind(cbind(Am+Cm+Em   ,0            ,.5%x%Am+Cm    ,0),
                           cbind(0          ,Af+Cf+Ef     ,0             ,.5%x%Af+Cf),
                           cbind(.5%x%Am+Cm ,0            ,Am+Cm+Em      ,0),
                           cbind(0          ,.5%x%Af+Cf   ,0             ,Af+Cf+Ef)), name="ThetaDZ")   
  
  
  
### Here we specify the total expected covariance, both agreement and disagreement for MZ's and then DZ's ###  
CovMZ <- mxAlgebra( expression=Ly%*%PsyMZ%*%t(Ly)+ThetaMZ,name="CovMZ")
CovDZ <- mxAlgebra( expression=Ly%*%PsyDZ%*%t(Ly)+ThetaDZ,name="CovDZ")

### ere we specify the FIML Objectives for the MZ and DZ model ###
MZObj <- mxFIMLObjective( covariance="CovMZ", means="MMZ",
                          dimnames=SelVarsMZ) 

DZObj <- mxFIMLObjective( covariance="CovDZ", means="MDZ",
                          dimnames=SelVarsDZ) 

### Here we combine the Algebra's Matrices and Objectives for the MZ and DZ models ###
MZModel <-  mxModel("MZModel",PsyMZ,MeansMZ, MxdataMZ, CovMZ, MZObj, Lambda,ThetaMZ,a,c,e,VA,VC,VE,am,cm,em,VAm,VCm,VEm,af,cf,ef,VAf,VCf,VEf)
DZModel <-  mxModel("DZModel",PsyDZ,MeansDZ, MxdataDZ, CovDZ, DZObj, Lambda,ThetaDZ,a,c,e,VA,VC,VE,am,cm,em,VAm,VCm,VEm,af,cf,ef,VAf,VCf,VEf)

### Here we combine the Objectives of the MZ and DZ part of the Model ###
CombAlg    <- mxAlgebra(MZModel.objective + DZModel.objective, name="minus2loglikelihood")
CombAlgObj <- mxAlgebraObjective("minus2loglikelihood")

### Build the final model ###
Psychometric <- mxModel("PsychometricModel",MZModel,DZModel,CombAlg,CombAlgObj)


### Run thew model ###
FittedPsych <- mxRun(Psychometric)


### Summerize the results ###
summary(FittedPsych)