glmmsel is an R package for generalised linear mixedmodel (GLMM) selection. Given observations on\(m\) clusters\((\mathbf{y}_i,\mathbf{X}_i)_{i=1}^m\),where\(\mathbf{y}_i\) and\(\mathbf{X}_i\) represent the responsevector and predictor matrix for cluster\(i\),glmmsel can fit a GLMM ofthe form
\[\operatorname{E}\left[\eta(\mathbf{y}_i)\right]=\mathbf{X}_i(\boldsymbol{\beta}+\mathbf{u}_i),\] where\(\boldsymbol{\beta}\)is a sparse vector of fixed effects (i.e., predictor effects that arethe same across clusters),\(\mathbf{u}_i\) is a sparse vector of randomeffects (i.e., predictor effects that differ across clusters), and\(\eta\) is a link function.glmmsel fits this model by solving the optimisation problem\[\underset{\boldsymbol{\beta},\boldsymbol{\gamma}}{\min}\;l(\mathbf{y},\mathbf{X};\boldsymbol{\beta},\boldsymbol{\gamma})+\lambda\alpha\|\boldsymbol{\beta}\|_0+\lambda(1-\alpha)\|\boldsymbol{\gamma}\|_0\quad\operatorname{s.t.}\;\beta_k=0\Rightarrow\gamma_k=0,\] where\(l\) is a negativelog-likelihood,\(\|\cdot\|_0\) is the\(\ell_0\)-norm (i.e., a count of thenumber of nonzeros), and\(\lambda\geq0\) and\(\alpha\in(0,1]\) are tuning parameters.Here,\(\boldsymbol{\gamma}\)characterises the variance of the random effects\(\mathbf{u}_i\), which we assume follow a\(N(\mathbf{0},\operatorname{diag}(\boldsymbol{\gamma}))\)distribution. Observe that if\(\gamma_k=0\) then\(u_{ik}\) is zero.
glmmsel operates on the hierarchy principle that arandom effect can only be selected if its corresponding fixed effect isalso selected; see the constraint\(\beta_k=0\Rightarrow\gamma_k=0\). Setting\(\alpha=1\) means there is no penaltyfor selecting a random effect if its fixed effect is also selected.Smaller values of\(\alpha\) encouragethe random effect to be selected only if it substantially improves thefit. The default value of\(\alpha=0.8\) works well in practice.
The two main functions provided by the package areglmmsel() andcv.glmmsel(), responsible formodel fitting and cross-validation, respectively.
Theglmmsel() function provides a convenient way offitting the model for a path of\(\lambda\) values. To demonstrate thisfunctionality, let’s simulate some clustered data.
set.seed(1234)n<-100# Number of observationsm<-4# Number of clustersp<-5# Number of predictorss.fix<-2# Number of nonzero fixed effectss.rand<-1# Number of nonzero random effectsx<-matrix(rnorm(n* p), n, p)# Predictor matrixbeta<-c(rep(1, s.fix),rep(0, p- s.fix))# True fixed effectsu<-cbind(matrix(rnorm(m* s.rand), m, s.rand),matrix(0, m, p- s.rand))# True random effectscluster<-sample(1:m, n,replace =TRUE)# Cluster labelsxb<-rowSums(x*sweep(u,2, beta,'+')[cluster, ])# x %*% (beta + u) matrixy<-rnorm(n, xb)# Response vectorOf the five candidate predictors, the first two have nonzero fixedeffects. Only the first predictor has a nonzero random effect.
The values of\(\lambda\) areautomatically computed from the data, providing a path of solutions fromthe null model (intercept only) to the full model (all predictorsincluded). The fixed effects and random effects from the path of fitscan be extracted using thefixef() andranef()functions.
fixef(fit)#> [,1] [,2] [,3] [,4] [,5]#> [1,] -0.2922322 -0.1167263 -0.11498018 -0.11484951 -0.10708102#> [2,] 0.0000000 1.1104188 1.11240541 1.12407682 1.12770012#> [3,] 0.0000000 1.1609182 1.16109489 1.17013772 1.17421123#> [4,] 0.0000000 0.0000000 0.00000000 0.00000000 -0.04771683#> [5,] 0.0000000 0.0000000 0.00000000 -0.05988147 -0.05850389#> [6,] 0.0000000 0.0000000 0.04929141 0.04396767 0.04550589ranef(fit)#> , , 1#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] 0 0 0 0 0 0#> [2,] 0 0 0 0 0 0#> [3,] 0 0 0 0 0 0#> [4,] 0 0 0 0 0 0#>#> , , 2#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] 0 -1.14669569 0 0 0 0#> [2,] 0 0.03435146 0 0 0 0#> [3,] 0 0.75530136 0 0 0 0#> [4,] 0 0.35703619 0 0 0 0#>#> , , 3#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] 0 -1.15635158 0 0 0 0#> [2,] 0 0.03829802 0 0 0 0#> [3,] 0 0.75971247 0 0 0 0#> [4,] 0 0.35850014 0 0 0 0#>#> , , 4#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] 0 -1.1514971 0 0 0 0#> [2,] 0 0.0327921 0 0 0 0#> [3,] 0 0.7502685 0 0 0 0#> [4,] 0 0.3687197 0 0 0 0#>#> , , 5#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] 0 -1.14421868 0 0 0 0#> [2,] 0 0.03769651 0 0 0 0#> [3,] 0 0.74853993 0 0 0 0#> [4,] 0 0.35795599 0 0 0 0Each column in the output offixef() corresponds to aset of fixed effects for a particular value of\(\lambda\), with the first row containingintercept terms. In the output ofranef(), each slicecorresponds to a set of random effects for a particular value of\(\lambda\), with each row containing therandom effects for a given cluster.
When making predictions, it is often useful to add the fixed andrandom effects to get the cluster-specific coefficients. Thecoef() function provides this functionality.
coef(fit)#> , , 1#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] -0.2922322 0 0 0 0 0#> [2,] -0.2922322 0 0 0 0 0#> [3,] -0.2922322 0 0 0 0 0#> [4,] -0.2922322 0 0 0 0 0#>#> , , 2#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] -0.1167263 -0.03627693 1.160918 0 0 0#> [2,] -0.1167263 1.14477022 1.160918 0 0 0#> [3,] -0.1167263 1.86572012 1.160918 0 0 0#> [4,] -0.1167263 1.46745495 1.160918 0 0 0#>#> , , 3#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] -0.1149802 -0.04394616 1.161095 0 0 0.04929141#> [2,] -0.1149802 1.15070343 1.161095 0 0 0.04929141#> [3,] -0.1149802 1.87211788 1.161095 0 0 0.04929141#> [4,] -0.1149802 1.47090555 1.161095 0 0 0.04929141#>#> , , 4#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] -0.1148495 -0.02742033 1.170138 0 -0.05988147 0.04396767#> [2,] -0.1148495 1.15686892 1.170138 0 -0.05988147 0.04396767#> [3,] -0.1148495 1.87434529 1.170138 0 -0.05988147 0.04396767#> [4,] -0.1148495 1.49279648 1.170138 0 -0.05988147 0.04396767#>#> , , 5#>#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] -0.107081 -0.01651856 1.174211 -0.04771683 -0.05850389 0.04550589#> [2,] -0.107081 1.16539663 1.174211 -0.04771683 -0.05850389 0.04550589#> [3,] -0.107081 1.87624005 1.174211 -0.04771683 -0.05850389 0.04550589#> [4,] -0.107081 1.48565611 1.174211 -0.04771683 -0.05850389 0.04550589Each row in each of these slices represents the fixed effects plusthe random effects for a given cluster, e.g., the second row represents\(\hat{\boldsymbol{\beta}}+\hat{\mathbf{u}}_2\).
Thepredict() function is available for makingpredictions on new data.
x.new<- x[1:3, ]cluster.new<- cluster[1:3]predict(fit, x.new, cluster.new)#> [,1] [,2] [,3] [,4] [,5]#> [1,] -0.2922322 0.40829027 0.3588954 0.3840867 0.35454510#> [2,] -0.2922322 -0.35024290 -0.3451526 -0.2907129 -0.31701760#> [3,] -0.2922322 -0.07945345 -0.1067836 -0.0751470 -0.06503485Again, the columns represent predictions for different values of\(\lambda\).
In practice,\(\lambda\) usuallyneeds to be cross-validated. Thecv.glmmsel() functionprovides a convenient way to perform cross-validation.
glmmsel() does not need to be run after usingcv.glmmsel(), as the latter calls the former and saves theresult asfit$fit.
Thecoef() andpredict() functions appliedto the output ofcv.glmmsel() return the resultcorresponding to the value of\(\lambda\) that minimises thecross-validation loss.
Currently,glmmsel supports Gaussian likelihoods(default) and binomial likelihoods. To use a binomial likelihood andperform a logistic linear mixed model fit, setfamily = 'binomial'.
The primary algorithm drivingglmmsel is coordinatedescent. Sometimes when the predictors are strongly correlated, themodels fit by coordinate descent can be improved using local search.This algorithm runs on top of coordinate descent. To use local search,setlocal.search = TRUE.
x<-0.2*matrix(rnorm(n* p), n, p)+0.8*matrix(rnorm(n), n, p)xb<-rowSums(x*sweep(u,2, beta,'+')[cluster, ])y<-rnorm(n, xb)fit<-cv.glmmsel(x, y, cluster)coef(fit)#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] 0.0772794 0.2402767 1.476905 0 -0.569285 0#> [2,] 0.0772794 0.8580355 1.476905 0 -0.569285 0#> [3,] 0.0772794 1.4901887 1.476905 0 -0.569285 0#> [4,] 0.0772794 1.8003791 1.476905 0 -0.569285 0fit<-cv.glmmsel(x, y, cluster,local.search =TRUE)coef(fit)#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] 0.04734797 -0.09505438 1.264529 0 0 0#> [2,] 0.04734797 0.47932137 1.264529 0 0 0#> [3,] 0.04734797 1.11757935 1.264529 0 0 0#> [4,] 0.04734797 1.48196842 1.264529 0 0 0The correct predictors are not selected without local search in thishigh-correlation example.