################################################
##        Introduction to OpenMx              ##
################################################
 
# Remember to change the address for your folder
setwd('/home/elizabeth/2020/IntroToOpenMx')
# First model - 1 Individual
# Going back to our Twins2a object made at 8am
Twins2a<- read.table(file='Twins2a.csv',header=T,sep=",")

# Running the regression using lm()
vA1conFit1 <- lm(vA1con_1 ~ AGE_1, data=Twins2a)
summary(vA1conFit1) # show results
coefficients(vA1conFit1) # model coefficients

plot(vA1con_1 ~ AGE_1, data=Twins2a)
abline(vA1conFit1, col="red", lwd=4)

################################################################
# Setting up the model 
# Y = beta0 + beta1*age 
################################################################
Twins2a_1 <- subset(x=Twins2a, subset= !is.na(AGE_1) , 
                    select=c(vA1con_1, AGE_1)) 

require (OpenMx)
mxOption(NULL,"Default optimizer","NPSOL")
depVar <- 'vA1con_1'

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

# Regression betas
b0			<-mxMatrix(type="Full", nrow=1, ncol=1, free=T, values=30,
                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_1", name="Age" )

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

# Specify the data
regData      <- mxData( observed=Twins2a, 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)



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

# using subset function create new dataset without missing data
Twins2a_2 <- subset(x=Twins2a, subset= !is.na(AGE_1) , 
                    select=c(vA1con_1, vA1con_2,AGE_1, AGE_2)) 


depVar <- c('vA1con_1','vA1con_2')

# Setting up a Variance/Covariance matrix
Variance     <- mxMatrix( type="Symm", nrow=2, ncol=2, free=TRUE, 
                          values=c(15,10,15), 
                          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_1", "data.AGE_2"), name="Age" )

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

# Specify the data
regData      <- mxData( observed=Twins2a_2, 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 ) )