Binomial Generalized Linear Mixed Models, or binomial GLMMs, are useful for modeling binary outcomes for repeated or clustered measures. For example, let’s say we design a study that tracks what college students eat over the course of 2 weeks, and we’re interested in whether or not they eat vegetables each day. For each student we’ll have 14 binary events: eat vegetables or not. Using a binomial GLMM we could model the probability of eating vegetables daily given various predictors such as sex of the student, race of the student, and/or some “treatment” we applied to a subset of the students, such as a nutrition class. Since each student is observed over the course of multiple days, we have repeated measures and thus the need for a mixed-effect model.

Before we show how to implement and interpret a binomial GLMM, we’ll first simulate some data that is appropriate for a binomial GLMM. If we know how to simulate data for a given model, then we have a better understanding of the model’s assumptions and coefficients. We’ll use R to simulate the data and perform the modeling.

To begin we simulate data for four predictor variables:

- id (indentifying number of subject)
- day (the day of observation)
- trt (treatment status: control or treat)
- sex (male or female)

We set “n” to 250 which means we’ll have 250 subjects in our data. Then we use `seq(n)`

to generate a vector of integers ranging 1 to 250 and assign to “id”. These will be our subject ids. We assign to “day” the integers 1 – 14 and then use `expand.grid`

to create all combinations of “id” and “day” and assign to “d”, which is a data frame. To simulate the “trt” and “sex” variables we simply sample 250 times from the possible values with replacement. The `set.seed(1)`

function ensures we always simulate the same values, which you’ll need to run if you want to replicate the results of this article. After we simulate the “trt” and “sex” values, we need to repeat them 14 times each for each subject id. We do that by using subsetting brackets and assigning the result to the data frame “d”. Finally we sort the data by “id” and “day”, reset the row numbers, and look at the first 15 records.

n <- 250 id <- seq(n) day <- 1:14 d <- expand.grid(id = id, day = day) set.seed(1) trt <- sample(c("control", "treat"), size = n, replace = TRUE) sex <- sample(c("female", "male"), size = n, replace = TRUE) d$trt <- trt[d$id] d$sex <- sex[d$id] d <- d[order(d$id, d$day),] rownames(d) <- NULL head(d, n = 15) ## id day trt sex ## 1 1 1 control female ## 2 1 2 control female ## 3 1 3 control female ## 4 1 4 control female ## 5 1 5 control female ## 6 1 6 control female ## 7 1 7 control female ## 8 1 8 control female ## 9 1 9 control female ## 10 1 10 control female ## 11 1 11 control female ## 12 1 12 control female ## 13 1 13 control female ## 14 1 14 control female ## 15 2 1 treat male

We see that subject 1 is a female in the control group, with 14 observations over 14 days.

Now we’re ready to use this data to simulate zeroes and ones for the event of eating vegetables or not, where a 0 means “did not eat vegetables” and a 1 means “ate vegetables”. We want to simulate the zeroes and ones using probabilities. We’ll use the `rbinom`

function to do this, which generates zeroes or ones from a binomial distribution. For example to generate 5 flips of a single coin with a probability of success (say, heads) of 0.5, we can do the following:

rbinom(n = 5, size = 1, prob = 0.5) ## [1] 1 1 1 1 0

But we don’t have a constant probability. We’ll have a probability that changes based on the sex of the subject and whether they were in the treatment group or not. Since we’re simulating the data, we get to pick the probabilities.

- female, control: p = 0.40
- female, treated: p = 0.85
- male, control: p = 0.30
- male, treated: p = 0.50

These are what we might call *fixed effects*. We’re saying the effects of treatment and sex are fixed, no matter who you are, where you’re from, how old you are, etc. This may not be a plausible assumption in real life, but that’s what we’re assuming when we simulate this particular data set. Also notice these effects *interact*. The effect of treatment increases the female probability by 0.45, but only increases the male probability by 0.20. The effect of treatment depends on sex, which implies they interact.

But recall we’re observing the same person 14 days in a row. This means we won’t have independent observations. Those 14 observations on the same person may be correlated in some way. The simplest way to induce this correlation is to add a *random effect* for each person. For example, maybe a male student grew up in a family that had a garden in the backyard and was raised eating homegrown vegetables. His random effect might be an additional 0.10 probability. So if he was in the control group, his probability might be 0.30 (fixed) + 0.10 (random) = 0.40. So now we have a *mix* of fixed effects and random effects.

Let’s add our fixed effect probabilities to the data frame. We use the `interaction`

function to create a 4-level factor called “trtsex” that indicates whether a subject is male/female and treated/control. Then we create a vector of probabilities and name the vector using the level names of the “trtsex” variable. This becomes a lookup table which we can use with “trtsex” to quickly determine a subject’s fixed effect probability.

d$trtsex <- interaction(d$trt, d$sex) probs <- c(0.40, 0.85, 0.30, 0.50) names(probs) <- levels(d$trtsex) d$p <- probs[d$trtsex]

Now let’s add random effect probabilities. We’ll do this by drawing `n`

random samples from a Normal distribution with a mean 0 and a standard deviation of 0.03. This will create random amounts of probability that we add or subtract from a subject’s fixed effect probability. Again we use “r_probs” as a lookup table and use subsetting brackets to pick the random probability according to subject id. Notice it gets repeated 14 times for each subject.

set.seed(3) r_probs <- rnorm(n = n, mean = 0, sd = 0.03) d$random_p <- r_probs[d$id] d$p <- d$p + d$random_p

And now we can use the probabilities to generate zeroes and ones with the `binom`

function.

d$y <- rbinom(n = nrow(d), size = 1, prob = d$p)

With that our data is simulated. Let’s inspect the first few records.

head(d[c("id", "day", "trt", "sex", "p", "y")]) ## id day trt sex p y ## 1 1 1 control female 0.371142 0 ## 2 1 2 control female 0.371142 0 ## 3 1 3 control female 0.371142 1 ## 4 1 4 control female 0.371142 0 ## 5 1 5 control female 0.371142 0 ## 6 1 6 control female 0.371142 0

Subject 1 (female/control) reported eating vegetables 1 time in the first 6 days. And we see her probability is 0.37. Her random effect is about 0.37 – 0.40 = -0.03.

Let’s look at subject 5.

head(d[d$id == 5, c("id", "day", "trt", "sex", "p", "y")]) ## id day trt sex p y ## 57 5 1 treat female 0.8558735 1 ## 58 5 2 treat female 0.8558735 1 ## 59 5 3 treat female 0.8558735 1 ## 60 5 4 treat female 0.8558735 1 ## 61 5 5 treat female 0.8558735 1 ## 62 5 6 treat female 0.8558735 1

Subject 5 (female/treat) reported eating vegetables everyday in the first 6 days. And we see her probability is about 0.86. Her random effect is around 0.86 – 0.85 = 0.01.

**Full Disclosure**: The way we just simulated our data is not technically correct! We could have easily generated a random effect that pushed a subject’s probability beyond 1 or below 0! However we thought this approach was a good place to start because it allows us to think in terms of probabilities. The correct way is demonstrated later in the article and, as we’ll see, involves log-odds which are much less intuitive.

Now let’s work backwards and pretend we don’t know the probabilities we defined above. We just have “id”, “day”, “trt”, “sex” and “y”. We want to fit a binomial GLMM model to see how “sex” and “trt” affect the probability of eating vegetables. In real life we won’t know if or how they affect the probability. We may have a hunch, but all we can really do is propose a model and see how well it fits the data. To do this we’ll use the `glmer`

function in the lme4 package.

Let’s make life easy for us and propose the “correct” model, which we know since we simulated the data. We model y as a function sex, trt, and their interaction: `y ~ trt * sex`

. We also include a random effect for each subject with `+ (1|id)`

, which is also known as a *random intercept* since the model estimates a different intercept for each subject. We can read `1|id`

as “the intercept is conditional on subject id.” (In R model syntax, `1`

represents the intercept.) We also specify `family = binomial`

to indicate we assume the response was drawn from a binomial distribution, which is correct in this case.

library(lme4) m <- glmer(y ~ trt * sex + (1|id), data = d, family = binomial) summary(m, corr = FALSE) ## Generalized linear mixed model fit by maximum likelihood (Laplace ## Approximation) [glmerMod] ## Family: binomial ( logit ) ## Formula: y ~ trt * sex + (1 | id) ## Data: d ## ## AIC BIC logLik deviance df.resid ## 4132.4 4163.3 -2061.2 4122.4 3495 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -2.6008 -0.8434 0.3774 0.9704 1.6878 ## ## Random effects: ## Groups Name Variance Std.Dev. ## id (Intercept) 0.03944 0.1986 ## Number of obs: 3500, groups: id, 250 ## ## Fixed effects: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -0.30840 0.07127 -4.327 1.51e-05 *** ## trttreat 2.18835 0.12423 17.615 < 2e-16 *** ## sexmale -0.63440 0.10859 -5.842 5.16e-09 *** ## trttreat:sexmale -1.21669 0.16556 -7.349 2.00e-13 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In the summary we see Random effects and Fixed effects. The fixed effect coefficients are not on the probability scale but on the log-odds, or *logit*, scale. The Logit transformation takes values ranging from 0 to 1 (probabilities) and transforms them to values ranging from `-Inf`

to `+Inf`

. This allows us to create additive linear models without worrying about going above 1 or below 0. To get probabilities out of our model, we need to use the *inverse logit*. There is function for this in base R called `plogis`

. For example, the predicted log-odds a female in the control group eats vegetables is the intercept: -0.30840. The inverse logit of that is:

plogis(-0.30840) ## [1] 0.4235053

The predicted probability of 0.42 is pretty close to the 0.40 value we used to generate the data.

The predicted log-odds a female in the treatment group eats vegetables is the intercept plus the coefficient for `trttreat`

: -0.30840 + 2.18835 = 1.87995. The inverse logit is:

plogis(1.87995) ## [1] 0.8676054

The predicted probability of 0.87 is close to the 0.85 value we used to generate the data.

The predicted log-odds a male in the control group eats vegetables is the intercept plus the coefficient for `sexmale`

: -0.30840 + -0.63440 = -0.9428. The inverse logit is:

plogis(-0.9428) ## [1] 0.2803351

The predicted probability of 0.28 is close to the 0.30 value we used to generate the data.

Finally the predicted log-odds a male in the treatment group eats vegetables is the sum of all the coefficients: -0.30840 + 2.18835 + -0.63440 + -1.21669 = 0.02886. The inverse logit is:

plogis(0.02886) ## [1] 0.5072145

The predicted probability of 0.51 is very close to the 0.50 value we used to generate the data.

Instead of using `plogis`

, we can simply use the `predict`

function. The `re.form = NA`

argument says to only estimate the fixed effects. The data frame given to the `newdata`

argument represents all possible combinations of subject type.

predict(m, type = "response", re.form = NA, newdata = data.frame(sex = c("female", "female", "male", "male"), trt = c("control", "treat", "control", "treat")) ) ## 1 2 3 4 ## 0.4235043 0.8676050 0.2803350 0.5072137

The predicted probabilities are the same as we calculated “by hand” above (minus some rounding error).

We could also interpret the model coefficients. To do this we exponentiate them to get an *odds ratio*. For example, exponentiating the treat coefficient yields 8.92.

exp(2.18835) ## [1] 8.920482

That says the odds of “success” for a female in the treatment group is about 8.9 times higher than the odds of success for a female in the control group.

We could also exponentiate the confidence interval to get a lower and upper bound on the odds ratio:

exp(confint(m, parm = "trttreat")) ## Computing profile confidence intervals ... ## 2.5 % 97.5 % ## trttreat 7.023558 11.44831

So we might conclude the odds of “success” are at least 7 times higher for the female treated versus the female control.

But odds ratios are relative and can sometimes sound alarmist. For example, let’s say the effect of “treat” caused the probability of success to increase from 0.001 to 0.009. Neither probability is very likely, but the odds ratio is about 9.

odds1 <- 0.001/(1 - 0.001) odds2 <- 0.009/(1 - 0.009) odds2/odds1 ## [1] 9.072654

In this case describing the treatment effect as making the odds of success 9 times more likely may suggest to an unsuspecting reader that it’s more efficacious than it really is. That’s why it’s a good idea to examine expected probabilities in addition to odds ratios.

In the Random Effects section we see the estimated Standard Deviation of the Intercept random effect, which we can extract from the model object with the `VarCorr`

function.

VarCorr(m) ## Groups Name Std.Dev. ## id (Intercept) 0.1986

0.20 is the estimate of 0.03, the standard deviation we used to simulate our random probability effects. This seems to be a rather poor estimate. However if we calculate a 95% confidence interval on this estimate, we see that 0.03 falls well within the interval. To do this we call the `confint`

function and specify `parm = "theta_"`

to indicate we want a confidence interval for the estimated standard deviation of the intercept random effect. Setting `oldNames = FALSE`

provides more informative output.

confint(m, parm = "theta_", oldNames = FALSE) ## Computing profile confidence intervals ... ## 2.5 % 97.5 % ## sd_(Intercept)|id 0 0.3453926

This confidence interval extends to 0, suggesting there may be no random variability between subjects. However we *know* we have variability because we generated the data. Plus the structure of the data dictates we should accommodate variance between subjects as well as within subjects. In other words, we don’t have independent observations. So in real life we wouldn’t seriously entertain dropping the random effect.

If we were to report this model using mathematical notation, we might write it in this form:

\[\text{Prob}(y_{j} = 1) = \text{logit}^{-1}(-0.31 + 2.19 \times \text{Treat} – 0.63 \times \text{Male} – 1.22 \times \text{Treat} \times \text{Male} + u_j)\]

where \(j\) represents the subject id and

\[u_j \sim \text{N}(0, 0.20)\]

The \(u_j\) is the random effect for each person. This says each subject’s random effect is assumed to be drawn from a Normal distribution with mean 0 and standard deviation 0.20.

\(\text{logit}^{-1}\) is the *inverse logit* function and it corresponds to the `plogis`

function we used to transform log-odds into probability. In the language of generalized linear models, this is the *link* function. The logit is the link between an additive linear model and a probability outcome.

**Side note:** Previously we mentioned that we had generated our fake data incorrectly because we ran the risk of simulating probabilities outside the range of [0,1]. The correct way would have been to use the logit transformation as just described. However that means working with log-odds instead of probabilities, which are less intuitive. It’s not impossible to do, just inconvenient. Here’s how we could have simulated our response probabilities using log-odds and a logit transformation. See this page for how we derived the log-odds.

d$p2 <- plogis(-0.4054651 + 2.140066*(d$trt == "treat") + -0.4418328*(d$sex == "male") + -1.292768*(d$trt == "treat")*(d$sex == "male") + rnorm(n = n, mean = 0, sd = 0.03)[d$id]) d$y2 <- rbinom(n = nrow(d), size = 1, prob = d$p2)

Returning to our estimated model, `m`

, we can use it to simulate response data with the `simulate`

function. If our model is “good”, it should simulate response data similar to what we observed. Here we use it to simulate 100 different data sets. It essentially takes our 3500 observed predictors, feeds them into the model, and generates a new series of ones and zeroes to indicate whether someone ate a vegetable or not. The result is a 100 column data frame which we save to an object called “sim”.

sim <- simulate(m, nsim = 100)

So what can we do with this data frame? One application is to use it to visualize how our estimated model performs compared to our observed data. We can look at our observed data by simply taking the mean of `y`

by `trt`

and `sex`

using the `aggregate`

function. (Recall the means of zeroes and ones is the proportion of ones.)

obs <- aggregate(d$y, by = list(trt = d$trt, sex = d$sex), mean) obs ## trt sex x ## 1 control female 0.4242424 ## 2 treat female 0.8659341 ## 3 control male 0.2820823 ## 4 treat male 0.5071429

Notice the observed probabilities are very close to our model predicted probabilities.

Now let’s do the same thing to all the columns in the “sim” data frame. We apply the `aggregate`

function to each column, passing to `aggregate`

the same arguments we used for the observed data. Then we combine the results into one data frame using the `do.call`

and `rbind`

functions.

s_p <- lapply(sim, aggregate, by = list(trt = d$trt,sex = d$sex), mean) s_df <- do.call(rbind, s_p)

Now we can plot the observed and simulated data to see how close the model-simulated data comes to the observed data. We use ggplot2 to create this plot. First we use `geom_point`

to plot the “obs” data frame, making the observed proportions appear with a bigger blue point. Then we plot the simulated data in the “s_df” data frame using `geom_jitter`

, which “jitters” the points sideways. The `alpha = 1/2`

argument makes the points transparent and the `shape = 1`

argument turns the points into open circles.

library(ggplot2) ggplot() + geom_point(aes(y = x, x = trt), data = obs, size = 4, color = "blue") + geom_jitter(aes(y = x, x = trt), data = s_df, width = 0.2, height = 0, shape = 1, alpha = 1/2) + facet_wrap(~sex) + labs(x = "Treatment", y = "Predicted Probability", title = "Model with interaction")

We can see our model-simulated data hovers very closely to the observed data, which is not surprising since we fit the “correct” model to the data.

Let’s fit a wrong model and recreate the plot. Below we fit a model without an interaction of “trt” and “sex”. This is not how we simulated the data so we know the model is wrong. After we fit the model we re-run the same code above to simulate responses from our model 100 times and plot with the observed data.

m2 <- glmer(y ~ trt + sex + (1|id), data = d, family = binomial) sim2 <- simulate(m2, nsim = 100) s_p2 <- lapply(sim2, aggregate, by = list(trt = d$trt, sex = d$sex), mean) s_df2 <- do.call(rbind, s_p2) ggplot() + geom_point(aes(y = x, x = trt), data = obs, size = 4, color = "blue") + geom_jitter(aes(y = x, x = trt), data = s_df2, width = 0.2, height = 0, shape = 1, alpha = 1/2) + facet_wrap(~sex) + labs(x = "Treatment", y = "Predicted Probability", title = "Wrong model: without interaction")

Our model is systematically under-predicting and over-predicting probabilities of success. If our model was “good” we would expect it to generate data similar to what we observed.

Logistic regression and mixed-effect modeling are massive topics and we have just touched on the basics. But hopefully you now have a better idea of how the two can be combined to allow us to model the probability of binary events when we have clustered or repeated measures.

For questions or clarifications regarding this article, contact the UVA Library StatLab: statlab@virginia.edu

View the entire collection of UVA Library StatLab articles.

Clay FordStatistical Research Consultant

University of Virginia Library

March 19, 2021