Perfect separation, perhaps? In binary outcome and repeated measure (random effect) with multiple independent variables (using R)

Question

(Using R) - this is my first time posting a stats question online, so please let me know if I'm on the wrong forum or haven't provided enough information and I'll do my best to fix it!

About the data and my goal here: Best analogy I can think of is that it's a language course and the final exam is a long conversation. Four times during the course I gather reports on student performance (for example, handwriting, speed of writing, reading ability). I want to know if I can predict pass or failure for the course based on these four reports. I've created a demo dataset here:

set.seed(22)
reportsdata <- structure(list(Student = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
                       3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 
                       6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 
                       9L, 9L, 9L, 9L), 
           TermReport = c("A", "B", "C", "D", "A", "B", "C", "D",
                     "A", "B", "C", "D", "A", "B", "C", "D", 
                     "A", "B", "C", "D",  "A", "B", "C", "D", 
                     "A", "B", "C", "D",  "A", "B", "C", "D", 
                     "A", "B", "C", "D"), 
           Handwriting = c(sample(x = 1:5, size = 36, replace = TRUE)), 
           Speedwriting = c(sample(x = 1:5, size = 36, replace = TRUE)), 
           Reading = c(sample(x = 1:5, size = 36, replace = TRUE)), 
           Loudness = c(sample(x = 1:5, size = 36, replace = TRUE)),  
           Enthusiasm = c("5", "5", "3", "5", "2", "4", "3", "NA", "1", "4",
                          "3", "3", "NA", "2", "1", "1", "1", "2", "2", "NA",
                          "3", "2", "4", "2", "4", "3","5", "2", "3", "1", 
                          "2", "3", "5", "4", "NA", "5"),  
           EndCoursePassFail = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
                                 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
                                 0L, 0L, 0L, 0L)), 
      class = "data.frame", row.names = c(NA, -36L))

Note that their score at the end of the course has been retroactively applied to all their entries, though whether they would pass (1) or fail (0) was not known at the time. My real dataset has the same structure, but contains a little over 600 observations and 30 variables (which have been filtered to those that contain less than 30% NA entries, e.g. sometimes could not get a score for enthusiasm).

So far I've been trying a mixed effects logistic regression with student and trial as random effects (Bobyqa & Nelder_Mead are the only optimisers that don't fail, I need to use ~ . syntax as there are too many variables to list and for reproducibility). E.g.:

model <- glmer(EndCoursePassFail 
                ~ . -Student -TermReport + (1|Student) + (1|TermReport),
                data = reportsdata, 
                family = binomial, 
                control = glmerControl
                (optimizer = "bobyqa", optCtrl=list(maxfun=1e6)), 
                nAGQ = 1, 
                na.action=na.exclude)

For both my original dataset and for the sample data provided above when the seed is set to 22, my model produces convergence errors:

Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : 
Model failed to converge with max|grad| = 0.0477113 (tol = 0.002, component 1)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
Model is nearly unidentifiable: very large eigenvalue - Rescale variables?

But, my sample dataset shows the following errors if seed is set to 1:

boundary (singular) fit: see help('isSingular')

I think my issue could be due to perfect separation between Student and Course Result since the result had been retroactively added.

The question is - where to go from here?

Some thoughts:

1. I could average student scores across term reports somehow, so that I no longer need repeated measures and therefore don't get perfect separation. But this seems crude and feels like looking at only the tip of the iceberg.

2. Looking at other answers for similar issues, I might need to switch to trying a penalised likelihood from blme (R package), however I don't understand it well enough yet to know whether (and if so, how) this sort of perfect separation can be dealt with using blme.

3. Or, I could pretend that there aren't any repeated measures and run the model as though there weren't any - but of course, this is also crude and ignores a lot of potentially useful information provided by the data.

Also, in case it is relevant - because there are so many scores to include in the full dataset, I want to later use stepAIC (or a loop, or equivalent) to roughly identify the 'best' model.

Answer 1

First, don't treat TermReport as a random effect. The isSingular warning comes from finding 0 variance for the attempt by glmer() to force the intercepts for TermReport into a Gaussian distribution. In many cases with 5 or fewer levels of a categorical predictor it's better to model as a fixed effect. When I replaced the random effect for TermReport with a fixed effect I had no errors or warnings with your sample data.

Second, be careful how you handle missing data. For your "Enthusiasm" you have specified "NA" as an actual level. That's not usually a wise choice. You might instead consider multiple imputation to deal with the missing data. See Section 1.3 of the linked reference for the problems with solutions like yours.

Third, watch out for overfitting with your data. You typically need about 15 cases in the minority outcome class per coefficient that you are evaluating. Even if half of the 600 students failed (300 in the minority outcome group), that would limit you to about 20 coefficients. I doubt that such a high proportion failed. If you had a continuous outcome (like a course grade) instead of pass/fail you could probably fit a more flexible model.

So 30 variables are already too many, and multi-level categorical predictors make it worse. For example, the way that you have coded Enthusiasm requires 5 coefficients, and would still require 4 if you did multiple imputation to get around the NA category.

Fourth, do NOT use stepAIC() or similar automated methods to find the "best" model. Yes, that cuts down on the number of predictors, but it does it in a way that is likely to lead you astray, to end up overfitting your data anyway, and to make inference and p-values uninterpretable. Any use of the outcomes to choose the model ends up posing problems. See this page among many others on this site.

Fifth, you will be better off using the data and your understanding of the subject matter to reduce the number of predictors in a way that doesn't involve the outcomes. Highly related predictors can be combined into single combined predictors. Continuous predictors might be reduced to a few principal components. See Chapter 4 of Harrell's Regression Modeling Strategies, in particular Section 4.7 dealing with "data reduction."

Sixth, penalized regression is certainly a possibility if you do end up with perfect separation. In the form of ridge regression, it also allows you to use all of your predictors in a way that, properly implemented, will minimize overfitting by reducing the magnitudes of all the regression coefficients. It's worth learning about in general even if you end up not needing it for this project.

Answer 2

You get a perfect separation because the random effects can seperate the results for each student. For n students you have 4n datarows, but within the students the result is each time the same.

So effectively you have no repetition within the groups for whih you compute the random effects.