# ---------------------------------------------------------------------
# Program: OneFactorModelMatrix_WELLBEING.R
#  Author: Marleen de Moor
#    Date: 22 02 2010
#
# Model: single group 1 factor model;
# Data: 6 observed variables (quality of life, happiness, life satisfaction, 
#       anxious depression, somatic complaints, social problems)
#       measured in 1000 individuals (1 twin per family)
#
# This script is an adaptation of the script OneFactorModelMatrix.R
#  Authors: Steve Boker/Hermine Maes -- 10 08 2009 updated & reformatted
# -----------------------------------------------------------------------

require(OpenMx)

# Prepare Data
# -----------------------------------------------------------------------
allData<-read.table("CFA_family_wellbeing.dat", header=TRUE, na.strings=-999)
cfaData<-allData[,c('qol','hap','sat','ad','soma','soc')]
colMeans(cfaData, na.rm=TRUE)
cov(cfaData[,c('qol','hap','sat','ad','soma','soc')],use="pairwise.complete.obs")
cor(cfaData[,c('qol','hap','sat','ad','soma','soc')],use="pairwise.complete.obs")

nvar<-6
nfac<-1

# Run single group 1 factor model - cov data input
# -----------------------------------------------------------------------
observedVars <- names(cfaData)
oneFactorModelcov <- mxModel("One Factor",
    mxMatrix(type="Full", nrow=nvar, ncol=nfac, values=0.2, free=TRUE, name="L"),
    mxMatrix(type="Symm", nrow=nfac, ncol=nfac, values=1, free=TRUE, name="P"),
    mxMatrix(type="Diag", nrow=nvar, ncol=nvar, values=1, free=TRUE, name="T"),
    mxAlgebra(expression=L %*% P %*% t(L) + T, name="expCov"),
    mxData(cov(cfaData, use="pairwise.complete.obs"),type="cov", numObs=1000),
    mxMLObjective(covariance="expCov", dimnames = observedVars))

oneFactorFitcov<-mxRun(oneFactorModelcov)
summary(oneFactorFitcov)