In another thread, I asked how to obtain a confidence interval of the difference in probability of success between two groups from a logit model and @Weihuang Wong showed me how to do this with a parametric bootstrap. I would like to get the p-value that corresponds to this bootstrap.
library(lme4)
library(pbkrtest)
# Example data:
set.seed(11)
treatment = c(rep("A",30), rep("B", 40))
site = rep(1:14, each = 5)
presence = c(rbinom(30, 1, 0.6),rbinom(40, 1, 0.8))
df = data.frame(presence, treatment, site)
# Likelihood ratio test
M0 = glmer(presence ~ 1 + (1|site), family = "binomial", data = df)
M1 = glmer(presence ~ treatment + (1|site), family = "binomial", data = df)
PBmodcomp(M1,M0, nsim = 500, seed = 11) # LRT-p = 0.3506; PB-p = 0.3936
# 3 samples were not used
195/500 # 0.39
ate <- function(.) {
treat_A <- treat_B <- df
treat_A$treatment <- "A"
treat_B$treatment <- "B"
c("ate" = mean(predict(., newdata = treat_B, type = "response") -
predict(., newdata = treat_A, type = "response")))
}
out <- bootMer(M1, ate, seed = 11, nsim = 500)
quantile(out$t, c(0.025, 0.5, 0.975))
p_value1 <- mean(abs(out$t - mean(out$t))>= abs(mean(out$t))) # p = 0.398
p_value2 <- mean(abs(out$t - mean(out$t))>= abs(out$t0)) # p = 0.356
p_value3 <- mean(abs(out$t - out$t0) >= abs(out$t0)) # p = 0.374
I calculated the p-value in 3 ways: With the first method, I tried to shift the distribution around zero. I haven’t seen a reference for this it just appeared somewhat logical to me before I found references on how to do this.
Thinking about it now, I also think the p-value calculation should include both the sample statistics "t" and the observed statistic "t0".
The second method I obtained from this thread and the third one from this one. (I am not sure though, I understood/did this correctly)
Which one is the correct method? None of them seems to fit the PBmodcomp p-value perfectly even if I removed the correction (i.e. unused samples that yielded negative values in the LRT reference distribution)
Perhaps, more importantly, I also realized in another case (not presented here) that the p-values also do not seem to perfectly match the CI.
Which method do you recommend/ regard as acceptable to calculate a p-value from a bootstrap distribution ?
