README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  fig.width = 8, fig.height = 8,
  cache.path = "man/cache/README-",
  out.width = "50%"
)
```
# lme4qtl

<!-- badges: start -->
[![travis-ci build status](https://travis-ci.org/variani/lme4qtl.svg?branch=master)](https://travis-ci.org/variani/lme4qtl)
<!-- badges: end -->

`lme4qtl` extends the [lme4](https://github.com/lme4/lme4) R package for quantitative trait locus (qtl) mapping. It is all about the covariance structure of random effects. `lme4qtl` supports user-defined matrices for that,
e.g. kinship or IBDs.

See slides [bit.ly/1UiTZvQ](http://bit.ly/1UiTZvQ) introducing the `lme4qtl` R package or read our [article](http://dx.doi.org/10.1186/s12859-018-2057-x) / [preprint](http://www.biorxiv.org/content/early/2017/08/31/139816).

|  Package | Continuous response |
|----------|---------------------|
| stats   | `lm(myTrait ~ myCovariate, myData)` |
| lme4    | `lmer(myTrait ~ myCovariate + (1|myID), myData)` |
| lme4qtl | `relmatLmer(myTrait ~ myCovariate + (1|myID), myData, relmat = list(myID = myMatrix))` |

|  Package | Binary response |
|----------|---------------------|
| stats    | `glm(myStatus ~ 1, myData, family = binomial)` |
| lme4    | `glmer(myStatus ~ (1|myID), myData, family = binomial)` |
| lme4qtl | `relmatGlmer(myStatus ~ (1|myID), myData, relmat = list(myID = myMatrix), family = binomial)` |

Note that rownames/colnames of `myMatrix` have to be values of `myID` variable, so matching between relationship matrix and grouping variable is possible. The order [doesn't matter](https://github.com/variani/lme4qtl/issues/3#issuecomment-346905978).

## Installation

You can install the development version from [GitHub](https://github.com/variani/lme4qtl) with:  

``` r
# install.packages("devtools")
devtools::install_github("variani/lme4qtl")
```

The official release on [CRAN](https://CRAN.R-project.org) is [pending](https://github.com/variani/lme4qtl/issues/9).

## Citation

To cite the `lme4qtl` package in publications use:

```
  Ziyatdinov et al., lme4qtl: linear mixed models with flexible
  covariance structure for genetic studies of related individuals, 
  BMC Bioinformatics (2018)
```


## Contact

You are welcome to submit suggestions and bug-reports at https://github.com/variani/lme4qtl/issues.

## Related projects

- fastGWA/GCTA tool for large-scale association studies (e.g. UKB) and based on sparse GRM; [bioarxiv](https://www.biorxiv.org/content/10.1101/598110v1) 
- [GENESIS](https://bioconductor.org/packages/release/bioc/html/GENESIS.html) Bioconductor R package

## Example

```{r ex_inc}
library(lme4)
library(lattice)
```

```{r ex_inc2, eval = FALSE}
library(lme4qtl)
```

```{r ex_inc3, echo = FALSE, message = F}
load_all()
```

```{r pkg_ver}
packageVersion("lme4qtl")
```

Load simulated data, phenotypes `dat40` and the kinship matrix `kin2`.

```{r ex_data}
data(dat40, package = "lme4qtl")

dim(dat40)
dim(kin2)

head(dat40)
kin2[1:5, 1:5] # nuclear family with 2 parents and 3 kids
```

Fit a model for continuous trait with two random effects, 
family-grouping `(1|FAM)` and additive genetic `(1|ID)`.

```{r ex_fit1}
m1 <- relmatLmer(trait1 ~ AGE + SEX + (1|FAMID) + (1|ID), dat40, relmat = list(ID = kin2))
m1
```

Get a point estimate of heritability (h2), the proportion of variance explained by `(1|ID)`.

```{r ex_h2}
lme4::VarCorr(m1)
lme4qtl::VarProp(m1)
```

Profile the variance components (h2) to get the 95% confidence intervals.
The method functions `profile` and `confint` are implemented in lme4.
Note that a different model `m2` is used,
because profiling is prone to errors/warnings if model fit is poor.

```{r ex_prof, cache = T}
m2 <- relmatLmer(trait2 ~ (1|ID), dat40, relmat = list(ID = kin2)) 
VarProp(m2)

prof <- profile(m2)
prof_prop <- lme4qtl::varpropProf(prof) # convert to proportions
confint(prof_prop)
```

```{r plot_prof, fig.width = 8, fig.height = 4}
densityplot(prof)
densityplot(prof_prop)
```

```{r plot_prof_splom}
try(splom(prof)) 

prof_clean <- na.omit(prof) # caution: NAs are indicators of poor fits
splom(prof_clean)
```

Fit a model with genetic and residual variances that differ by gender (sex-specificity model).
The formula syntax with `dummy` (see `?lme4::dummy`) is applied to the residual variance `(1|RID)`
to cancel the residual correlation. 

```{r gxe, cache = T}
dat40 <- within(dat40, RID <- ID) # replicate ID 

m4 <- relmatLmer(trait2 ~ SEX + (0 + SEX|ID) + (0 + dummy(SEX)|RID), dat40, relmat = list(ID = kin2)) 
VarCorr(m4)
```

An example of parameter constraints that make the genetic variance between genders equal.

```{r gex_eq, cache = T}
m4_vareq <- relmatLmer(trait2 ~ SEX + (0 + SEX|ID) + (0 + dummy(SEX)|RID), dat40, relmat = list(ID = kin2),
  vcControl = list(vareq = list(id = c(1, 2, 3)))) 
VarCorr(m4_vareq)
```

Another example of parameter constraint that implies the genetic correlation between genders equal to 1.

```{r gxe_rho1, cache = T}
m4_rhog1 <- relmatLmer(trait2 ~ SEX + (0 + SEX|ID) + (0 + dummy(SEX)|RID), dat40, relmat = list(ID = kin2),
  vcControl = list(rho1 = list(id = 3))) 
VarCorr(m4_rhog1)
```

Fit a model for binary trait.

```{r ex_fit3, cache = T}
m3 <- relmatGlmer(trait1bin ~ (1|ID), dat40, relmat = list(ID = kin2), family = binomial(probit))
m3
```