# 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)
# random effects
head(ranef(fm1))
# detailed output
summary(fm1)
# fitted values for the 'mean subject', i.e., with
# random effects values equal to 0
head(fitted(fm1, type = "mean_subject"))
# fitted values for the conditioning on the estimated random effects
head(fitted(fm1, type = "subject_specific"))
##############
# \donttest{
fm2 <- mixed_model(fixed = y ~ year, random = ~ 1 | id, data = DF,
family = binomial())
# likelihood ratio test between the two models
anova(fm2, fm1)
# the same hypothesis but with a Wald test
anova(fm1, L = rbind(c(0, 0, 1, 0), c(0, 0, 0, 1)))
##############
# 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)
# }
Run the code above in your browser using DataLab