2
$\begingroup$

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 ?

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.