Porównywanie poziomów czynników po GLM w R

25

Oto kilka informacji na temat mojej sytuacji: moje dane odnoszą się do liczby ofiar, które z powodzeniem zjadły drapieżniki. Ponieważ liczba ofiar jest ograniczona (25 dostępnych) w każdej próbie, miałem kolumnę „Próbka” reprezentującą liczbę dostępnej ofiary (czyli 25 w każdej próbie), a drugą o nazwie „Liczba”, która była liczbą powodzenia ( ile ofiar padło). Swoją analizę oparłem na przykładzie z książki R na danych dotyczących proporcji (strona 578). Zmienne objaśniające to Temperatura (4 poziomy, które traktowałem jako czynnik) i Płeć drapieżnika (oczywiście mężczyzna lub kobieta). Skończyłem z tym modelem:

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

Po uzyskaniu tabeli analizy dewiacji okazuje się, że temperatura i płeć (ale nie interakcja) mają znaczący wpływ na konsumpcję zdobyczy. Teraz mój problem: muszę wiedzieć, które temperatury się różnią, tzn. Muszę porównać 4 temperatury ze sobą. Gdybym miał model liniowy, użyłbym funkcji TukeyHSD, ale ponieważ używam GLM, nie mogę. Przeglądałem pakiet MASS i próbuję skonfigurować matrycę kontrastu, ale z jakiegoś powodu to nie działa. Wszelkie sugestie lub referencje?

Oto podsumowanie, które otrzymuję z mojego modelu, jeśli to pomaga uczynić go bardziej przejrzystym ...

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
źródło
2
glhtmultcompglht(my.glm, mcp(Temperature="Tukey"))model<-glm(y ~ Temperature*Sex data=predator, family=quasibinomial)
Cześć, dziękuję za szybką odpowiedź! Jednak muszę robić coś złego, ponieważ pojawia się tylko komunikat o błędzie ... Zakładam, że my.glm to glm, który wykonałem wcześniej (dlatego „model” w tym przypadku). Do czego odnosi się mcp? Otrzymuję komunikat o błędzie informujący, że wymiary współczynników i macierzy kowariancji nie pasują ...?
Anne
Przydałoby się edytować pytanie i dołączyć wyniki modelu.
COOLSerdash
3
Dlaczego modelowałeś Temperaturejako czynnik? Nie masz rzeczywistych wartości liczbowych? Użyłbym ich jako zmiennej ciągłej, a wtedy cały ten problem jest dyskusyjny.
gung - Przywróć Monikę
3
Zasadniczo rozsądnie jest chcieć wiedzieć, jak to zrobić w ogóle; twoje pytanie pozostaje bez odpowiedzi. Jednak w / uwzględniając twoją konkretną sytuację, użyłbym temp jako zmiennej ciągłej, nawet jeśli początkowo myślałeś o tym jako o czynniku. Odkładając na bok problemy z wieloma porównaniami, temperatura modelowania jako czynnik jest nieefektywnym wykorzystaniem posiadanych informacji.
gung - Przywróć Monikę

Odpowiedzi:

15

Anne, krótko wyjaśnię, jak przeprowadzić takie wielokrotne porównania w ogóle. Dlaczego to nie działa w twoim konkretnym przypadku, nie wiem; Przepraszam.

Ale zwykle można to zrobić za pomocą multcomppakietu i funkcji glht. Oto przykład:

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

Jeśli chcesz obliczyć porównania parami między rankużyciem HSD Tukeya, możesz to zrobić w następujący sposób:

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

Uwaga: Jak zauważył @gung w komentarzach, należy - w miarę możliwości - uwzględniać temperaturę jako zmienną ciągłą, a nie kategoryczną. Odnośnie interakcji: możesz wykonać test współczynnika wiarygodności, aby sprawdzić, czy termin interakcji znacząco poprawia dopasowanie modelu. W twoim przypadku kod wyglądałby mniej więcej tak:

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

Jeśli ten test nie jest znaczący, możesz usunąć interakcję ze swojego modelu. Może glhtwtedy zadziała?

COOLSerdash
źródło
1
O Boże, dziękuję WIĘCEJ !! Tym razem udało mi się poprawnie napisać polecenie i zadziałało! Dzięki jeszcze raz !
Anne
1
Dodatkowe pytanie: czy istnieje sposób na uzyskanie wielu porównań interakcji? Mam podobne dane, w których interakcja (z początkowego pytania, tj. Temperatura * Płeć) jest znacząca i zastanawiałem się, czy można je porównać razem
Anne
1
Czy masz na myśli wielokrotne porównanie dla każdego poziomu interakcji? Jeśli tak, ta witryna może być dla Ciebie interesująca (ostatni akapit pokazuje, jak przetestować wszystkie możliwe kombinacje par).
COOLSerdash
możesz utworzyć zmienną, która odpowiada interakcjom dla zmiennej i użyć tej zmiennej do przeprowadzenia MCCP. Robisz to w ten sposób. mydata $ gparank <- interakcja (mydata $ gpa, mydata $ rank)
Notquitesure
1
@Nova, który link masz na myśli? Ten w komentarzach? Oto nowy link do tej strony.
COOLSerdash