# Rater Bias Model, code by Michel Nivard
#

# code by Michel Nivard
#

# specified as desctibed 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])
require(OpenMx)


MxdataMZ <- mxData( observed=dataMZ[3:6], type="raw")
MxdataDZ <- mxData( observed=dataDZ[3:6], type="raw")

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")

### Below we estimate path loadings a,c and e used to model the twin correltions a the latent level ########
a <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.3, label="a", name="X")
c <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.2, label="c", name="Y")
e <- mxMatrix("Full", nrow=1, ncol=1, free=TRUE, values=.5, label="e", name="Z")

### Here we square the path loadings a, c and e, they then result in variance components ###
VA <- mxAlgebra(X * t(X), name="A")
VC <- mxAlgebra(Y * t(Y), name="C")
VE <- mxAlgebra(Z * t(Z), name="E")  


#### below we specify the matrix for the MZ factor loadings  ###
LambdaMZ <- mxMatrix(type="Full", nrow=4, ncol=2, free= c(F,F
                                                          ,T,F
                                                          ,F,F
                                                          ,F,T),values=c(1,0
                                                                         ,1,0
                                                                         ,0,1
                                                                         ,0,1),labels=c("lm","n"
                                                                                        ,"lf","n"
                                                                                        ,"n","lm"
                                                                                        ,"n","lf"), byrow=T, name="LyMZ")
#### below we specify the matrix for the DZ factor loadings  ###
LambdaDZ <- mxMatrix(type="Full", nrow=4, ncol=2, free= c(F,F
                                                          ,T,F
                                                          ,F,F
                                                          ,F,T),values=c(1,0
                                                                         ,1,0
                                                                         ,0,1
                                                                         ,0,1),labels=c("lm","n"
                                                                                        ,"lf","n"
                                                                                        ,"n","lm"
                                                                                        ,"n","lf"), byrow=T, name="LyDZ")


### below we specify the algebra for the covariance between the latent factors of MZ twins ###  
PsyMZ <-  mxAlgebra(rbind(cbind(A+C+E   , A+C),
                          cbind(A+C     , A+C+E)), name="PsyMZ")


### below we specify the algebra for the covariance between the latent factors of DZ twins ###
PsyDZ <-  mxAlgebra(rbind(cbind(A+C+E , .5%x%A+C),
                          cbind(.5%x%A+C        , A+C+E)), name="PsyDZ")


### below we specify the matrix containing the residuals and the male and female bias parameters for MZ twins ### 
ThetaMZ <- mxMatrix("Full", nrow=4, ncol=4,value= c( 1,0,.5,0
                                                     ,0,1,0,.5
                                                     ,.5,0,1,0
                                                     ,0,.5,0,1), free=as.logical(c( 1,0,1,0
                                                                                    ,0,1,0,1
                                                                                    ,1,0,1,0
                                                                                    ,0,1,0,1)), labels=c( "ResM1","n","BiasM","n"
                                                                                                          ,"n","ResF1","n","BiasF"
                                                                                                          ,"BiasM","n","ResM2","n"
                                                                                                          ,"n","BiasF","n","ResF2"),byrow=T, name="ResMZ")

### below we specify the matrix containing the residuals and the male and female bias parameters for DZ twins ### 
ThetaDZ <- mxMatrix("Full", nrow=4, ncol=4,value= c( 1,0,.5,0
                                                     ,0,1,0,.5
                                                     ,.5,0,1,0
                                                     ,0,.5,0,1), free=as.logical(c( 1,0,1,0
                                                                                    ,0,1,0,1
                                                                                    ,1,0,1,0
                                                                                    ,0,1,0,1)), labels=c( "ResM1","n","BiasM","n"
                                                                                                          ,"n","ResF1","n","BiasF"
                                                                                                          ,"BiasM","n","ResM2","n"
                                                                                                          ,"n","BiasF","n","ResF2"),byrow=T, name="ResDZ")
### below we specify the expected covariance between our 4 observed variables ###
CovMZ <- mxAlgebra( expression=LyMZ%*%PsyMZ%*%t(LyMZ)+ResMZ,name="CovMZ")
CovDZ <- mxAlgebra( expression=LyDZ%*%PsyDZ%*%t(LyDZ)+ResDZ,name="CovDZ")


### below we specify FIML objectives for MZ's and DZ's ###
MZObj <- mxFIMLObjective( covariance="CovMZ", means="MMZ",
                          dimnames=SelVarsMZ) 

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


### below we combine the matrices and algebras for MZ's and DZ's into  mxModels ##
MZModel <-  mxModel("MZModel",PsyMZ,MeansMZ, MxdataMZ, CovMZ, MZObj, LambdaMZ,ThetaMZ,a,c,e,VA,VC,VE)
DZModel <-  mxModel("DZModel",PsyDZ,MeansDZ, MxdataDZ, CovDZ, DZObj, LambdaDZ,ThetaDZ,a,c,e,VA,VC,VE)

CombAlg    <- mxAlgebra(MZModel.objective + DZModel.objective, name="minus2loglikelihood")
CombAlgObj <- mxAlgebraObjective("minus2loglikelihood")


### below we combine the MZ and DZ model ###
RaterBias <- mxModel("RaterBiasModel",MZModel,DZModel,CombAlg,CombAlgObj)


### below We run the model ###
FittedRB <- mxRun(RaterBias)
### below we look at our results ###
summary(FittedRB)