logistic.mm.aireml: Logistic mixed model fitting with Penalized Quasi-Likelihood...
In gaston: Genetic Data Handling (QC, GRM, LD, PCA) & Linear Mixed Models

logistic.mm.aireml

R Documentation

Logistic mixed model fitting with Penalized Quasi-Likelihood / AIREML

Description

Estimate the parameters of a logistic linear mixed model using the Penalized Quasi-Likelihood with an AIREML step for the linear model.

Usage

logistic.mm.aireml(Y, X = matrix(1, nrow = length(Y)), K, 
                   min_tau, tau, beta, constraint = TRUE, max.iter = 50L, eps = 1e-5,
                   verbose = getOption("gaston.verbose",TRUE), get.P = FALSE, EM = FALSE)

Arguments

`Y`	Binary phenotype vector
`X`	Covariable matrix. By default, a column of ones to include an intercept in the model
`K`	A positive definite matrix or a `list` of such matrices
`min_tau`	Minimal value for model parameter tau (if missing, will be set to 1e-6)
`tau`	(Optional) Optimization starting point for variance component(s) `tau`
`beta`	(Optional) Optimization starting point for fixed effect(s) `beta`
`constraint`	If `TRUE`, the model parameters respect the contraints given by `min_tau`
`max.iter`	Maximum number of iterations
`eps`	The algorithm stops when the gradient norm is lower than this parameter
`verbose`	If `TRUE`, display information on the algorithm progress
`get.P`	If `TRUE`, the function sends back the last matrix P computed in the optimization process
`EM`	If `TRUE`, the AIREML step is replaced by an EM step

Details

Estimate the parameters of the following logistic mixed model:

\mbox{logit}(P[Y=1|X,ω_1,…,ω_k]) = Xβ + ω_1 + … + ω_k

logit P(Y=1|X,omega_1,...,omega_k) = X beta + omega_1 + ... + omega_k with omega_i ~ N(0, tau_i K_i) for i \in 1, …,k .

The estimation is based on the Penalized Quasi-Likelihood with an AIREML step for the linear model (the algorithm is similar to the algorithm described in Chen et al 2016). If EM = TRUE the AIREML step is replaced by an EM step. In this case the convergence will be much slower, you're advised to use a large value of max.iter.

The variance matrices K_1, ..., K_k, are specified through the parameter K.

After convergence, the function also compute Best Linear Unbiased Predictors (BLUPs) for beta and omega.

Value

A named list with members:

`tau`	Estimate(s) of the model parameter(s) tau_1, ..., tau_k
`niter`	Number of iterations done
`P`	Last computed value of matrix P (see reference)
`BLUP_omega`	BLUPs of random effects
`BLUP_beta`	BLUPs of fixed effects beta
`varbeta`	Variance matrix for beta estimates

If get.P = TRUE, there is an additional member:

`P`	The last matrix P computed in the AIREML step

References

Gilmour, A. R., Thompson, R., & Cullis, B. R. (1995), Average information REML: an efficient algorithm for variance parameter estimation in linear mixed models, Biometrics, 1440-1450

Chen, Han et al. (2016), Control for Population Structure and Relatedness for Binary Traits in Genetic Association Studies via Logistic Mixed Models, The American Journal of Human Genetics, 653–666

Examples

# Load data
data(AGT)
x <- as.bed.matrix(AGT.gen, AGT.fam, AGT.bim)

# Compute Genetic Relationship Matrix
standardize(x) <- "p"
K <- GRM(x)

# Simulate a quantitative genotype under the LMM
set.seed(1)
mu <- 1 + x %*% rnorm(ncol(x), sd = 2)/sqrt(ncol(x))
pi <- 1/(1+exp(-mu))
y <- 1*( runif(length(pi))<pi )

# Estimates
estimates <- logistic.mm.aireml(y, K = K, verbose = FALSE)
str(estimates)

gaston documentation built on Dec. 28, 2022, 1:30 a.m.