Description Usage Arguments Details Value Author(s) References Examples

Fit a generalized linear mixed model with a random intercept, or a random intercept and an optional random slope of time effect for longitudinal data. The covariance matrix of the random intercept is proportional to a known relationship matrix (e.g. kinship matrix in genetic association studies). Alternatively, it can be a variance components model with multiple random effects, and each component has a known relationship matrix.

1 2 3 4 |

`fixed` |
an object of class |

`data` |
a data frame or list (or object coercible by |

`kins` |
a known positive semi-definite relationship matrix (e.g. kinship matrix in genetic association studies) or a list of known positive semi-definite relationship matrices. The rownames and colnames of these matrices must at least include all samples as specified in the |

`id` |
a column in the data frame |

`random.slope` |
an optional column indicating the random slope for time effect used in a mixed effects model for longitudinal data. It must be included in the names of |

`groups` |
an optional categorical variable indicating the groups used in a heteroscedastic linear mixed model (allowing residual variances in different groups to be different). This variable must be included in the names of |

`family` |
a description of the error distribution and link function to be used in the model. This can be a character string naming a family function, a family function or the result of a call to a family function. (See |

`method` |
method of fitting the generalized linear mixed model. Either "REML" or "ML" (default = "REML"). |

`method.optim` |
optimization method of fitting the generalized linear mixed model. Either "AI", "Brent" or "Nelder-Mead" (default = "AI"). |

`maxiter` |
a positive integer specifying the maximum number of iterations when fitting the generalized linear mixed model (default = 500). |

`tol` |
a positive number specifying tolerance, the difference threshold for parameter estimates below which iterations should be stopped (default = 1e-5). |

`taumin` |
the lower bound of search space for the variance component parameter |

`taumax` |
the upper bound of search space for the variance component parameter |

`tauregion` |
the number of search intervals for the REML or ML estimate of the variance component parameter |

`verbose` |
a logical switch for printing detailed information (parameter estimates in each iteration) for testing and debugging purpose (default = FALSE). |

`...` |
additional arguments that could be passed to |

Generalized linear mixed models (GLMM) are fitted using the penalized quasi-likelihood (PQL) method proposed by Breslow and Clayton (1993). Generally, fitting a GLMM is computationally expensive, and by default we use the Average Information REML algorithm (Gilmour, Thompson and Cullis, 1995; Yang et al., 2011) to fit the model. If only one relationship matrix is specified (`kins`

is a matrix), iterations may be accelerated using the algorithm proposed by Zhou and Stephens (2012) for linear mixed models. An eigendecomposition is performed in each outer iteration and the estimate of the variance component parameter *τ* is obtained by maximizing the profiled log restricted likelihood (or likelihood) in a search space from `taumin`

to `taumax`

, equally divided into `tauregion`

intervals on the log scale, using Brent's method (1973). If `kins`

is a list of matrices and `method = "Nelder-Mead"`

, iterations are performed as a multi-dimensional maximization problem solved by Nelder and Mead's method (1965). It can be very slow, and we do not recommend using this method unless the likelihood function is badly behaved. Both Brent's method and Nelder and Mead's method are derivative-free. When the Average Information REML algorithm fails to converge, a warning message is given and the algorithm is default to derivative-free approaches: Brent's method if only one relationship matrix is specified, Nelder and Mead's method if more than one relationship matrix is specified.

For longitudinal data (with duplicated `id`

), two types of models can be applied: random intercept only models, and random intercept and random slope models. The random intercept only model is appropriate for analyzing repeated measures with no time trends, and observations for the same individual are assumed to be exchangeable. The random intercept and random slope model is appropriate for analyzing longitudinal data with individual-specific time trends (therefore, a random slope for time effect). Typically, the time effect should be included in the model as a fixed effect covariate as well. Covariances of the random intercept and the random slope are estimated.

For multiple phenotype analysis, `formula`

recognized by `lm`

, such as `cbind(y1, y2, y3) ~ x1 + x2`

, can be used in `fixed`

as fixed effects. For each matrix in `kins`

, variance components corresponding to each phenotype, as well as their covariance components, will be estimated. Currently, `family`

must be "gaussian" and `method.optim`

must be "AI".

`theta` |
a vector or a list of variance component parameter estimates. See below. For cross-sectional data, if For longitudinal data (with duplicated For longitudinal data (with duplicated For multiple phenotype analysis, |

`n.pheno` |
an integer indicating the number of phenotypes in multiple phenotype analysis (for single phenotype analysis, |

`n.groups` |
an integer indicating the number of distinct residual variance groups in heteroscedastic linear mixed models (for other models, |

`coefficients` |
a vector or a matrix for the fixed effects parameter estimates (including the intercept). |

`linear.predictors` |
a vector or a matrix for the linear predictors. |

`fitted.values` |
a vector or a matrix for the fitted mean values on the original scale. |

`Y` |
a vector or a matrix for the final working vector. |

`X` |
model matrix for the fixed effects. |

`P` |
the projection matrix with dimensions equal to the sample size multiplied by |

`residuals` |
a vector or a matrix for the residuals on the original scale. NOT rescaled by the dispersion parameter. |

`scaled.residuals` |
a vector or a matrix for the scaled residuals, calculated as the original residuals divided by the dispersion parameter (in heteroscedastic linear mixed models, corresponding residual variance estimates by each group). |

`cov` |
covariance matrix for the fixed effects (including the intercept). |

`Sigma_i` |
the inverse of the estimated covariance matrix for samples, with dimensions equal to the sample size multiplied by |

`Sigma_iX` |
Sigma_i multiplied by X. Used in |

`converged` |
a logical indicator for convergence. |

`call` |
the matched call. |

`id_include` |
a vector indicating the |

Han Chen, Matthew P. Conomos

Brent, R.P. (1973) "Chapter 4: An Algorithm with Guaranteed Convergence for Finding a Zero of a Function", Algorithms for Minimization without Derivatives, Englewood Cliffs, NJ: Prentice-Hall, ISBN 0-13-022335-2.

Breslow, N.E. and Clayton, D.G. (1993) Approximate Inference in Generalized Linear Mixed Models. Journal of the American Statistical Association 88, 9-25.

Chen, H., Wang, C., Conomos, M.P., Stilp, A.M., Li, Z., Sofer, T., Szpiro, A.A., Chen, W., Brehm, J.M., Celedón, J.C., Redline, S., Papanicolaou, G.J., Thornton, T.A., Laurie, C.C., Rice, K. and Lin, X. (2016) Control for population structure and relatedness for binary traits in genetic association studies via logistic mixed models. The American Journal of Human Genetics 98, 653-666.

Gilmour, A.R., Thompson, R. and Cullis, B.R. (1995) Average Information REML: An Efficient Algorithm for Variance Parameter Estimation in Linear Mixed Models. Biometrics 51, 1440-1450.

Nelder, J.A. and Mead, R. (1965) A simplex algorithm for function minimization. Computer Journal 7, 308-313.

Yang, J., Lee, S.H., Goddard, M.E. and Visscher, P.M. (2011) GCTA: A Tool for Genome-wide Complex Trait Analysis. The American Journal of Human Genetics 88, 76-82.

Zhou, X. and Stephens, M. (2012) Genome-wide efficient mixed-model analysis for association studies. Nature Genetics 44, 821-824.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | ```
data(example)
attach(example)
model0 <- glmmkin(disease ~ age + sex, data = pheno, kins = GRM, id = "id",
family = binomial(link = "logit"))
model0$theta
model0$coefficients
model0$cov
model1 <- glmmkin(y.repeated ~ sex, data = pheno2, kins = GRM, id = "id",
family = gaussian(link = "identity"))
model1$theta
model1$coefficients
model1$cov
model2 <- glmmkin(y.trend ~ sex + time, data = pheno2, kins = GRM, id = "id",
random.slope = "time", family = gaussian(link = "identity"))
model2$theta
model2$coefficients
model2$cov
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.