How to show that a random factor is not needed in the model?

Question

I collected data from an experiment where I showed one of four videos (condition) to a person and asked them to predict how it ended / assign one of three labels to the video (prediction). I asked each person to make a prediction after they had completed watching 25%, 50%, 75%, and 100% of the video (completed). I collected 330 responses total.

From a more technical perspective, I have a dataset with two independent variables, and one dependent variable. condition is categorical and has 4 levels measured between subject; completion is ordinal with 4 levels (although the underlying construct is continuous - should I model it as a continuous variable?) and measured within-subject. The dependent variable prediction is categorical with 3 levels. I could make it dichotomous, since I have ground truth on the correct prediction, so I could code it as correct/incorrect.

My Hypothesis is that slope of correct predictions over time differs significantly between conditions, i.e., for some videos people make better predictions earlier compared to other videos.

Since this is a repeated measures experiment with a categorical dependent variable, I am using a Generalized Linear Mixed-Effects Model (glmer in R) to fit my data. Here are the models I'm currently fitting:

interaction <- glmer("prediction ~ condition*completion + (1|id)", data=data, family="binomial")
main_effects <- glmer("prediction ~ condition + completion + (1|id)", data=data, family="binomial")
baseline <- glmer("prediction ~ completion + (1|id)", data=data, family="binomial")

I then compare main_effects to baseline and interaction to main_effects using anova(). Am I actually testing my hypothesis like this? (I still get confused easily with linear models)

My main question, however, is: I am getting the error: boundary (singular) fit: see ?isSingular, for main_effects and interaction. Am I interpreting this correctly as: The data doesn't support participants as a random factor, and I should fit a glm instead?

Answer 1

I have fitted a model to your data using the mixed_model function from the GLMMadaptive package without receiving any error or warning, and compared the output to that produced by glmer(with the singular fit):

(data is available here )

> library(lme4)
> library(GLMMadaptive)

> dt <- read.csv("PilotDataStacked-results.csv")

> dt$id <- as.factor(dt$id)
> m0_glmer <- glmer(isCorrect ~ condition + completion + (1|id), data = dt, nAGQ = 20, family="binomial"(link = logit))

> m0_GLMMadaptive <- mixed_model(fixed = isCorrect ~ condition + completion, random = ~ 1 | id, data = dt,
               family = binomial())
> summary(m0_glmer)

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0        0       
Number of obs: 1160, groups:  id, 290

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      2.171244   0.233476   9.300  < 2e-16 ***
conditiongreen  -0.831986   0.222863  -3.733 0.000189 ***
conditiongreen2 -1.139616   0.224589  -5.074 3.89e-07 ***
conditionred    -0.753830   0.215395  -3.500 0.000466 ***
completion      -0.050959   0.003675 -13.867  < 2e-16 ***

> summary(m0_GLMMadaptive) 

Random effects covariance matrix:
               StdDev
(Intercept) 0.1508784

Fixed effects:
                Estimate Std.Err  z-value    p-value
(Intercept)       2.1786  0.2349   9.2745    < 1e-04
conditiongreen   -0.8342  0.2248  -3.7110 0.00020646
conditiongreen2  -1.1428  0.2265  -5.0454    < 1e-04
conditionred     -0.7560  0.2173  -3.4792 0.00050283
completion       -0.0512  0.0037 -13.8441    < 1e-04

As you can see, the results are essentially the same. It is a mystery why glmer obtains a singular fit. The variance of the random intercept estimated by mixed_model is quite small, so that is my only idea, however it is not extremely small, so if I were you I would want to investigate a bit more before proceeding with either model. Ben Bolker is one of the authors of lme4 so he might see this message if I ping him in a comment, which I will do, but if not you might want to ask for input from lme4 specialists on the R-sig-ME mailing list: https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

Also Dimitris Rizopoulos (author of GLMMAdaptive) may have some input. I will also ping him.

Edit: As to your other questions:

completion is ordinal with 4 levels (although the underlying construct is continuous - should I model it as a continuous variable?)

Yes, you are losing information by modelling it as ordinal.

Am I actually testing my hypothesis like this? (I still get confused easily with linear models)

Since your hypothesis is whether the slope of completion changes at different levels of condition then your interest should center on the interaction between them. This will be easier to interpret if you switch to continuous completion.