Wie Profil Vertrauensintervalle der Unterschied in der Wahrscheinlichkeit des Erfolgs zwischen zwei Gruppen von einem Logit-Modell (glmer) erhalten?

Question

Ich habe Schwierigkeiten, die Log-Odds-Ratio-Profil-Konfidenzintervalle aus einem Logit-Modell in Wahrscheinlichkeiten zu transformieren. Ich würde gerne wissen, wie man die Konfidenzintervalle der Differenz zwischen zwei Gruppen berechnet.

Wenn der p-Wert> 0,05 ist, sollte das 95% CI der Differenz von unter Null bis über Null reichen. Ich weiß jedoch nicht, wie negative Werte erhalten werden können, wenn die logarithmischen Verhältnisse potenziert werden müssen. Deshalb habe ich versucht, das CI einer der Gruppen (B) zu berechnen und zu sehen, was der Unterschied zwischen dem unteren und dem oberen Ende des CI zur Schätzung der Gruppe A ist. Ich glaube, dass dies nicht der richtige Weg ist, das CI der Differenz zu berechnen, da die Schätzung von A ebenfalls unsicher ist.

Ich würde mich freuen, wenn mir jemand helfen könnte.

library(lme4)  
# 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) 
anova(M1, M0) 

# Calculating confidence intervals 
cc <- confint(M1, parm = "beta_") 
ctab <- cbind(est = fixef(M1), cc) 
cdat = as.data.frame(ctab) 

# Function to back-transform to probability (0-1) 
unlogit = function(y){ 
    y_retransfromed = exp(y)/(1+exp(y)) 
    y_retransfromed 
} 

# Getting estimates 
A_est = unlogit(cdat$est[1]) 
B_est = unlogit(cdat$est[1] + cdat$est[2]) 
B_lwr = unlogit(cdat$est[1] + cdat[2,2]) 
B_upr = unlogit(cdat$est[1] + cdat[2,3]) 

Difference_est = B_est - A_est 

# This is how I tried to calculate the CI of the difference 
Difference_lwr = B_lwr - A_est 
Difference_upr = B_upr - A_est 

# However, I believe this is wrong because A_est is also “uncertain” 

Wie das Konfidenzintervall der Differenz der Wahrscheinlichkeit des Vorhandenseins bekommen?

Answer 1

Wir können den durchschnittlichen Behandlungseffekt wie folgt berechnen. Erstellen Sie aus den Originaldaten zwei neue Datensätze, einen, in dem alle Einheiten die Behandlung A erhalten, und einen, in dem alle Einheiten die Behandlung B erhalten. Basierend auf Ihren Modellschätzungen (in Ihrem Fall M1) berechnen wir vorhergesagte Ergebnisse für Einheiten in jedem dieser zwei Datensätze. Wir berechnen dann den mittleren Unterschied in den Ergebnissen zwischen den beiden Datensätzen, um unseren geschätzten durchschnittlichen Behandlungseffekt zu erhalten. Hier können wir eine Funktion schreiben, die ein glmer Objekt nimmt und berechnet den durchschnittlichen Behandlungseffekt:

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"))) 
} 
ate(M1) 
#  ate 
# 0.09478276 

Wie wir das Unsicherheitsintervall bekommen? Wir können den Bootstrap verwenden, d. H. Das Modell viele Male neu berechnen, indem wir zufällig erzeugte Stichproben aus Ihren Originaldaten verwenden, wobei wir jedes Mal den durchschnittlichen Behandlungseffekt berechnen. Wir können dann die Verteilung der durchschnittlichen Bootstrapper-Behandlungseffekte verwenden, um unser Unsicherheitsintervall zu berechnen. Hier erzeugen wir 100 Simulationen mit der bootMer Funktion

out <- bootMer(M1, ate, seed = 1234, nsim = 100) 

und untersuchen die Verteilung der Wirkung:

quantile(out$t, c(0.025, 0.5, 0.975)) 
#  2.5%   50%  97.5% 
# -0.06761338 0.10508751 0.26907504