## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
is_CRAN <- !identical(Sys.getenv("NOT_CRAN"), "true")
knitr::opts_chunk$set(eval = !is_CRAN)

## -----------------------------------------------------------------------------
location <- 'https://jpritikin.github.io/gwsem/gwsemCustomExample'

#
# This model specification was automatically generated by Onyx
# 
library("OpenMx")                                                       # Edited
modelData <- read.table(file.path(location, "lgcData.txt"), header = TRUE)                     # Cut
manifests<-c("time1","time2","time3","time4","snp")                      # Kept
latents<-c("intercept","slope")                                          # Kept
model <- mxModel("LGC_GWAS",                                             # Edited
type="RAM",                                                              # Cut
manifestVars = manifests,                                                # Cut
latentVars = latents,                                                    # Cut
mxPath(from="intercept",to=c("time1","time2","time3","time4"), free=c(FALSE,FALSE,FALSE,FALSE), value=c(1.0,1.0,1.0,1.0) , arrows=1, label=c("intercept__time1","intercept__time2","intercept__time3","intercept__time4") ),   # Edited
mxPath(from="slope",to=c("time2","time3","time4"), free=c(FALSE,FALSE,FALSE), value=c(1.0,2.0,3.0) , arrows=1, label=c("slope__time2","slope__time3","slope__time4") ),                                  # Edited
mxPath(from="one",to=c("intercept","slope"), free=c(TRUE,TRUE), value=c(1.0,1.0) , arrows=1, label=c("const__intercept","const__slope") ),                                             # Edited
mxPath(from="snp",to=c("intercept","slope"), free=c(TRUE,TRUE), value=c(0.5,0.1) , arrows=1, label=c("snp__intercept","snp__slope") ),                                                 # Edited
mxPath(from="intercept",to=c("intercept","slope"), free=c(TRUE,TRUE), value=c(1.0,0.1) , arrows=2, label=c("sigma_i","COV_intercept_slope") ),                                               # Edited 
mxPath(from="slope",to=c("slope"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("sigma_s") ),    # Edited
mxPath(from="time1",to=c("time1"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("res") ),        # Edited 
mxPath(from="time2",to=c("time2"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("res") ),        # Edited  
mxPath(from="time3",to=c("time3"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("res") ),        # Edited 
mxPath(from="time4",to=c("time4"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("res") ),        # Edited 
mxPath(from="snp",to=c("snp"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_snp") ),        # Edited 
mxPath(from="one",to=c("time1","time2","time3","time4","snp"), free=F, value=0, arrows=1),         # Edited 
mxData(modelData, type = "raw")                     # Cut
);                                                  # Cut

result <- mxRun(model)                              # Cut
summary(result)                                     # Cut


## -----------------------------------------------------------------------------
# This model that we will build on for GW-SEM

library(gwsem)

# load the LGC phenotypic data
lgcData <- read.table(file.path(location, "lgcData.txt"), header = TRUE)

# Use the Onyx Mode to edit the necessary components for the GWAS model
manifests<-c("time1","time2","time3","time4","snp")
latents<-c("intercept","slope")

path <- list(mxPath(from="intercept",to=c("time1","time2","time3","time4"), free=c(FALSE,FALSE,FALSE,FALSE), value=c(1.0,1.0,1.0,1.0) , arrows=1, label=c("intercept__time1","intercept__time2","intercept__time3","intercept__time4") ),
             mxPath(from="slope",to=c("time2","time3","time4"), free=c(FALSE,FALSE,FALSE), value=c(1.0,2.0,3.0) , arrows=1, label=c("slope__time2","slope__time3","slope__time4") ),
             mxPath(from="one",to=c("intercept","slope"), free=c(TRUE,TRUE), value=c(1.0,1.0) , arrows=1, label=c("const__intercept","const__slope") ),
             mxPath(from="snp",to=c("intercept","slope"), free=c(TRUE,TRUE), value=c(0.5,0.1) , arrows=1, label=c("snp__intercept","snp__slope") ),
             mxPath(from="intercept",to=c("intercept","slope"), free=c(TRUE,TRUE), value=c(1.0,0.1) , arrows=2, label=c("sigma_i","COV_intercept_slope") ),
             mxPath(from="slope",to=c("slope"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("sigma_s") ),
             mxPath(from="time1",to=c("time1"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("res") ),
             mxPath(from="time2",to=c("time2"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("res") ),
             mxPath(from="time3",to=c("time3"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("res") ),
             mxPath(from="time4",to=c("time4"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("res") ),
             mxPath(from="snp",to=c("snp"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_snp") ),
             mxPath(from="one",to=c("time1","time2","time3","time4","snp"), free=F, value=0, arrows=1))

# Package the components into a single GW-SEM Model

lgcGWAS <- mxModel("lgc", type="RAM",
        manifestVars = manifests,
        latentVars = c(latents, paste0('pc', 1:5)),
        path,
        mxExpectationRAM(M="M"),
        mxFitFunctionWLS(allContinuousMethod="marginals"),
        mxData(observed=lgcData, type="raw", minVariance=0.1, warnNPDacov=FALSE))

# Add the necessary covariates
lgcGWAS <- setupExogenousCovariates(lgcGWAS, paste0('pc', 1:5), paste0('time',1:4))

# Make sure that the model you specified runs properly and gives sensible parameter estimates
LGCtest <- mxRun(lgcGWAS)
summary(LGCtest)

library(curl)
curl_download(file.path(location, 'example.pgen'),
              file.path(tempdir(),'example.pgen'))
curl_download(file.path(location, 'example.pvar'),
              file.path(tempdir(),'example.pvar'))

# Run the GWAS Model
GWAS(lgcGWAS,
     file.path(tempdir(), "example.pgen"),
     file.path(tempdir(), "lgc.log"))


## -----------------------------------------------------------------------------
LGCresInt <- loadResults(file.path(tempdir(), 'lgc.log'), 'snp__intercept')
LGCresSlo <- loadResults(file.path(tempdir(), 'lgc.log'), 'snp__slope')


## -----------------------------------------------------------------------------
head(LGCresInt[order(LGCresInt$Z, decreasing = T),])

head(LGCresSlo[order(abs(LGCresSlo$Z), decreasing = T),])


