R package wlm

The wlm R package supports linear models, which are more complex than OLS, but the output is still that as returned by lm. That means you can use the formula interface, handy print, coef, vcov and other functions. Missing values are handled the best way possible, but it is better to provide clean data.

See Fitting linear mixed models for QTL mapping blog post by Karl Broman for examplanation of the eigen-decomposition trick possible for LMM1 model (1 random genetic + 1 random noise effects).

Installation

library(devtools)
install_github("variani/wlm")

Example code

library(wlm)

# data
data(mtcars)

mtcars <- within(mtcars, {
  cyl <- factor(cyl)
  weight_cyl = 1/sqrt(as.numeric(cyl))
})
  
# OLS
m1 <- lm(mpg ~ disp + cyl, mtcars)

# WLS
m2 <- lm(mpg ~ disp, mtcars, weights = weight_cyl)

# GLS
varcov_cyl <- with(mtcars, sapply(cyl, function(x) as.numeric(x) * as.numeric(x == cyl)))

m3 <- wlm(mpg ~ disp, mtcars, varcov = varcov_cyl)

# LMM1
varcov_cyl <- with(mtcars, sapply(cyl, function(x) as.numeric(x == cyl)))

m4 <- lmm1(mpg ~ disp, mtcars, varcov = varcov_cyl)

#> m4$lmm$r2
#[1] 0.1049433
# ~10% of variance explained by the random effect with `varcov = varcov_cyl`

# variance-covariance matrix of the outcome 
var_mpg <- m4$lmm$r2 * varcov_cyl + (1 - m4$lmm$r2) * diag(nrow(varcov_cyl))

# `m4` is equivalent to the following GLS
m4_wlm <- wlm(mpg ~ disp, mtcars, varcov = varcov_mpg)

# LMM1 + precomputed eigendecomposition 
# - useful when you need to test many predictors such as in GWAS
evd_varcov_cyl <- eigen(varcov_cyl)

m5 <- lmm1(mpg ~ disp, mtcars, varcov = evd_varcov_cyl)

# Computation time of LMMs
system.time(lmm1(mpg ~ disp, mtcars, varcov = varcov_cyl))
#   user  system elapsed
#  0.005   0.000   0.005

library(lme4qtl)
system.time(lmer(mpg ~ disp + (1|cyl), mtcars))
#   user  system elapsed
#  0.014   0.000   0.014