require(psych)
require(OpenMx)
source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R")

# Make a little labeling function
labFun <- function(name="matrix",nrow=1,ncol=1,lower=FALSE)
            {
                matlab <- matrix(paste(rep(name, each=nrow*ncol), rep(rep(1:nrow),ncol), rep(1:ncol,each=nrow),sep="_"),nrow,ncol)
                if (lower) {for (i in 1:(nrow-1)){for (j in (i+1):ncol){matlab[i,j]<-NA}}}
                return(matlab)
            }


# Import the data into R

data <- read.csv("conGrowth.csv")


head(data) 
describe(data)
table(data$zyg)


# Report the correlations within the ages for BPD and the other phenotypes

cor(c(data$suc12_1,data$suc12_2), c(data$suc13_1,data$suc13_2), use = "complete")
cor(c(data$suc13_1,data$suc13_2), c(data$suc14_1,data$suc14_2), use = "complete")
cor(c(data$suc14_1,data$suc14_2), c(data$suc15_1,data$suc15_2), use = "complete")
cor(c(data$suc15_1,data$suc15_2), c(data$suc16_1,data$suc16_2), use = "complete")
cor(c(data$suc16_1,data$suc16_2), c(data$suc17_1,data$suc17_2), use = "complete")

## Plotting Conduct Disorder Symptom Counts
error.bar.se <- function(x, y, upper, lower=upper, length=0.1,...){
if(length(x) != length(y) | length(y) !=length(lower) | length(lower)
!= length(upper))
stop("vectors must be same length")
arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}

MEAN <- describe(data[,1:6])[3]
STDV <- describe(data[,1:6])[13]

age   <- c(12,13,14,15,16,17)

# Graph the mean BPD Symptom Counts 
plot(MEAN[,1], type="o", col="black", xaxt = "n", xlab = "Age"
, ylab = "Mean Current Substance Use", ylim=c(0,1), lwd = 2 )
error.bar.se(c(1:6),MEAN[,1], 1.96*STDV[,1], col = "grey")
axis(1, at=1:6, lab= age)
# Create a title with a bold/italic font
title(main="Means Substance Use Counts by Age", col.main="black", font.main=4)

# Derive the Contrast Codes for the Factor Loadings
orthog <- function(waves){
int <- rep(1, waves)
lin <- (1:waves)-(mean(1:waves))
qua <- ((1:waves)-(mean(1:waves)))^2 - mean(((1:waves)-(mean(1:waves)))^2)

lam <- cbind(int, lin, qua)
out <- t(lam) %*% lam

lambda <- matrix(NA, waves, 3)
for(i in 1:3){lambda[,i] <- lam[,i]/ sqrt(out[i,i])}
	lambda
}


selVars <- colnames(data)[1:12]
mzData <- subset(data, zyg==1, select = selVars)
dzData <- subset(data, zyg==2, select = selVars)

head(mzData)
head(dzData)



cor(mzData[,c(1:12)], use="complete")
cor(dzData[,c(1:12)], use="complete")


#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#
#@#             Quadratic Model                     #@#
#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#@#

nv     <- 6            # 6 measurement occations per person
nvar   <- 1            # 1 set of variables
npeop  <- 2            # 2 twins 
nf     <- 3            # Latent Intercept & Linear & Quadratic Slope
ntv    <- nv * npeop   # 6 measurement occations * 2 twins
nth    <- 3


FacLoad <- orthog(nv)




#######

	 # Matrices store a, c, and e path coefficients for latent growth factors
aa <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=T, labels = labFun("aChol", nrow = nf, ncol = nf, lower = TRUE), values=.2, name="a" ) # Gen effects of Growth Factors
cc <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=T, labels = labFun("cChol", nrow = nf, ncol = nf, lower = TRUE), values=.2, name="c" ) # Com effects of Growth Factors
ee <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=T, labels = labFun("eChol", nrow = nf, ncol = nf, lower = TRUE), values=.2, name="e" ) # Ini effects of Growth Factors 

	 # Matrices T, U, and V to store a, c, and e path coefficients for residuals of observed phenotypes
as <- mxMatrix( type="Diag",  nrow=nv, ncol=nv, free=T, values=0.3, labels = labFun("as", nrow = nv, ncol = nv), name="as" ) # Specific genetic effects on obs phenotypes
cs <- mxMatrix( type="Diag",  nrow=nv, ncol=nv, free=T, values=0.3, labels = labFun("cs", nrow = nv, ncol = nv), name="cs" ) # Specific common env effects on obs phenotypes
es <- mxMatrix( type="Diag",  nrow=nv, ncol=nv, free=T, values=0.3, labels = labFun("es", nrow = nv, ncol = nv), name="es" ) # Specific non-shared env' effects

	 # Factor loading matrix of Intercept and Slope on observed phenotypes
load <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=F, values=FacLoad, name="load" )         # Int & Slope to obs phenotypes         
It   <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="It") # Identity matrix (removes off diagonals)
SDt  <- mxAlgebra( expression=solve(sqrt(It*Rt)), name="SDt") # Standardized deviations (inverse)

LatentA <- mxAlgebra( expression= load %&% (a %*% t(a)), name="LatentA" ) # Total genetic variance/covariance Int & Slope
LatentC <- mxAlgebra( expression= load %&% (c %*% t(c)), name="LatentC" ) # Total shared env' variance/covariance Int & Slope 
LatentE <- mxAlgebra( expression= load %&% (e %*% t(e)), name="LatentE" ) # Total non-shared env'  variance/covariance Int & Slope  
ResA    <- mxAlgebra( expression= as %*% t(as), name="ResA" ) # Total genetic variance/covariance Residuals
ResC    <- mxAlgebra( expression= cs %*% t(cs), name="ResC" ) # Total shared env' variance/covariance Residuals 
ResE    <- mxAlgebra( expression= es %*% t(es), name="ResE" ) # Total non-shared env'  variance/covariance Residuals 

	 # Total A, C, and E variance components (latents factors + residuals)
AA <- mxAlgebra( expression= LatentA + ResA, name="A" ) # Total genetic variance/covariance
CC <- mxAlgebra( expression= LatentC + ResC, name="C" ) # Total shared env' variance/covariance 
EE <- mxAlgebra( expression= LatentE + ResE, name="E" ) # Total non-shared env'  variance/covariance  
Rt <- mxAlgebra( expression= A+C+E, name="Rt" ) # Total variance


	 # Slope & Intercept (co)variance + standard deviations
Al <- mxAlgebra( expression=a %*% t(a), name="Al" ) # Slope & Intercept genetic (co)variance
Cl <- mxAlgebra( expression=c %*% t(c), name="Cl" )
El <- mxAlgebra( expression=e %*% t(e), name="El" )
Rf <- mxAlgebra( expression= Al+Cl+El , name="Rf" ) # Total Slope & Intercept  (co)variance
CovProp <- mxAlgebra( expression= cbind(Al/Rf,Cl/Rf,El/Rf), name = "CovProp" )
Ident <- mxMatrix( type="Iden", nrow=nf, ncol=nf, name="Ident")
SDf <- mxAlgebra( expression=(solve(sqrt(Ident*Rf)) ), name="SDf")  # Inverse of the standard deviations 

     # Means & Betas for Intercept and Slope     
Mean <- mxMatrix( type="Full", nrow=nf, ncol=1, free=T, values=c(1,.1,.1), labels=c("Ibpd","Lbpd","Qbpd"), name="Mean" ) # Int & Slope means
expMean <-  mxAlgebra(expression=cbind( t(load %*% Mean ), t(load %*% Mean) ), name="expMean")   # Fix This so that it is only in the model once

# Calculating the thresholds

# Inc	      <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=F, values=1, name="Inc" )
# thre      <- mxMatrix( type="Full", nrow=nth, ncol=nv, free=c(F, F, T), labels=paste (c("th1","th2","th3"), rep(c("_12", "_13", "_14", "_15", "_16", "_17") , each = 3), sep = ""), values=c(0,1, .5), name="thre" )
# expThre   <- mxAlgebra( expression= cbind( Inc %*% thre, Inc %*% thre ), name="expThre")   


     # Algebra for expected variance/covariance matrice for MZs and DZs 
expCovMZ <- mxAlgebra(expression=(rbind(cbind(A+C+E, A+C),
   					                    cbind(A+C, A+C+E))), name="expCovMZ" )
expCovDZ <- mxAlgebra(expression=(rbind(cbind(A+C+E, 0.5%x%A+C),
   					                    cbind(0.5%x%A+C, A+C+E))), name="expCovDZ" )
										

pars <- list(aa, cc, ee, as, cs, es, load, It, SDt,  
             LatentA, LatentC, LatentE, ResA, ResC, ResE,
             AA, CC, EE, Rt, Al, Cl, El, Rf, CovProp, Ident, SDf, 
             Mean, expMean, expCovMZ, expCovDZ)

   		
    # Algebra for expected Int & Slope means + sex effects (betas) for MZs
MZdata <- 	mxData(mzData, type="raw" )
DZdata <- mxData(dzData, type="raw" )

# mzObj <- mxFIMLObjective( covariance="expCovMZ", means="expMean", thresholds = "expThre", dimnames=selVars ) 
# dzObj <- mxFIMLObjective( covariance="expCovDZ", means="expMean", thresholds = "expThre", dimnames=selVars )

mzObj <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars ) 
dzObj <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars )


MZ     <- mxModel("MZ", MZdata, pars, MZdata, mzObj)           
DZ     <- mxModel("DZ", DZdata, pars, DZdata, dzObj)           

sumll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="sumll" )
totObj <- mxAlgebraObjective("sumll")

lgcACEModel <- mxModel("lgc", MZ, DZ, sumll, totObj)
lgcACEFit  <- mxRun(lgcACEModel) 		# Run Latent Growth Curve ACE model

summary(lgcACEFit)



# No Correlation Between Slope and Intercept
M2 <- omxSetParameters(lgcACEFit, labels = c("aChol_2_1","aChol_3_1","aChol_3_2","cChol_2_1","cChol_3_1","cChol_3_2","eChol_2_1","eChol_3_1","eChol_3_2"), newlabels = c("cor"), values = 0, free= FALSE, name = "NoGrR")
F2 <- mxRun(M2);S2 <- summary(F2)

# No General A
M3 <- omxSetParameters(lgcACEFit, labels = c("aChol_1_1","aChol_2_1","aChol_3_1","aChol_2_2","aChol_3_2"), newlabels = c("cor"), values = 0, free= FALSE, name = "NoGenA")
F3 <- mxRun(M3); S3 <- summary(F3)

# No General C
M4 <- omxSetParameters(lgcACEFit, labels = c("cChol_1_1","cChol_2_1","cChol_3_1","cChol_2_2","cChol_3_2"), newlabels = c("cor"), values = 0, free= FALSE, name = "GenC")
F4 <- mxRun(M4);S4 <- summary(F4)

# No General A or C
M5 <- omxSetParameters(lgcACEFit, labels = c("aChol_1_1","aChol_2_1","aChol_3_1","aChol_2_2","aChol_3_2","cChol_1_1","cChol_2_1","cChol_3_1","cChol_2_2","cChol_3_2"), newlabels = c("cor"), values = 0, free= FALSE, name = "GenAC")
F5 <- mxRun(M5);S5 <- summary(F5)

# No General E
M6 <- omxSetParameters(lgcACEFit, labels = c("eChol_1_1","eChol_2_1","eChol_3_1","eChol_2_2","eChol_3_2"), newlabels = c("cor"), values = 0, free= FALSE, name = "GenE")
F6 <- mxRun(M6);S6 <- summary(F6)

# No Specific A
M7 <- omxSetParameters(lgcACEFit, labels = c("as_1_1","as_2_2","as_3_3","as_4_4","as_5_5","as_6_6"), newlabels = c("cor"), values = 0, free= FALSE, name = "NoSpA")
F7 <- mxRun(M7);S7 <- summary(F7)

# No Specific C
M8 <- omxSetParameters(lgcACEFit, labels = c("cs_1_1","cs_2_2","cs_3_3","cs_4_4","cs_5_5","cs_6_6"), newlabels = c("cor"), values = 0, free= FALSE, name = "NoSpC")
F8 <- mxRun(M8);S8 <- summary(F8)

# No Specific A or C
M9 <- omxSetParameters(lgcACEFit, labels = c("as_1_1","as_2_2","as_3_3","as_4_4","as_5_5","as_6_6","cs_1_1","cs_2_2","cs_3_3","cs_4_4","cs_5_5","cs_6_6"), newlabels = c("cor"), values = 0, free= FALSE, name = "NoSpAC")
F9 <- mxRun(M9);S9 <- summary(F9)


tableFitStatistics(lgcACEFit,list(F2,F3,F4,F5,F6,F7,F8,F9))

# Get the Proportion of the Variance due to genetic factors
lgcACEFit@output$algebras$MZ.CovProp

lgcACEFit@output$matrices$MZ.a %*% t(lgcACEFit@output$matrices$MZ.a)
lgcACEFit@output$matrices$MZ.c %*% t(lgcACEFit@output$matrices$MZ.c)
lgcACEFit@output$matrices$MZ.e %*% t(lgcACEFit@output$matrices$MZ.e)

# Genetic and Enviornmental Cholesky Matrices
F8@output$matrices$MZ.a
F8@output$matrices$MZ.c
F8@output$matrices$MZ.e

# Genetic and Enviornmental Covaraince matrices
F8@output$algebras$MZ.Al
F8@output$algebras$MZ.Cl
F8@output$algebras$MZ.El
F8@output$algebras$MZ.Rf


cov2cor(F8@output$algebras$MZ.Al)
cov2cor(F8@output$algebras$MZ.Cl)
cov2cor(F8@output$algebras$MZ.El)
cov2cor(F8@output$algebras$MZ.Rf)

F8@output$matrices$MZ.as
F8@output$matrices$MZ.cs
F8@output$matrices$MZ.es


# The Best Model for our Analysis


#Grab the Mean Growth parameter effects
means <- t(F8@output$matrices$MZ.Mean)
means

head(data)

raw <- cbind(data$suc12_1,data$suc13_1,data$suc14_1,data$suc15_1,data$suc16_1,data$suc17_1)
traj <- load@values %*% t(means)

# Graph the Effects

plot(load@values[,2],traj, ylab = "Substance Use Propensity", xlab = "Age", ylim = c(-2,5), type = "b", xaxt = "n" )
axis(side = 1, at = load@values[,2], labels = c(12,13,14,15,16,17))
for(i in 1:400){
	lines(load@values[,2],raw[sample(1:nrow(raw),1),], ylab = " ", xlab = " ", ylim = c(0,2), type = "l", xaxt = "n", col = i , lwd = .5)
}
lines(load@values[,2],traj, ylab = " ", xlab = " ", ylim = c(0,1), type = "b", xaxt = "n" , lwd = 3)