`mixed_model.Rd`

Fits generalized linear mixed effects models under maximum likelihood using adaptive Gaussian quadrature.

```
mixed_model(fixed, random, data, family, weights = NULL,
na.action = na.exclude, zi_fixed = NULL, zi_random = NULL,
penalized = FALSE, n_phis = NULL, initial_values = NULL,
control = list(), ...)
```

- fixed
a formula for the fixed-effects part of the model, including the outcome.

- random
a formula for the random-effects part of the model. This should only contain the right-hand side part, e.g.,

`~ time | id`

, where`time`

is a variable, and`id`

the grouping factor. When the symbol`||`

is used in the definition of this argument (instead of`|`

), then the covariance matrix of the random effects is assumed to be diagonal.- data
a data.frame containing the variables required in

`fixed`

and`random`

.- family
a

`family`

object specifying the type of the repeatedly measured response variable, e.g.,`binomial()`

or`poisson()`

. The function also allows for user-defined family objects, but with specific extra components; see the example is`negative.binomial`

for more details. Contrary to the standard practice in model fitting R functions with a`family`

argument (e.g.,`glm`

) in which the default family is`gaussian()`

, in`mixed_model()`

no default is provided. If the users wish to fit a mixed model for a Gaussian outcome, this could be done with function`lme()`

from the**nlme**package or function`lmer()`

from the**lme4**package.- weights
a numeric vector of weights. These are simple multipliers on the log-likelihood contributions of each group/cluster, i.e., we presume that there are multiple replicates of each group/cluster denoted by the weights. The length of `weights` need to be equal to the number of independent groups/clusters in the data.

- na.action
what to do with missing values in

`data`

.- zi_fixed, zi_random
formulas for the fixed and random effects of the zero inflated part.

- penalized
logical or a list. If logical and equal to

`FALSE`

, then no penalty is used. If logical and equal to`TRUE`

, for the fixed effects a Student's-t penalty/prior with mean 0, scale equal to 1 and 3 degrees of freedom is used. If a list, then it is expected to have the components`pen_mu`

,`pen_sigma`

and`pen_df`

, denoting the mean, scale and degrees of freedom of the Student's-t penalty/prior for the fixed effects.- n_phis
a numeric scalar; in case the family corresponds to a distribution that has extra (dispersion/shape) parameters, you need to specify how many extra parameters are needed.

- initial_values
a list of initial values. This can have up to three components, namely,

- betas
a numeric vector of fixed effects. This can also be

`family`

object. In this case initial values for the fixed effects will be calculated by using`glm`

to the data ignoring the correlations in the repeated measurements. For example, for a negative binomial response outcome, we could set`betas = poisson()`

.- D
a numeric matrix denoting the covariance matrix of the random effects.

- phis
a numeric vector for the extra (dispersion/shape) parameters.

- control
a list with the following components:

- iter_EM
numeric scalar denoting the number of EM iterations; default is 30.

- iter_qN_outer
numeric scalar denoting the number of outer iterations during the quasi-Newton phase; default is 15. In each outer iteration the locations of the quadrature points are updated.

- iter_qN
numeric scalar denoting the starting number of iterations for the quasi-Newton; default is 10.

- iter_qN_incr
numeric scalar denoting the increment in

`iter_qN`

for each outer iteration; default is 10.- optimizer
character string denoting the optimizer to be used; available options are

`"optim"`

(default),`"nlminb"`

and`"optimParallel"`

, the last option implemented in the**optimParallel**package.- optim_method
character string denoting the type of

`optim`

algorithm to be used when`optimizer = "optim"`

; default is the BFGS algorithm.- parscale_betas
the control argument

`parscale`

of`optim`

for the fixed-effects; default is 0.1.- parscale_D
the control argument

`parscale`

of`optim`

for the unique element of the covariance matrix of the random effects; default is 0.01.- parscale_phis
the control argument

`parscale`

of`optim`

for the extra (dispersion/shape) parameters; default is 0.01.- tol1, tol2, tol3
numeric scalars controlling tolerances for declaring convergence;

`tol1`

and`tol2`

are for checking convergence in successive parameter values;`tol3`

is similar to`reltop`

of`optim`

; default values are`1e-03`

,`1e-04`

, and`1e-08`

, respectively.- numeric_deriv
character string denoting the type of numerical derivatives to be used. Options are

`"fd"`

for forward differences, and`cd`

for central difference; default is`"fd"`

.- nAGQ
numeric scalar denoting the number of quadrature points; default is 11 when the number of random effects is one or two, and 7 otherwise.

- update_GH_every
numeric scalar denoting every how many iterations to update the quadrature points during the EM-phase; default is 10.

- max_coef_value
numeric scalar denoting the maximum allowable value for the fixed effects coefficients during the optimization; default is 10.

- max_phis_value
numeric scalar denoting the maximum allowable value for the shape/dispersion parameter of the negative binomial distribution during the optimization; default is

`exp(10)`

.- verbose
logical; print information during the optimization phase; default is

`FALSE`

.

- ...
arguments passed to

`control`

.

**General:** The `mixed_model()`

function fits mixed effects models in which the
integrals over the random effects in the definition of the marginal log-likelihood cannot
be solved analytically and need to be approximated. The function works under the
assumption of normally distributed random effects with mean zero and variance-covariance
matrix \(D\). These integrals are approximated numerically using an adaptive
Gauss-Hermite quadrature rule. Using the control argument `nAGQ`

, the user can
specify the number of quadrature points used in the approximation.

**User-defined family:** The user can define its own family object; for an example,
see the help page of `negative.binomial`

.

**Optimization:** A hybrid approach is used, starting with `iter_EM`

iterations
and unless convergence was achieved it continuous with a direct optimization of the
log-likelihood using function `optim`

and the algorithm specified by
`optim_method`

. For stability and speed, the derivative of the log-likelihood with
respect to the parameters are internally programmed.

An object of class `"MixMod"`

with components:

- coefficients
a numeric vector with the estimated fixed effects.

- phis
a numeric vector with the estimated extra parameters.

- D
a numeric matrix denoting the estimated covariance matrix of the random effects.

- post_modes
a numeric matrix with the empirical Bayes estimates of the random effects.

- post_vars
a list of numeric matrices with the posterior variances of the random effects.

- logLik
a numeric scalar denoting the log-likelihood value at the end of the optimization procedure.

- Hessian
a numeric matrix denoting the Hessian matrix at the end of the optimization procedure.

- converged
a logical indicating whether convergence was attained.

- data
a copy of the

`data`

argument.- id
a copy of the grouping variable from

`data`

.- id_name
a character string with the name of the grouping variable.

- Terms
a list with two terms components,

`termsX`

derived from the`fixed`

argument, and`termsZ`

derived from the`random`

argument.- model_frames
a list with two model.frame components,

`mfX`

derived from the`fixed`

argument, and`mfZ`

derived from the`random`

argument.- control
a copy of the (user-specific)

`control`

argument.- Funs
a list of functions used in the optimization procedure.

- family
a copy of the

`family`

argument.- call
the matched call.

```
# simulate some data
set.seed(123L)
n <- 200
K <- 15
t.max <- 25
betas <- c(-2.13, -0.25, 0.24, -0.05)
D <- matrix(0, 2, 2)
D[1:2, 1:2] <- c(0.48, -0.08, -0.08, 0.18)
times <- c(replicate(n, c(0, sort(runif(K-1, 0, t.max)))))
group <- sample(rep(0:1, each = n/2))
DF <- data.frame(year = times, group = factor(rep(group, each = K)))
X <- model.matrix(~ group * year, data = DF)
Z <- model.matrix(~ year, data = DF)
b <- cbind(rnorm(n, sd = sqrt(D[1, 1])), rnorm(n, sd = sqrt(D[2, 2])))
id <- rep(1:n, each = K)
eta.y <- as.vector(X %*% betas + rowSums(Z * b[id, ]))
DF$y <- rbinom(n * K, 1, plogis(eta.y))
DF$id <- factor(id)
################################################
fm1 <- mixed_model(fixed = y ~ year * group, random = ~ 1 | id, data = DF,
family = binomial())
# fixed effects
fixef(fm1)
#> (Intercept) year group1 year:group1
#> -2.51135823 0.23304276 -0.88784195 -0.05174851
# random effects
head(ranef(fm1))
#> (Intercept)
#> 1 0.2428149
#> 2 2.1943191
#> 3 1.2898486
#> 4 -1.2667384
#> 5 -0.7294567
#> 6 3.1614903
# detailed output
summary(fm1)
#>
#> Call:
#> mixed_model(fixed = y ~ year * group, random = ~1 | id, data = DF,
#> family = binomial())
#>
#> Data Descriptives:
#> Number of Observations: 3000
#> Number of Groups: 200
#>
#> Model:
#> family: binomial
#> link: logit
#>
#> Fit statistics:
#> log.Lik AIC BIC
#> -1225.973 2461.947 2478.438
#>
#> Random effects covariance matrix:
#> StdDev
#> (Intercept) 2.851543
#>
#> Fixed effects:
#> Estimate Std.Err z-value p-value
#> (Intercept) -2.5114 0.3430 -7.3216 < 1e-04
#> year 0.2330 0.0155 15.0812 < 1e-04
#> group1 -0.8878 0.4816 -1.8437 0.0652272
#> year:group1 -0.0517 0.0195 -2.6534 0.0079694
#>
#> Integration:
#> method: adaptive Gauss-Hermite quadrature rule
#> quadrature points: 11
#>
#> Optimization:
#> method: hybrid EM and quasi-Newton
#> converged: TRUE
# fitted values for the 'mean subject', i.e., with
# random effects values equal to 0
head(fitted(fm1, type = "mean_subject"))
#> 1 2 3 4 5 6
#> 0.03232047 0.03944053 0.10950605 0.17572685 0.20676636 0.20921567
# fitted values for the conditioning on the estimated random effects
head(fitted(fm1, type = "subject_specific"))
#> 1 2 3 4 5 6
#> 0.04084042 0.04974091 0.13552321 0.21370151 0.24941898 0.25221291
##############
# \donttest{
fm2 <- mixed_model(fixed = y ~ year, random = ~ 1 | id, data = DF,
family = binomial())
# likelihood ratio test between the two models
anova(fm2, fm1)
#>
#> AIC BIC log.Lik LRT df p.value
#> fm2 2477.04 2486.94 -1235.52
#> fm1 2461.95 2478.44 -1225.97 19.1 2 1e-04
#>
# the same hypothesis but with a Wald test
anova(fm1, L = rbind(c(0, 0, 1, 0), c(0, 0, 0, 1)))
#>
#> Marginal Wald Tests Table
#>
#> User-defined contrasts matrix:
#> (Intr) year group1 yr:gr1
#> 0 0 1 0
#> 0 0 0 1
#>
#> Chisq df Pr(>|Chi|)
#> 18.7743 2 1e-04
#>
##############
# An effects plot for the mean subject (i.e., with random effects equal to 0)
nDF <- with(DF, expand.grid(year = seq(min(year), max(year), length.out = 15),
group = levels(group)))
plot_data <- effectPlotData(fm2, nDF)
require("lattice")
xyplot(pred + low + upp ~ year | group, data = plot_data,
type = "l", lty = c(1, 2, 2), col = c(2, 1, 1), lwd = 2,
xlab = "Follow-up time", ylab = "log odds")
expit <- function (x) exp(x) / (1 + exp(x))
xyplot(expit(pred) + expit(low) + expit(upp) ~ year | group, data = plot_data,
type = "l", lty = c(1, 2, 2), col = c(2, 1, 1), lwd = 2,
xlab = "Follow-up time", ylab = "Probabilities")
# An effects plots for the marginal probabilities
plot_data_m <- effectPlotData(fm2, nDF, marginal = TRUE, cores = 1L)
expit <- function (x) exp(x) / (1 + exp(x))
xyplot(expit(pred) + expit(low) + expit(upp) ~ year | group, data = plot_data_m,
type = "l", lty = c(1, 2, 2), col = c(2, 1, 1), lwd = 2,
xlab = "Follow-up time", ylab = "Probabilities")
##############
# include random slopes
fm1_slp <- update(fm1, random = ~ year | id)
# increase the number of quadrature points to 15
fm1_slp_q15 <- update(fm1_slp, nAGQ = 15)
# a diagonal covariance matrix for the random effects
fm1_slp_diag <- update(fm1, random = ~ year || id)
anova(fm1_slp_diag, fm1_slp)
#>
#> AIC BIC log.Lik LRT df p.value
#> fm1_slp_diag 2112.08 2131.87 -1050.04
#> fm1_slp 2113.83 2136.91 -1049.91 0.25 1 0.6141
#>
# }
```