I have no explanation for this. Note how "school" isn't significant, yet the model with only that main effect has a much lower AIC (and BIC) than the one with time, intervention, and intervention*time (2 of those are significant).

Response is binary: 0=no, 1=yes

There are three predictors:

Intervention: 0=control, 1=received intervention

School: 4 levels (2 schools got treatment, 2 didn't)

Time: baseline, timepoint 1, timepoint 2

There are ~200 observations per treatment per time point (200*2 treatments*3 times = 1200 total n).

Random effect is SID (student ID#)

Mixed modelling:

model1 <- glmer(y ~ (1 | studentid), 
data = dat, family = binomial, na.action=na.exclude)
AIC: 1314
Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.0759     0.1316  -8.178  2.9e-16 ***

model2 <- glmer(y ~ school.factor +  
              (1 | studentid), data = dat, family = binomial, na.action=na.exclude)
AIC: 462
Fixed effects:
                                    Estimate Std. Error z value Pr(>|z|) 

(Intercept)                             -0.71668    0.26635  -2.691  0.00713 **
school.factorHorace Elementary School   -0.46895    0.43541  -1.077  0.28147  
school.factorLegacy Elementary School    0.01046    0.33193   0.032  0.97486 
school.factorTea Area Legacy Elementary -0.18038    0.32310  -0.558  0.57665  

model3 <- glmer(y ~ intervention + time + time*intervention +
              (1 | studentid), data = dat, family = binomial)
AIC: 1314
Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -1.26999    0.25739  -4.934 8.05e-07 ***
intervention       -0.04564    0.33241  -0.137  0.89080    
time2               0.31367    0.28942   1.084  0.27846    
time3               0.78005    0.29532   2.641  0.00826 ** 
intervention:time2 -0.07410    0.39286  -0.189  0.85040    
intervention:time3 -0.77871    0.39605  -1.966  0.04928 *  

For m3 the AIC is slightly worse if we drop any of the 3 main effects. Also, if I try to include both school.factor and any other main effect variable I get a warning that the fixed-effect matrix is rank deficient so they drop the additional columns/coefficients, always just leaving me with school.factor.

Any idea what's going on here?

  • $\begingroup$ Can you try plotting your data? I suspect that there might be difference due to school ID that are very strong and potentially tower over the treatment effect. For example, school will be much better in capturing socio-economic indicators that might affect a student-associated outcome. $\endgroup$
    – usεr11852
    Aug 7, 2018 at 20:50

1 Answer 1


There are several issues with the data here that I'm going to list out in bullet points. They may or may not directly answer your question, but I feel they are still important:

  • Treating time as continuous may not be appropriate, especially if the timespan between baseline/time1 and time1/time2 are not equal
  • This analysis approach may not be appropriate for this structure of data. The within-school and between-school correlations are not properly accounted for.
  • The error that you get when trying to include 'school.factor' and 'intervention' is likely because there is complete separation and redundancy in these two variables, as intervention is nested within schools. 'schools.factor' is retained because this variable contains more information than 'intervention' (this is my guess).

I don't know this study's design, but my suggestion for you would be to look up design and analyses of cluster randomized trials, as this data structure is very similar to what you'll find in a typical cluster RCT. Here is one paper I found that may be useful for you. This is assuming that you are interested in estimating the effect of intervention, and not building a predictive model.

  • $\begingroup$ I edited the question to show time as a factor, thanks. Indeed I just want to see the effect of intervention on the outcome (I have ~10 outcomes). That paper gives their implementation in STATA, which I have unfortunately never used. $\endgroup$
    – StatsNTats
    Aug 7, 2018 at 21:29
  • $\begingroup$ @StatsNTats I think you should be able to determine the model formulation from the referenced article in the difference in differences section, and easily implement this in R. I am not familiar with the parameterization of a difference-in-differences model myself, but there should be lots of literature on this. $\endgroup$
    – dwhdai
    Aug 8, 2018 at 21:11
  • $\begingroup$ The other issue is they are talking about "randomized" clusters in that difference of difference section. Whereas for my study of only 4 schools (clusters) they purposefully chose the pairs of 2 to be try and be "balanced". I am yet to find anything resembling my situation (it's for sure not an ideal study, but it's the best they could do). $\endgroup$
    – StatsNTats
    Aug 8, 2018 at 21:17
  • $\begingroup$ @StatsNTats That's fine. For the purpose of this analysis, the randomized vs. non-randomized nature is not important. What's important is the nesting structure of the data. Due to its parallels in the structure of data that you would get from a cluster RCT, the analytical methods for a cluster RCT can be used. Where you will need to adapt to the non-randomized nature is in the interpretation of the results, in which case a randomized study would be of higher quality (but in this case with 4 clusters, it's unlikely to be much better anyway) $\endgroup$
    – dwhdai
    Aug 9, 2018 at 14:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.