Comparación de niveles de factores después de un GLM en R

Aquí hay algunos antecedentes sobre mi situación: mis datos se refieren al número de presas que un depredador ha comido con éxito. Como la cantidad de presas es limitada (25 disponibles) en cada prueba, tenía una columna "Muestra" que representaba la cantidad de presas disponibles (25 en cada prueba), y otra llamada "Cuenta", que era la cantidad de éxito ( cuántas presas se comieron). Basé mi análisis en el ejemplo del libro 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, hombre o mujer). Así que termino con este modelo:

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

Después de obtener la tabla 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, usaría la función TukeyHSD, pero como estoy usando un GLM no puedo. He estado revisando el paquete MASS e intentando configurar una matriz de contraste, pero por alguna razón no funciona. ¿Alguna sugerencia o referencia?

Aquí está 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

Anne, en breve explicaré cómo hacer comparaciones múltiples en general. Por qué esto no funciona en su caso específico, no lo sé; Lo siento.

Pero normalmente, puede hacerlo con el multcomppaquete y la función glht. Aquí hay un ejemplo:

mydata      <- read.csv("http://www.ats.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 desea calcular las comparaciones por pares entre el rankuso de HSD de Tukey, puede hacerlo 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

$p$

Nota: Como @gung señaló en los comentarios, siempre que sea posible, debe incluir la temperatura como una variable continua en lugar de una categórica. Con respecto a la interacción: puede realizar una prueba de razón de probabilidad para verificar si el término de interacción mejora significativamente el ajuste del modelo. En su caso, el código se vería 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 glhtfuncione entonces?