Bayesian Variable Selection in Generalized Linear Mixed Models

Abstract

Repeated measures and longitudinal data are commonly collected for analysis in epidemiology, clinical trials, biology, sociology, and economic sciences. In such studies, a response is measured repeatedly over time for each subject under study, and the number and timing of the observations often varies among subjects. In contrast to cross-sectional studies that collect a single measurement per subject, longitudinal studies have the extra complication of within-subject dependence in the repeated measures. Such dependence can be thought to arise due to the impact of unmeasured predictors. Main effects of unmeasured predictors lead to variation in the average level of response among subjects, while interactions with measured predictors lead to heterogeneity in the regression coefficients. This justification has motivated random effects models, which allow the intercept and slopes in a regression model to be subject-specific. Random effects models are broadly useful for modeling of dependence not only for longitudinal data but also in multicenter studies, meta analysis and functional data analysis.

Appendix

The normal linear mixed model of Laird and Ware (1982) is a special case of a GLMM having g(μ _ij ) = μ _ij = η _ij = x′ _ij β + z′ _ij ζ _i , ϕ = σ² and $b\left({\theta _{ij} } \right) = {{\eta _{ij}^2 } /2}$ In this case,$\frac{{\partial l\left({\beta,\varphi,\zeta \left| y \right.} \right)}}{{\partial \eta }}\frac{{\partial l\left({\beta,\varphi,\zeta \left| y \right.} \right)}}{{\partial \eta '}} = \left({y - \eta } \right){{\left({y - \eta } \right)^\prime }/{\sigma ^2 }} and \frac{{\partial l\left({\beta,\varphi,\zeta \left| y \right.} \right)}}{{\partial \eta \partial \eta '}} = - 1_N {{1'_N }/{\sigma ^2 }}$. Therefore we have

$$ \begin{array}{*{20}c} {B_{i,k}^{\left(1 \right)} = \left\{ {\left({{\text{y}}_i - X_i \beta } \right)^\prime {\text{Z}}_{ik} } \right\}^2 - {\text{Z'}}_{ik} {\text{Z}}_{ik} } \\ {B_{i,m,k}^{\left(2 \right)} = \left({{\text{y}}_i - X_i \beta } \right)^\prime {\text{Z}}_{im} {\text{Z'}}_{ik} \left({{\text{y}}_i - X_i \beta } \right){\text{Z'}}_{ik} {\text{Z}}_{im},} \\ \end{array} $$

where Z _ik denotes the kth column of Z _i , and

$$ L_0 = \exp \left\{ { - \frac{1} {{2\sigma ^2 }}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^{n_i } {\left({y_{ij} - {\text{x'}}_{ij} \beta } \right)^2 } } } \right\}. $$

When y _ij are 0–1 random variables, the logistic regression model can be obtained by the canonical link function $g\left({\pi _{ij} } \right) = \log \frac{{\pi _{ij} }}{{1 - \pi _{ij} }} = \eta _{ij} = x'_{ij} \beta + z'_{ij} \zeta _i,\varphi = 1$, b(θ _ij ) = log(1 + e^{η ij}) = − log(1 − π _ij ), hence $\frac{{\partial l\left({\beta,\varphi,\zeta \left| y \right.} \right)}}{{\partial \eta }}\frac{{\partial l\left({\beta,\varphi,\zeta \left| y \right.} \right)}}{{\partial \eta '}} = \left({y - \pi } \right)\left({y - \pi } \right)^\prime and \frac{{\partial ^2 l\left({\beta,\varphi,\zeta \left| y \right.} \right)}}{{\partial \eta \partial \eta '}} = - \pi \left({1_N - \pi } \right)^\prime$. Then

$$ \begin{array}{*{20}c} {B_{i,k}^{\left(1 \right)} = \left\{ {\left({{\text{y}}_i - \pi _i } \right)^\prime {\text{Z}}_{ik} } \right\}^2 - \pi '_i {\text{DG}}\left({{\text{Z}}_{ik} {\text{Z'}}_{ik} } \right)\left({1_{ni} - \pi _i } \right)} \\ {B_{i,m,k}^{\left(2 \right)} = \left({{\text{y}}_i - \pi _i } \right)^\prime {\text{Z}}_{im} {\text{Z'}}_{ik} \left({{\text{y}}_i - \pi _i } \right) - \pi '_i {\text{DG}}\left({{\text{Z}}_{im} {\text{Z'}}_{ik} } \right)\left({1_{ni} - \pi _i } \right),} \\ \end{array} $$

where $\pi _i = \left({\pi _{i1}, \ldots,\pi _{in_i } } \right)^\prime$ with π _ij = exp(x′ _ij β)/ 1 + exp(x′ _ij β), and

$$ L_0 = \exp \left\{ {y_{ij} \log \frac{{\pi _{ij} }} {{1 - \pi _{ij} }} + \log \left({1 - \pi _{ij} } \right)} \right\}. $$

Similarly, when y _ij are counts with mean λ _ij , the Poisson regression model can be obtained by the canonical link function g(λ _ij ) = logλ _ij = η _ij = x′ _ij β + z′ _ij ζ _i , ϕ = 1, b (θ _ij ) = e ^ηij = λ _ij , $\frac{{\partial l\left({\beta,\varphi,\zeta \left| y \right.} \right)}}{{\partial \eta }}\frac{{\partial l\left({\beta,\varphi,\zeta \left| y \right.} \right)}}{{\partial \eta '}} = \left({y - \lambda } \right)\left({y - \lambda } \right)^\prime and \frac{{\partial ^2 l\left({\beta,\varphi,\zeta \left| y \right.} \right)}}{{\partial \eta \partial \eta '}} = - \lambda 1'_N$. Then we obtain that

$$ \begin{array}{*{20}c} {B_{i,k}^{\left(1 \right)} = \left\{ {\left({{\text{y}}_i - \lambda _i } \right)^\prime {\text{Z}}_{ik} } \right\}^2 - \lambda '_i {\text{DG}}\left({{\text{Z}}_{ik} {\text{Z'}}_{ik} } \right)1_{ni} } \\ {B_{i,m,k}^{\left(2 \right)} = \left({{\text{y}}_i - \lambda _i } \right)^\prime {\text{Z}}_{im} {\text{Z'}}_{ik} \left({{\text{y}}_i - \lambda _i } \right) - \lambda '_i {\text{DG}}\left({{\text{Z}}_{im} {\text{Z'}}_{ik} } \right)1_{ni},} \\ \end{array} $$

where $\lambda _i = \left({\lambda _{i1}, \ldots,\lambda _{in_i } } \right)^\prime$ with λ _ij = exp(x′ _ij β), and L ₀ = exp y _ij logλ _ij − λ _ij − logy _ij !.