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?

  • $\begingroup$ You don't seem to have time included in your model. It may not fix the singular fit problem, bur I would have thought that time plays a part in the data generation process. Have you tried including it as a fixed effect ? Moreover, does the response (correct/incorrect) vary for each participant ? $\endgroup$ – Robert Long Feb 5 '20 at 8:18
  • $\begingroup$ @RobertLong completion is time in this case. Initially a person observes the first 25% of the video and makes a prediction, then the same person observes the first 50%, makes a prediction, and so forth. Also yes, the response varies between and within subject for each level of the independent variables (some change their response as more of the video is revealed), with exception for the 100% level. Predicting the end after you have observed the end doesn't seem to cause much variance ... $\endgroup$ – FirefoxMetzger Feb 5 '20 at 8:29
  • $\begingroup$ Ahh OK ! So are you using the 2 level or 3 level variable as the response ? It is difficult to say much more without actually seeing the data. Can you share it ? $\endgroup$ – Robert Long Feb 5 '20 at 8:34
  • $\begingroup$ @RobertLong Yes I have two 4 level variables as inputs and one 3 level variable as response. I think it will be okay to share the data, but I will have to double-check with my supervisor :) How would I share it? $\endgroup$ – FirefoxMetzger Feb 5 '20 at 8:43
  • $\begingroup$ Wait, so the problems with singularity are with a 3-level response variable ? This may be causing the problems - that would be a multinomial model and glmer doesn't fit those.. What happens when you convert it to binary. $\endgroup$ – Robert Long Feb 5 '20 at 9:00

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:
(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.

  • $\begingroup$ @DimitrisRizopoulos do you have any insight ? $\endgroup$ – Robert Long Feb 5 '20 at 13:58
  • $\begingroup$ @BenBolker do you have any insight about why glmer is having trouble with these data ? $\endgroup$ – Robert Long Feb 5 '20 at 13:59
  • $\begingroup$ Perhaps worth considering the blme package as well? It imposes a reasonable prior on the random intercept that pushes it away from 0. cran.r-project.org/web/packages/blme $\endgroup$ – Erik Ruzek Feb 5 '20 at 17:31
  • $\begingroup$ Another question: Does @FirefoxMetzger actually have the continuous version of the outcome or are they saying that there is an underlying (latent) continuous distribution, suggesting that the generalized mixed model approach is correct? $\endgroup$ – Erik Ruzek Feb 5 '20 at 17:33
  • 1
    $\begingroup$ A very quick exploration shows that glmer gets a fit with log-likelihood 0.32 units greater than GLMMadaptive's fit - a suggestion that there may (?) be multiple modes and glmer may (?) be finding the better one. Next thing to try would be to start each method at the other's optimum. (This assumes that the log-likelihood computations themselves match.) $\endgroup$ – Ben Bolker Feb 6 '20 at 15:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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