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!
 A: 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"
  ) %>%
  read_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)

