#Set the working directory
getwd()
setwd("???")

################################################################################
# First OpenMx model - 1 Individual
################################################################################

#Read in the data
data <-read.table("ozbmi2.txt", header=T, na.strings = "NA")
head(data)
# using subset function create new dataset without missing data
OZbmi <- subset(data, age !="NA" , select=c(bmi1, age)) 

#Running the regression using lm
BMIfit <- lm(bmi1 ~ age, data=OZbmi)
summary(BMIfit) # show results
coefficients(BMIfit) # model coefficients

plot(bmi1 ~ age, data=OZbmi)
abline(BMIfit, col="red", lwd=4)


#Set up the regression in OpenMx
#4 matrices and 3 estimated parameters


################################################################################
# Setting up the model 
# Y = beta0 + beta1*age 
################################################################################

require (OpenMx)
mxOption(NULL,"Default optimizer","NPSOL")

depVar <- 'bmi1'

# Variance/Covariance matrix
Variance     <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, 
				values=11, labels='resid', name="residualVar" )

# Regression betas
b0			<-mxMatrix(type="Full", nrow=1, ncol=1, free=T, values=22,
			labels="beta0", name="Intercept" )
b1			<-mxMatrix(type="Full", nrow=1, ncol=1, free=T, values=0,
			labels="beta1", name="bAge" )

# Independent variable
x			<-mxMatrix(type="Full", nrow=1, ncol=1, free=F,
			labels="data.age", name="Age" )

# Building the model  (y=b0+b1*Age)
expMean      <- mxAlgebra(Intercept + bAge%x%Age, name="regress")

# Specify the data
regData      <- mxData( observed=OZbmi, type="raw" )

# List objects referenced in the model
# Remember these are objects to be included - not matrix names
inclusions   <- list(Variance, b0,  b1, x, expMean)

# Finally, we call up the results of the algebras as the arguments for the expectation function. 
# The dimnames map the data to the model.  
# NOTE- The matrix names defined in mxMatrix() statements are used here NOT the objects that are used to store the matrices.

exp         <- mxExpectationNormal( covariance="residualVar", means="regress",
                                    dimnames=depVar )
#The fit function declares that the model is fit using maximum likelihood.
# When combined with raw data this means full information maximum likelihood (FIML) is optimized.

funML       <- mxFitFunctionML()

#Build the model 
#specify the name of the model, the objects referenced, the data and the objective
regModel     <- mxModel( "Regression101", inclusions, regData, exp, funML )
regFit       <- mxRun( regModel, intervals=FALSE)

( regSum    <- summary( regFit ) )

# looking at the output in more depth
regFit

( regVerbose    <- summary( regFit, verbose=TRUE ) )


 
# ReRunning to look at optimization
regModel     <- mxOption(regModel,"Checkpoint Units", "iterations")
regModel     <- mxOption(regModel,"Checkpoint Count", 1)
regFit       <- mxRun( regModel, intervals=FALSE, checkpoint=T )

 
 
#####################################################################################
# SUBMODELS
#####################################################################################

# Go back and pickup the model so that we can run significance tests
# check for sig effects of age
  ageEfModel    <- regFit
  ageEfModel    <- omxSetParameters( ageEfModel, label="beta1", 
                                   free=FALSE, values=0 )
  ageEfFit      <- mxRun(ageEfModel, intervals=FALSE)
  (ageEfSumm     <- summary(ageEfFit))

# difference in fit
(deltaLL <-ageEfSumm$Minus2LogLikelihood - regSum$Minus2LogLikelihood)
# difference in df
(deltaDF <-ageEfSumm$degreesOfFreedom - regSum$degreesOfFreedom)
# significance test
(pchisq(deltaLL, lower.tail=F, deltaDF))
#R2
(R2<-ageEfFit$residualVar@values-regFit$residualVar@values)


m<-(mean(OZbmi$bmi1, na.rm=T))
n<-length(OZbmi$bmi1)-133
popV<-sum((na.omit(OZbmi$bmi1)-m)^2)/n 
sampleV<-sum((na.omit(OZbmi$bmi1)-m)^2)/(n -1)
builtin<-var(OZbmi$bmi1,na.rm=T)
openMx<-ageEfFit$residualVar@values
v<-(c(popV,sampleV,builtin,openMx))
names(v)<-c("PopV","SampV","R","OpenMx")

################################################################################
# Second OpenMx model - Add a cotwin
################################################################################

# using subset function create new dataset without missing data
OZbmi2 <- subset(data, age !="NA" , select=c(bmi1,bmi2,age)) 


depVar <- c('bmi1','bmi2')

# Variance/Covariance matrix
Variance     <- mxMatrix( type="Symm", nrow=2, ncol=2, free=TRUE, 
				values=c(11,5,11), labels=c('varA','covAB','varB'), name="residualVar" )

# Regression betas
b0			<-mxMatrix(type="Full", nrow=1, ncol=2, free=T, values=22,
			labels=c("beta0A","beta0B"), name="Intercept" )
b1			<-mxMatrix(type="Full", nrow=1, ncol=1, free=T, values=0,
			labels="beta1", name="bAge" )

# Independent variable
x			<-mxMatrix(type="Full", nrow=1, ncol=2, free=F,
			labels=c("data.age","data.age"), name="Age" )

# Building the model  (y=b0+b1*Age)
expMean      <- mxAlgebra(Intercept + bAge%x%Age, name="regress")

# Specify the data
regData      <- mxData( observed=OZbmi2, type="raw" )

# List objects referenced in the model
# Remember these are objects to be included - not matrix names
inclusions   <- list(Variance, b0,  b1, x, expMean)

# Finally, we call up the results of the algebras as the arguments for the expectation function. 
# The dimnames map the data to the model.  
# NOTE- The matrix names defined in mxMatrix() statements are used here NOT the objects that are used to store the matrices.

exp         <- mxExpectationNormal( covariance="residualVar", means="regress",
                                    dimnames=depVar )
#The fit function declares that the model is fit using maximum likelihood.
# When combined with raw data this means full information maximum likelihood (FIML) is optimized.

funML       <- mxFitFunctionML()

#Build the model 
#specify the name of the model, the objects referenced, the data and the objective
regModel2     <- mxModel( "Regression201", inclusions, regData, exp, funML )
regFit2       <- mxRun( regModel2, intervals=FALSE)

( regSum2    <- summary( regFit2 ) )