# lme4 mixed model analysis based on sum of count data

I have a dataset with pupils (1000+) that are nested within school classes (100+). Pupils were either in group1 or group2 (different treatments) and had two measurement time points (time0, time2). The outcome variable is a summary score based on the sum of 5 items that were answered (yes=1, no=0).

ID  group   school_class    time    riskOverall_val
1   1       a               0       3
1   1       a               2       2
2   2       a               0       2
2   2       a               2       1
3   1       b               0       3
3   1       b               2       3
4   2       c               0       2
4   2       c               2       4
5   1       c               0       3
5   1       c               2       2


Questions:

a) I have coded the outcome variable (riskOverall_val) as numeric since I was previously getting an error due to selecting family=poisson. Is that appropriate given that the sum is based on count data?

b) Does poisson as family make sense or would an lmer model without any specification be more appropriate since the outcome variable, even though based on count data from previous items, can take values from 0 to 5?

c) Is the group x time interaction appropriate even though I have only two time points?

My current model:

M = glmer(riskOverall_val ~ group + time + group:time + (1|school_class), data=c, family=poisson)


Normally I would also add a random intercept for the pupils (1 | ID) but that produces the following error: boundary (singular) fit: see help('isSingular'). Which I assume means the model is too complex.

Does this strategy make sense or would you handle the data differently or modify the chosen model?

I'd appreciate any help and suggestions!

• Is the time0 measurement taken before the students are "treated" with their assigned treatment? Commented Sep 28, 2022 at 18:55
• Hi, yes, the time0 measurement is taken before the treatment. Commented Sep 29, 2022 at 9:25

This is a before-after (pre-post) experimental design: the outcome is recorded at baseline (at time 0), then the assigned treatment is applied, then the outcome is measured again (at time 2).

Note: I'll use your notation, but it's customary to label time points consecutively: time 0, time 1, time 2, and so on.

So to start with, you can improve the model by including the baseline risk as a covariate in the model. As a rule of thumb, characteristics measured before treatment are not affected by the treatment and are best modeled as covariates.

You can implement this by pivoting the data, so that there is one row per student.

data_wide
#>      ID group school_class time0 time2
#> 1     1     1 a                3     2
#> 2     2     2 a                2     1
#> 3     3     1 b                3     3
#> 4     4     2 c                2     4
#> 5     5     1 c                3     2


The right hand-side of the model formula becomes:

~ time0 + group + (1 | school_class)


There is no need to include an interaction between baseline and group if students were assigned to a treatment group randomly.

The second improvement is to model the outcome (number of "yes" answers out of 5 items) as a binomial random variable, ie., as number of successes in 5 trials. The full model formula is:

items <- 5L
cbind(time2, items - time2) ~ time0 + group + (1 | school_class)


You can fit this model with lme4::glmer to the data in wide (one row per student) format.

model <-
glmer(
cbind(time2, items - time2) ~ time0 + group + (1 | school_class),
family = binomial,
data = data_wide
)


Here is R code to simulate a binomial response and fit a binomial regression.

library("lme4")
library("tidyverse")

data_long <-
here::here(
"data", "590413.csv"
) %>%
data_long

data_wide <- data_long %>%
pivot_wider(
names_from = "time",
names_glue = "time{time}",
values_from = "riskOverall_val"
)
# One row per student
data_wide

n <- 100
N <- 1000
items <- 5L

set.seed(1234)

# The only purpose of this simulation is to illustrate how to fit a binomial regression.
data_wide <- tibble(
ID = seq(N),
school_class = sample(seq(n), N, replace = TRUE),
group = sample(c("1", "2"), N, replace = TRUE),
a0 = rnorm(N, sd = 0.1),
a2 = if_else(group == "1", 0.0, 0.3) + rnorm(N, sd = 0.1),
time0 = rbinom(N, items, 0.5),
time2 = rbinom(N, items, plogis(a2))
)

model <-
glmer(
cbind(time2, items - time2) ~ time0 + group + (1 | school_class),
data = data_wide,
family = binomial
)
summary(model)

• Hi, thank you very much for your response! Could you please help me understand the second part of your answer? I have to create two new variable displaying the number of successful “trials” from 0 to 5 based on the individual dichotomous item scores. This will be called OverallRisk_T0 and OverallRisk_T2. The outcome are two factor variables with 5 levels. m<-glmer(riskOverall_T2 ~ riskOverall_T0 + group+ (1|school_class), family=binomial,data=data_wide) Is this model correct? Commented Oct 3, 2022 at 12:13
• The output shows me the intercept (which I ignore), then it shows different slopes: riskOverall_T01, riskOverall_T02 etc until 05, which I also ignore. Finally, I look at the values of the group slope. This group coefficient tells me whether my treatment produced significant results. Do you agree? Commented Oct 3, 2022 at 12:14
• m<-glmer(riskOverall_T2 ~ riskOverall_T0 + group+ (1|school_class), family=binomial,data=data_wide) is not correct because the outcome riskOverall_T2 is not 0s and 1s. It's integer counts between 0 and 5. (Note: Integer counts, not factors). That's where the notation cbind(time2, items - time2) comes in. It says (number of successes, number of failures). Commented Oct 3, 2022 at 12:35
• Thank you so much dipetkov. I really appreciate your help. So the outcome variable is being converted to a dichotomous variable. You used the concept of success for your solution. Would this still be the best way to summarize items that represent risks? The items in my dataset don't represent success but whether a participant fulfills a risk category. A value of 1 in the 1st item represents a present risk eg. problem drinking, while 0 would mean there is no risk (no problem drinking). Item 2 would be internet addiction, Item 3 Smoking etc. Do you still recommend this way of handling the data? Commented Oct 3, 2022 at 14:09
• I'm also asking because I initially thought Poisson as a family would be appropriate for such count data. Could you perhaps let me know why you don't use Poisson? Commented Oct 3, 2022 at 14:13