Salida inexacta de los coeficientes cuando se utiliza el glm binomial en R

Question

Tengo un conjunto de datos que tiene 8 tratamientos, junto con los datos de éxito. Utilicé un glm binomial para analizar los datos, pero generó algunos valores de coeficiente inesperados para algunos de los tratamientos y no estoy seguro de qué hacer al respecto.

Aquí están los datos y el resumen:

treatment = as.factor(c("A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "C", "C", 
"D", "D", "D", "D", "E", "E", "E", "E", "F", "F", "F", "F", "G", 
"G", "G", "G", "H", "H", "H", "H"))
rep = c(1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 
1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4)
success = c(1, 1, 1, 2, 14, 17, 15, 18, 0, 0, 0, 0, 18, 18, 17, 18, 4, 
4, 2, 4, 2, 4, 1, 1, 1, 0, 0, 1, 8, 6, 6, 2)
total = c(20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 
20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 
20)

data = data.frame(treatment,rep,success,total)
data$perc = data$success/data$total

library(tidyverse)
data %>% group_by(treatment) %>% summarize(mean = mean(perc))

Podemos ver que la media de los tratamientos B y D es de 0,8 y 0,88 respectivamente.

Ahora realiza un glm:

model.glm = glm(cbind(success,total) ~ treatment-1,data = data,family="binomial")

logit2prob <- function(logit){
  odds <- exp(logit)
  prob <- odds / (1 + odds)
  return(prob)
}

SuccessProb = logit2prob(coef(model.glm))
SuccessProb = round(logit2prob(coef(model.glm)),2)
SuccessProb

Podemos ver que las estimaciones, utilizando un glm, para B y D son de 0,44 y 0,47, respectivamente. Estas no se acercan a las estimaciones de resumen.

Si utilizamos un anova, los resultados son mejores.

model.aov = aov(perc ~ treatment-1,data=data)
SuccessProb.aov = coef(model.aov)
SuccessProb.aov

Aquí, las estimaciones para B y D son de 0,8 y 0,89. Mucho mejor que el glm.

¿Alguien sabe qué estoy haciendo mal?

Answer 1

Puedes estar confundiendo cómo R glm() funciona con la forma en que SAS PROC LOGISTIC funciona. En R, cuando se tienen datos binomiales con $>1$ En el ensayo de Bernoulli, hay que utilizar las probabilidades (eventos a no eventos), no la probabilidad (eventos a total). También es posible utilizar los eventos/total con el total como weights argumento. Considera:

agg.mean  = aggregate(perc                   ~treatment,   data, mean)
model.glm = glm(cbind(success,total)         ~treatment-1, data, family="binomial")
mod.glm.o = glm(cbind(success,total-success) ~treatment-1, data, family="binomial")
mod.glm.w = glm(      success/total          ~treatment-1, data, family="binomial",
                                                           weights=data$total)
output     = data.frame(agg.mean, 
                        round(logit2prob(coef(mod.glm.o)),4),
                        round(logit2prob(coef(mod.glm.w)),4),
                        round(logit2prob(coef(model.glm)),4) )
names(output)[3:5] = c("odds type", "weighted % type", "prob type")
output
#            treatment   perc odds type weighted % type prob type
# treatmentA         A 0.0625    0.0625          0.0625    0.0588
# treatmentB         B 0.8000    0.8000          0.8000    0.4444
# treatmentC         C 0.0000    0.0000          0.0000    0.0000
# treatmentD         D 0.8875    0.8875          0.8875    0.4702
# treatmentE         E 0.1750    0.1750          0.1750    0.1489
# treatmentF         F 0.1000    0.1000          0.1000    0.0909
# treatmentG         G 0.0250    0.0250          0.0250    0.0244
# treatmentH         H 0.2750    0.2750          0.2750    0.2157

Tal vez le ayude leer mi respuesta aquí: Diferencia en la salida entre proc genmod de SAS y glm de R .