Comparación de los niveles de los factores tras un GLM en R

Question

Aquí hay un poco de información sobre mi situación: mis datos se refieren al número de presas comidas con éxito por un depredador. Como el número de presas es limitado (25 disponibles) en cada ensayo, tenía una columna "Muestra" que representaba el número de presas disponibles (por lo tanto, 25 en cada ensayo), y otra llamada "Recuento" que era el número de éxitos (cuántas presas fueron comidas). Basé mi análisis en el ejemplo del libro de R sobre datos de proporción (página 578). Las variables explicativas son Temperatura (4 niveles, que traté como factor), y Sexo del depredador (obviamente, macho o hembra). Así que termino con este modelo:

model <- glm(y ~ Temperature+Sex+Temperature*Sex, 
          data=predator, family=quasibinomial) 

Tras obtener la tabla de análisis de desviación, resulta que la temperatura y el sexo (pero no la interacción) tienen un efecto significativo en el consumo de presas. Ahora, mi problema: necesito saber qué temperaturas difieren, es decir, tengo que comparar las 4 temperaturas entre sí. Si tuviera un modelo lineal, utilizaría el TukeyHSD pero como estoy usando un GLM no puedo. He estado buscando en el paquete MASS y tratando de configurar una matriz de contraste pero por alguna razón no funciona. ¿Alguna sugerencia o referencia?

Este es el resumen que obtengo de mi modelo, si eso ayuda a aclararlo...

y <- cbind(data$Count, data$Sample-data$Count)
model <- glm(y ~ Temperature+Sex+Temperature*Sex, 
             data=predator, family=quasibinomial) 
> summary(model)

# Call:
# glm(formula = y ~ Temperature + Sex + Temperature * Sex, 
       family=quasibinomial, data=data)

# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -3.7926  -1.4308  -0.3098   0.9438   3.6831  

# Coefficients:
#                                        Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                             -1.6094     0.2672  -6.024 3.86e-08 ***
# Temperature8                             0.3438     0.3594   0.957   0.3414    
# Temperature11                           -1.0296     0.4803  -2.144   0.0348 *  
# Temperature15                           -1.2669     0.5174  -2.449   0.0163 *  
# SexMale                                    0.3822     0.3577   1.069   0.2882    
# Temperature8:SexMale                    -0.2152     0.4884  -0.441   0.6606    
# Temperature11:SexMale                    0.4136     0.6093   0.679   0.4990    
# Temperature15:SexMale                    0.4370     0.6503   0.672   0.5033    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

# (Dispersion parameter for quasibinomial family taken to be 2.97372)    
#     Null deviance: 384.54  on 95  degrees of freedom
# Residual deviance: 289.45  on 88  degrees of freedom
# AIC: NA   
# Number of Fisher Scoring iterations: 5

Answer 1

Anne, voy a explicar brevemente cómo hacer esas comparaciones múltiples en general. Por qué no funciona en tu caso concreto, no lo sé; lo siento.

Pero normalmente, se puede hacer con el multcomp y la función glht . He aquí un ejemplo:

mydata      <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
mydata$rank <- factor(mydata$rank)
my.mod      <- glm(admit~gre+gpa*rank, data=mydata, family=quasibinomial)

summary(my.mod)
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)  
# (Intercept) -4.985768   2.498395  -1.996   0.0467 *
# gre          0.002287   0.001110   2.060   0.0400 *
# gpa          1.089088   0.731319   1.489   0.1372  
# rank2        0.503294   2.982966   0.169   0.8661  
# rank3        0.450796   3.266665   0.138   0.8903  
# rank4       -1.508472   4.202000  -0.359   0.7198  
# gpa:rank2   -0.342951   0.864575  -0.397   0.6918  
# gpa:rank3   -0.515245   0.935922  -0.551   0.5823  
# gpa:rank4   -0.009246   1.220757  -0.008   0.9940  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Si quieres calcular las comparaciones por pares entre rank utilizando el HSD de Tukey, se podría hacer de esta manera:

library(multcomp)
summary(glht(my.mod, mcp(rank="Tukey")))
# 
#    Simultaneous Tests for General Linear Hypotheses
# 
# Multiple Comparisons of Means: Tukey Contrasts
# 
# Fit: glm(formula = admit ~ gre + gpa * rank, family = quasibinomial, data = mydata)   
# 
# Linear Hypotheses:
#            Estimate Std. Error z value Pr(>|z|)
# 2 - 1 == 0   0.5033     2.9830   0.169    0.998
# 3 - 1 == 0   0.4508     3.2667   0.138    0.999
# 4 - 1 == 0  -1.5085     4.2020  -0.359    0.984
# 3 - 2 == 0  -0.0525     2.6880  -0.020    1.000
# 4 - 2 == 0  -2.0118     3.7540  -0.536    0.949
# 4 - 3 == 0  -1.9593     3.9972  -0.490    0.960
# (Adjusted p values reported -- single-step method)
# 
# Warning message:
# In mcp2matrix(model, linfct = linfct) :
#   covariate interactions found -- default contrast might be inappropriate

Todas las comparaciones por pares se dan junto con un $p$ -valor. Advertencia: Estas comparaciones sólo se refieren a los efectos principales. Las interacciones se ignoran. Si las interacciones están presentes, se dará una advertencia (como en la salida anterior). Para un tutorial más extenso sobre cómo proceder cuando las interacciones están presentes, vea estos adicionales ejemplos de multcomp .

Nota: Como ha señalado @gung en los comentarios, siempre que sea posible, debería incluir la temperatura como una variable continua en lugar de categórica. En cuanto a la interacción: podría realizar una prueba de razón de verosimilitud para comprobar si el término de interacción mejora significativamente el ajuste del modelo. En tu caso, el código sería algo así:

# Original model
model <- glm(y ~ Temperature+Sex+Temperature*Sex, data=predator, family=quasibinomial) 

# Model without an interaction
model2 <- glm(y ~ Temperature+Sex data=predator, family=quasibinomial) 

# Likelihood ratio test
anova(model, model2, test="LRT")

Si esta prueba no es significativa, puede eliminar la interacción de su modelo. Quizás glht ¿funcionará entonces?