marginal-conditional.R

library(lme4)
library(ggplot2)
n.item <- 1000
n.trial <- 20
sd.item <- 1.5
intercept <- -.25
trend <- 2

dataset <- expand.grid(
  covariate = seq(-2, 2, length = 21),
  item = seq_len(n.item)
)
re <- rnorm(n.item, mean = 0, sd = sd.item)
dataset$eta <- with(dataset, intercept + trend * covariate + re[item])
dataset$mu <- plogis(dataset$eta)
dataset$success <- rbinom(nrow(dataset), size = n.trial, prob = dataset$mu)
dataset$failure <- n.trial - dataset$success

ggplot(dataset, aes(x = covariate, y = mu)) + geom_line(aes(group = item), alpha = 0.1) + stat_summary(colour = "red", geom = "line", fun.y = mean)

model.marginal <- glm(cbind(success, failure) ~ covariate, data = dataset, family = binomial)
model.conditional <- glmer(cbind(success, failure) ~ covariate + (1|item), data = dataset, family = binomial)

coef(model.marginal)
fixef(model.conditional)

newdata <- unique(dataset[, "covariate", drop = FALSE])
newdata$marginal <- predict(model.marginal, newdata = newdata, type = "response")
newdata$conditional <- predict(model.conditional, newdata = newdata, type = "response", re.form = ~ 0)

ggplot(dataset, aes(x = covariate, y = mu)) + geom_line(aes(group = item), alpha = 0.1) + stat_summary(colour = "red", geom = "line", fun.y = mean) + geom_line(data = newdata, aes(y = marginal), colour = "blue") + geom_line(data = newdata, aes(y = conditional), colour = "green")