# -----------------------------------------------------------------------
# Program: Att8P.R  
#  Author: Hermine Maes
#    Date: 03 08 2011 
#
# Multivariate model to estimate means and (co)variances
# Raw Data -Matrix Style
#
# -----------------------------------------------------------------------

#----------------------------------------------------------------------------------------------
#  Reading in Phenotypic Data 
#----------------------------------------------------------------------------------------------

mzmPData <- read.table("mzmP.dat",sep="",na.strings=".", header=T)
dzmPData <- read.table("dzmP.dat",sep="",na.strings=".", header=T)
mzfPData <- read.table("mzfP.dat",sep="",na.strings=".", header=T)
dzfPData <- read.table("dzfP.dat",sep="",na.strings=".", header=T)
dzoPData <- read.table("dzoP.dat",sep="",na.strings=".", header=T)


#----------------------------------------------------------------------------------------------
#  Start of the Model 
#----------------------------------------------------------------------------------------------
require(OpenMx)
manifestVars  <- c("mo3t1","fa3t1","mo7t1","fa7t1","mo10t1","fa10t1","mo12t1","fa12t1","mo3t2","fa3t2","mo7t2","fa7t2","mo10t2","fa10t2","mo12t2","fa12t2")
nv <- 8		# number of variables

factorModel <- mxModel("factor",

  mxModel("All",
#  Means, freely estimated
		mxMatrix( type="Full", nrow=1, ncol=nv, free=T, 
			values=0, labels=c("meanm3","meanf3","meanm7","meanf7","meanm10","meanf10","meanm12","meanf12"), name="Mean"),
#  Regressions on sex, freely estimated
		mxMatrix( type="Full", nrow=1, ncol=nv, free=T, 
			values=0.1, labels=c("beta1","beta2","beta3","beta4","beta5","beta6","beta7","beta8"), name="Beta"),
#  Factor loadings, freely estimated
		mxMatrix( type="Full", nrow=nv, ncol=1, free=T, 
			values=1.5, name="facL"),  
#  Residuals, freely estimated
		mxMatrix( type="Diag", nrow=nv, ncol=nv, free=T, 
			values=2, name="res"),
#  Correlated factor loadings and correlated residuals, freely estimated
    	mxMatrix( type="Full", nrow=1, ncol=1, free=T, 
    		values=.6, labels="corMZ", name="facCorMZ"),
    	mxMatrix( type="Full", nrow=1, ncol=1, free=T, 
    		values=.1, labels="corDZ", name="facCorDZ"),
		mxMatrix( type="Diag", nrow=nv, ncol=nv, free=T, 
			values=.6, name="resCorMZ"),
		mxMatrix( type="Diag", nrow=nv, ncol=nv, free=T, 
			values=.1, name="resCorDZ"),
#  Within-person and Across-person covariances
 		mxAlgebra( expression= facL %*% t(facL), name="facVar"),
    	mxAlgebra( expression= facL %*% facCorMZ %*% t(facL), name="facCovMZ"),
    	mxAlgebra( expression= facL %*% facCorDZ %*% t(facL), name="facCovDZ"),
		mxAlgebra( expression= res %*% t(res), name="resVar"), 
		mxAlgebra( expression= res %*% resCorMZ %*% t(res), name="resCovMZ"),
		mxAlgebra( expression= res %*% resCorDZ %*% t(res), name="resCovDZ"),
#  Expected covariance matrices
       mxAlgebra( expression= rbind(cbind(facVar + resVar,     facCovMZ + resCovMZ),
                                    cbind(facCovMZ + resCovMZ, facVar + resVar)), name="covMZ"),
       mxAlgebra( expression= rbind(cbind(facVar + resVar,     facCovDZ + resCovDZ),
                                    cbind(facCovDZ + resCovDZ, facVar + resVar)), name="covDZ"),
#  Derived matrices
       mxAlgebra( expression= facVar + resVar, name="Var" ),        # phenotypic (co)variances
       mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),          #
       mxAlgebra( expression= solve(sqrt(I * Var)), name="invSD"),  # inverse of standard deviations
       mxAlgebra( expression= invSD %*% facL, name="stFacL"),       # standardized factor loadings
       mxAlgebra( expression= invSD %*% res, name="stRes"),         # standardized residuals
       mxAlgebra( expression= invSD %&% Var, name="stVar")          # standardized (co)variances
  ),

#----------------------------------------------------------------------------------------------
  mxModel("MZM",
       mxData(mzmPData, type="raw"),
       mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex"),
       mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex %x% All.Beta, name="mean"),
       mxFIMLObjective( covariance="All.covMZ", means="mean", dimnames=manifestVars)
  ),

#----------------------------------------------------------------------------------------------
  mxModel("DZM",
       mxData(dzmPData, type="raw"),
       mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex"),
       mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex %x% All.Beta, name="mean"),
       mxFIMLObjective( covariance="All.covDZ", means="mean", dimnames=manifestVars)
  ),

#----------------------------------------------------------------------------------------------
  mxModel("MZF",
       mxData(mzfPData, type="raw"),
       mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex"),
       mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex %x% All.Beta, name="mean"),
       mxFIMLObjective( covariance="All.covMZ", means="mean", dimnames=manifestVars)
  ),

#----------------------------------------------------------------------------------------------
  mxModel("DZF",
       mxData(dzfPData, type="raw"),
       mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex"),
       mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex %x% All.Beta, name="mean"),
       mxFIMLObjective( covariance="All.covDZ", means="mean", dimnames=manifestVars)
  ),

#----------------------------------------------------------------------------------------------
  mxModel("DOS",
       mxData(dzoPData, type="raw"),
       mxMatrix( type="Full", nrow=1, ncol=2, free=F, labels=c("data.SEXt1","data.SEXt2"), name="sex"),
       mxAlgebra( expression= cbind(All.Mean,All.Mean) + sex%x%All.Beta, name="mean"),
       mxFIMLObjective( covariance="All.covDZ", means="mean", dimnames=manifestVars)
  ),

#----------------------------------------------------------------------------------------------
    mxAlgebra( MZM.objective + DZM.objective + MZF.objective + DZF.objective + DOS.objective, name="minus2sumloglikelihood"),
    mxAlgebraObjective("minus2sumloglikelihood")
)

#----------------------------------------------------------------------------------------------
#  Run the Model 
#----------------------------------------------------------------------------------------------
factorFit <- mxRun(factorModel)

#----------------------------------------------------------------------------------------------
#  Summarize/Generate Output
#----------------------------------------------------------------------------------------------
factorSumm <- summary(factorFit)
factorSumm