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 fixedeffects part of the model, including the outcome. 

random  a formula for the randomeffects part of the model. This should only contain
the righthand side part, e.g., 
data  a data.frame containing the variables required in 
family  a 
weights  a numeric vector of weights. These are simple multipliers on the loglikelihood 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 
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 
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,

control  a list with the following components:

...  arguments passed to 
General: The mixed_model()
function fits mixed effects models in which the
integrals over the random effects in the definition of the marginal loglikelihood cannot
be solved analytically and need to be approximated. The function works under the
assumption of normally distributed random effects with mean zero and variancecovariance
matrix \(D\). These integrals are approximated numerically using an adaptive
GaussHermite quadrature rule. Using the control argument nAGQ
, the user can
specify the number of quadrature points used in the approximation.
Userdefined 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
loglikelihood using function optim
and the algorithm specified by
optim_method
. For stability and speed, the derivative of the loglikelihood with
respect to the parameters are internally programmed.
An object of class "MixMod"
with components:
a numeric vector with the estimated fixed effects.
a numeric vector with the estimated extra parameters.
a numeric matrix denoting the estimated covariance matrix of the random effects.
a numeric matrix with the empirical Bayes estimates of the random effects.
a list of numeric matrices with the posterior variances of the random effects.
a numeric scalar denoting the loglikelihood value at the end of the optimization procedure.
a numeric matrix denoting the Hessian matrix at the end of the optimization procedure.
a logical indicating whether convergence was attained.
a copy of the data
argument.
a copy of the grouping variable from data
.
a character string with the name of the grouping variable.
a list with two terms components, termsX
derived from the fixed
argument, and termsZ
derived from the random
argument.
a list with two model.frame components, mfX
derived from the
fixed
argument, and mfZ
derived from the random
argument.
a copy of the (userspecific) control
argument.
a list of functions used in the optimization procedure.
a copy of the family
argument.
the matched call.
Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl
# simulate some data set.seed(123L) n < 500 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(K1, 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.875293208 0.204144801 0.648693468 0.006044763#> (Intercept) #> 1 2.407943 #> 2 2.395551 #> 3 3.878526 #> 4 1.686797 #> 5 2.109013 #> 6 2.951583#> #> Call: #> mixed_model(fixed = y ~ year * group, random = ~1  id, data = DF, #> family = binomial()) #> #> Data Descriptives: #> Number of Observations: 7500 #> Number of Groups: 500 #> #> Model: #> family: binomial #> link: logit #> #> Fit statistics: #> log.Lik AIC BIC #> 3062.201 6134.403 6155.476 #> #> Random effects covariance matrix: #> StdDev #> (Intercept) 2.954122 #> #> Fixed effects: #> Estimate Std.Err zvalue pvalue #> (Intercept) 2.8753 0.2246 12.8013 < 1e04 #> year 0.2041 0.0087 23.3726 < 1e04 #> group1 0.6487 0.3172 2.0448 0.040871 #> year:group1 0.0060 0.0122 0.4970 0.619191 #> #> Integration: #> method: adaptive GaussHermite quadrature rule #> quadrature points: 11 #> #> Optimization: #> method: hybrid EM and quasiNewton #> 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.05338851 0.06643478 0.19661180 0.31259174 0.36316625 0.36704727# fitted values for the conditioning on the estimated random effects head(fitted(fm1, type = "subject_specific"))#> 1 2 3 4 5 6 #> 0.3852436 0.4415565 0.7311246 0.8347831 0.8636915 0.8656506############## # \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 6137.21 6149.86 3065.61 #> fm1 6134.40 6155.48 3062.20 6.81 2 0.0332 #>#> #> Marginal Wald Tests Table #> #> Userdefined contrasts matrix: #> (Intr) year group1 yr:gr1 #> 0 0 1 0 #> 0 0 0 1 #> #> Chisq df Pr(>Chi) #> 6.9211 2 0.0314 #>############## # 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 = "Followup 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 = "Followup 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 = "Followup 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 5200.27 5225.56 2594.13 #> fm1_slp 5201.62 5231.12 2593.81 0.65 1 0.419 #># }