In variani/lme4qtl: Linear mixed models with custom covariance matrices and restrictions on model parameters.

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

lme4qtl extends the 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 introducing the lme4qtl R package or read our article / preprint.

| 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.

Installation

You can install the development version from GitHub with:

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

The official release on CRAN is pending.

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
• GENESIS Bioconductor R package

Example

library(lme4)
library(lattice)

library(lme4qtl)

load_all()

packageVersion("lme4qtl")

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

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).

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).

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.

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)

densityplot(prof)
densityplot(prof_prop)

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.

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.

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.

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.

m3 <- relmatGlmer(trait1bin ~ (1|ID), dat40, relmat = list(ID = kin2), family = binomial(probit))
m3

variani/lme4qtl documentation built on Jan. 10, 2021, 2:45 p.m.