I am having trouble accounting for interdependence in data that I would like to model with binomial regression. I am looking at school attendance during two months (43 school days), counting whether a student attended the school on a given day or not (0/1). I would like to model these data with individual-level predictors (e.g. socio-economic status, religiosity, or sex), and I am not interested in predicting attendance on a specific day. This looks like an ideal situation for aggregated binomial regression with 43 trials and binary success/failure outcomes.
For example, in R, I could fit:
glm(cbind(attend, not_attend) ~ sex, family = 'binomial')
where 'attend' is the number of days a student attended the school while 'not_attend' is the number of days the student were not at school. Note that I need to use ratio (instead of counts) because some students were recorded only for 40 days.
However, I am not sure how to take into account the fact that the 43 trials are not independent. E.g. if a student has fever and fails to attend the school one day, it's more likely that the student will not attend also the following day. I considered using beta-binomial regression, but the probability of attending is not random, it is interdependent (although see Ben Bolker's answer here).
The only reasonable solution I could come up with was to nest the data within individual using a multi-level model. For instance:
glmer(cbind(attend, not_attend) ~ sex + (1|id), family = 'binomial')
where 'id' is a student's identification number. If I understand the model correctly, 'sex' models the varying intercept for each individual rather than the probability of attendance on a particular day. This is getting close to what I would like to do, but I am worried that this approach still does not account for the temporal interdependence within the data.
This looks like a simple question, yet I am failing to find a satisfactory answer.