Vergleich der Faktoren nach einem GLM in R

Hier ist ein kleiner Hintergrund zu meiner Situation: Meine Daten beziehen sich auf die Anzahl der Beutetiere, die ein Raubtier erfolgreich gefressen hat. Da die Anzahl der Beutetiere in jedem Versuch begrenzt ist (25 verfügbar), hatte ich eine Spalte "Stichprobe", die die Anzahl der verfügbaren Beutetiere darstellt (also 25 in jedem Versuch), und eine andere Spalte namens "Anzahl", die die Anzahl der Erfolge war ( Wie viele Beute wurden gefressen? Ich habe meine Analyse auf das Beispiel aus dem R-Book für Anteilsdaten (Seite 578) gestützt. Die erklärenden Variablen sind Temperatur (4 Ebenen, die ich als Faktor behandelt habe) und Geschlecht des Raubtiers (offensichtlich männlich oder weiblich). So lande ich mit diesem Modell:

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

Nach dem Abrufen der Abweichungsanalyse-Tabelle hat sich herausgestellt, dass Temperatur und Geschlecht (jedoch nicht die Interaktion) einen signifikanten Einfluss auf den Beutekonsum haben. Nun mein Problem: Ich muss wissen, welche Temperaturen unterschiedlich sind, dh ich muss die 4 Temperaturen miteinander vergleichen. Wenn ich ein lineares Modell hätte, würde ich die TukeyHSD-Funktion verwenden, aber da ich ein GLM verwende, kann ich nicht. Ich habe das Paket MASS durchgesehen und versucht, eine Kontrastmatrix einzurichten, aber aus irgendeinem Grund funktioniert es nicht. Anregungen oder Hinweise?

Hier ist die Zusammenfassung, die ich von meinem Modell bekomme, wenn das hilft, es klarer zu machen ...

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, ich werde kurz erklären, wie solche mehrfachen Vergleiche im Allgemeinen durchgeführt werden. Warum dies in Ihrem speziellen Fall nicht funktioniert, weiß ich nicht; Es tut mir Leid.

Normalerweise können Sie dies jedoch mit dem multcompPaket und der Funktion tun glht. Hier ist ein Beispiel:

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

Wenn Sie die paarweisen Vergleiche zwischen der rankVerwendung von Tukeys HSD berechnen möchten, können Sie dies folgendermaßen tun:

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$

Hinweis: Wie in den Kommentaren unter @gung vermerkt, sollten Sie die Temperatur möglichst als stetige und nicht als kategoriale Variable angeben. In Bezug auf die Interaktion: Sie können einen Likelihood-Ratio-Test durchführen, um zu überprüfen, ob der Interaktionsterm die Modellanpassung signifikant verbessert. In Ihrem Fall würde der Code ungefähr so ​​aussehen:

# 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")

Wenn dieser Test nicht signifikant ist, können Sie die Interaktion aus Ihrem Modell entfernen. Vielleicht klappt glhtdas dann?