Trust p-values from glmer?

IDEA:

• Use sleepstudy data-structure
• to use the logistic regression consider the binary response is.slow <- Reaction > 300
  
• to simulate under the NULL
• test if the p-values are uniform.

Conclusion that in this setting, the coverage is correct.

# plot sleepstudy data w ggplot

library(ggplot2)
library(lmerTest)
data(sleepstudy)

ggplot(sleepstudy, aes(x=Days, y=Reaction)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE) +
  theme_minimal() +
  labs(title="Reaction time by days of sleep deprivation",
       x="Days",
       y="Reaction time (ms)") + facet_wrap(~Subject, scales="free_y")

sleepstudy$is.slow <- sleepstudy$Reaction > 300

fit <- nullfit <- glmer(is.slow ~ Days + (1|Subject), data=sleepstudy, family=binomial)
lfit <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudy)
# plot(lfit)  # TA-plot

# make nullfit actually null
library(simr)
fixef(nullfit)[] <- 0
sim <- simulate(nullfit, nsim=500)

get_pval <- function(is.slow){
  sleepstudy$is.slow <- is.slow
  fit <- glmer(is.slow ~ Days + (1|Subject), data=sleepstudy, family=binomial)
  coef(summary(fit))["Days", "Pr(>|z|)"]
}

get_pval_nlmer <- function(is.slow){
  sleepstudy$is.slow <- is.slow
  fit <- 
  coef(summary(fit))["Days", "Pr(>|z|)"]
}

pvals <- sapply(sim, get_pval)
summary(pvals<0.05)

   Mode   FALSE    TRUE 
logical     470      30 

ks.test(pvals, "punif", 0, 1)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  pvals
D = 0.02557, p-value = 0.8993
alternative hypothesis: two-sided

# hist(pvals, breaks=20, col="grey", border="white", xlab="p-value", main="Histogram of p-values")