---
title: "Single Factor Residuals Genome Wide Association Study"
output:
  rmarkdown::html_vignette:
    self_contained: false
vignette: >
  %\VignetteIndexEntry{Single Factor Residuals GWAS}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

The Residuals GWAS Model is a variant of the One Factor GWAS Model where the SNP directly predicts one or more of the item-specific residuals, after accounting for the latent variable.  The Figure below presents a schematic representation of the Residuals Model. The identification requirements are the same as for the one factor GW-SEM model. The notable difference is that the individual items are regressed on each SNP (&gamma;<sub>k</sub>). The residuals model allows researchers to regress a SNP on the residuals of the items that contribute to a latent factor. It is also possible to fit the residuals model so that the latent factor and a subset of the items are regressed on the SNPs. If you do this, in order to identify the SEM it is only possible to regress a subset of the items and the latent factor on the SNP or the model will not be identified. 

In the Figure below, the circle (F<sub>1</sub>) is the latent factor, and the squares (x<sub>k</sub>) are observed indicators of the latent factor. The latent factor and observed indicators are connected by the factor loadings (λ<sub>k</sub>). The residual variances (δ<sub>k</sub>) represent the remaining variance in (x<sub>k</sub>) that is not shared with the latent factor. The regression of the observed items on the SNP (for all SNPs in the analysis) is depicted by (γ<sub>k</sub>). 

![A schematic depiction of the Residuals Model](https://jpritikin.github.io/gwsem/exampleFigs/resid.jpg)

## Example

To demonstrate how to use the Residuals model function, we use the same simulated data that we used for the one factor GWAS model example. Specifically, we simulated data for 6,000 individuals consisting of 2,000 SNPs, three correlated items and 6 covariates (that are a proxy for age, sex, ancestry principle components or other confounds). The simulated data are loosely based upon the data used in the demonstration of Pritikin et al (under review). In the paper, we examined the association a latent substance use frequency variable that used tobacco, cannabis, and alcohol as indicators. The latent variable was regressed on each SNP with a minor allele frequency larger than 1%.

## Running the One-Factor Residuals GWAS

**Note: Many of these steps are similar to those that you would take to run the One Factor GWAS.**

The first step is opening R and loading the GW-SEM package into the R computing environment (which will also load all of the dependencies, such as OpenMx), which can be achieved by using the command below:

```{r}
library(gwsem)
```

The next step is to load the phenotypic data. This will include any items that define the latent variable or necessary covariates that you will include in the analysis, such as age, sex, or ancestry principle components. Once the phenotypic data has been loaded, it is advisable to make sure the data has been loaded properly. This can be done in a number of ways. For example, you can use the `head` command to inspect the first few lines of data, look at the means and variances, construct tables of ordinal variables, etc.  If your data didn't load properly, it will likely cause problems for the analysis and for the interpretation of the results.

```{r}
location <- 'https://jpritikin.github.io/gwsem/gwsemOneFacExample'
phenoData <- read.table(file.path(location, "oneFacphenoData.txt"), header=TRUE)
head(phenoData)
```

Once the data are loaded into R, you can recode the data, transform it, and tell R if we have ordinal or binary indicators using `mxFactor()`. The data for this example are simulated and continuous, and therefore, we will not be doing anything now, but if you would like to chop the indicators up into binary or ordinal variables, this would be when you would do it.


After you are satisfied that the data are in the appropriate shape we can build a one factor GWAS model with the following command: 

```{r}
                                                                                       # You must tell GW-SEM:
addFacRes <- buildOneFacRes(phenoData,                                                 # what the data object is (which you read in above)
                     itemNames = c("tobacco", "cannabis", "alcohol"),                  # what the items of the latent factor are
                     res = c("tobacco", "cannabis", "alcohol"),                        # the items you want to regress onto the SNP (the default is all of the items)
                     factor = F,                                                       # factor is set to F by default
                     covariates=c('pc1','pc2','pc3','pc4','pc5'),                      # what covariates that you want to include in the analysis
                     fitfun = "WLS")                                                   # and the fit function that you would like to use (WLS is much faster than ML)

```

**Note about regressing the latent factor and SOME of the items on the SNPs:** The model is not identified if both the latent variable and all of the items are simultaneously regressed on the SNP. In some cases, however, it may be important to understand the effect of the SNP on both the latent factor and the items. In such cases, you must remove at least path from the SNP to the item residuals, and then it is possible to toggle the _factor_ argument from FALSE to TRUE.

You can take this object (`addFacRes`) which is technically an OpenMx model, and simply fit it using `mxRun()`

```{r}
addFacResFit <- mxRun(addFacRes)
summary(addFacResFit)
```
This is strongly advised, as it is a great time to test whether the model is being specified the way you want it to be and that you are not getting unexpected estimates.

Provided that the model looks reasonable, you can plug the model that you have built into the `GWAS` function using the command below:

```{r}
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'))

GWAS(model = addFacRes,                                                                 # what model object you would like to fit
	snpData = file.path(tempdir(), 'example.pgen'),                                                        # that path to the snpData file.
	out=file.path(tempdir(), "FacRes.log"),                                                                # the file that you would like to save the full results into
	SNP=1:200)                                                                       # the index of the snps (how many) you would like to fit
```

While the `GWAS` function will take a while and frequently be executed on a computing cluster, it is very useful to run a few SNPs (say 10 or 50) in an interactive session to ensure that: all of your relative file paths to the genotypes are correct, the model is taking a reasonable amount of time (i.e. 1-2 seconds/snp), the SNPs are giving sensible estimates, etc., as the results from a few SNPs can often tell you if there is a problem and that you are running a nonsensical model.

Refer to local documentation to learn how to submit jobs to your cluster queue submission systems. The R script that you provide will typically be similar to the one below:

```{r}
library(gwsem)

phenoData <- read.table(file.path(location, "oneFacphenoData.txt"), header=TRUE)
head(phenoData)

addFacRes <- buildOneFacRes(phenoData,                                                 # what the data object is (which you read in above)
                     c("tobacco", "cannabis", "alcohol"),                              # what the items of the latent factor are
                     covariates=c('pc1','pc2','pc3','pc4','pc5'),                      # what covariates that you want to include in the analysis
                     fitfun = "WLS")                                                   # and the fit function that you would like to use (WLS is much faster than ML)

GWAS(model = addFacRes,                                                                 # what model object you would like to fit
	snpData = file.path(tempdir(), 'example.pgen'),                                                        # that path to the snpData file.
	out=file.path(tempdir(), "facRes.log"))                                                                # the file that you would like to save the full results into
```

The next step is to read the results into R. The output from the `GWAS` function contains all of the estimates and standard errors for all parameters in the model (for each SNP), as well as other general model fitting information.  While you are unlikely to want to read all of this output into R, it is useful to do this with a small test GWAS analysis to ensure that you are getting reasonable results and to familiarize yourself with what information is available, in case you need to explore a specific parameter or model in more detail. This can be done with,

```{r}
FullResult <- read.delim(file.path(tempdir(), "facRes.log"))
```

More likely, you are going to want to focus on only the SNP regression coefficients. The `loadResults` function takes takes two arguments: the path to the data and the column name from the results file for the parameter that you want to examine.

Now that we have run the residual's model, we have three parameters that we want to load into R. The function to do this is the same as for one parameter, but you need to do it three times. This obtains several R objects,

```{r}
succinctTob <- loadResults(path = file.path(tempdir(), "facRes.log"), focus =  "snp_to_tobacco")
succinctCan <- loadResults(path = file.path(tempdir(), "facRes.log"), focus =  "snp_to_cannabis")
succinctAlc <- loadResults(path = file.path(tempdir(), "facRes.log"), focus =  "snp_to_alcohol")
```

# Now we can plot all the residual manhattan plots

We can construct a Manhattan plot using the GW-SEM `plot()` function. 

```{r}
plot(succinctTob)
plot(succinctCan)
plot(succinctAlc)
```

The `loadResults()` function formats the summary statistics from the analyses so that the data object can be seamlessly used by `qqman`, LDhub, and other post GWAS processing software.
