Comparando níveis de fatores após um GLM em R

25

Aqui está um pouco da minha situação: meus dados se referem ao número de presas comidas com sucesso por um predador. Como o número de presas é limitado (25 disponíveis) em cada tentativa, eu tive uma coluna "Amostra" representando o número de presas disponíveis (portanto, 25 em cada tentativa) e outra chamada "Contagem", que foi o número de sucessos ( quantas presas foram comidas). Baseei minha análise no exemplo do livro R em dados de proporção (página 578). As variáveis ​​explicativas são Temperatura (4 níveis, que eu tratei como fator), e Sexo do predador (obviamente, homem ou mulher). Então, acabo com este modelo:

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

Após obter a tabela Análise de desvio, verifica-se que a temperatura e o sexo (mas não a interação) têm um efeito significativo no consumo de presas. Agora, meu problema: preciso saber quais temperaturas diferem, ou seja, tenho que comparar as 4 temperaturas entre si. Se eu tivesse um modelo linear, usaria a função TukeyHSD, mas como estou usando um GLM, não posso. Eu estive olhando o pacote MASS e tentando configurar uma matriz de contraste, mas por algum motivo isso não funciona. Alguma sugestão ou referência?

Aqui está o resumo que recebo do meu modelo, se isso ajudar a torná-lo mais claro ...

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
fonte
2
Oi @ Anne e bem-vindo. Você pode tentar usar a glhtfunção no multcomppacote . Para realizar testes de temperatura do TukeyHSD, use-o assim glht(my.glm, mcp(Temperature="Tukey")). E btw: Sua fórmula modelo pode ser abreviado para: model<-glm(y ~ Temperature*Sex data=predator, family=quasibinomial). Com o asterisco ( ), as interações e os principais efeitos são ajustados.
COOLSerdash
Olá, obrigado pela sua resposta rápida! No entanto, devo estar fazendo algo errado, porque só recebo uma mensagem de erro ... Suponho que my.glm seja a glm que eu executei anteriormente (portanto, "model" no caso). A que mcp se refere? Recebo uma mensagem de erro dizendo que as dimensões dos coeficientes e da matriz de covariância não correspondem ...?
Anne
Seria útil se você editasse sua pergunta e incluísse a saída do modelo.
COOLSerdash
3
Por que você modelou Temperaturecomo um fator? Você não tem os valores numéricos reais? Eu os usaria como uma variável contínua e todo esse problema é discutível.
gung - Restabelece Monica
3
É perfeitamente razoável querer saber como fazer isso em geral; sua pergunta permanece. No entanto, em relação à sua situação específica, eu usaria temp como uma variável contínua, mesmo se você tivesse pensado nisso originalmente como um fator. Deixando de lado os problemas com várias comparações, modelar a temperatura como fator é um uso ineficiente das informações que você possui.
gung - Restabelece Monica

Respostas:

15

Anne, vou explicar rapidamente como fazer comparações múltiplas em geral. Por que isso não funciona no seu caso específico, eu não sei; Eu sinto Muito.

Mas normalmente, você pode fazer isso com o multcomppacote e a função glht. Aqui está um exemplo:

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

Se você deseja calcular as comparações aos pares entre o rankuso do HSD de Tukey, você pode fazer isso desta maneira:

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 o @gung observou nos comentários, você deve - sempre que possível - incluir a temperatura como uma variável contínua e não categórica. Em relação à interação: você pode executar um teste de razão de verossimilhança para verificar se o termo de interação melhora significativamente o ajuste do modelo. No seu caso, o código seria algo assim:

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

Se esse teste não for significativo, você poderá remover a interação do seu modelo. Talvez glhtfuncione então?

COOLSerdash
fonte
11
Oh Deus, muito obrigado !! Consegui escrever o comando corretamente desta vez e funcionou! Obrigado novamente !
Anne
11
Pergunta adicional: existe uma maneira de obter várias comparações sobre a interação? Eu tenho dados semelhantes, onde a interação (a partir da pergunta inicial, que seria Temperatura * Sexo) é significativa, e eu queria saber se é possível comparar esses dados juntos ...
Anne
11
Você quer dizer comparação múltipla para cada nível da interação? Se sim, você pode achar este site interessante (o último parágrafo mostra como testar todas as combinações possíveis em pares).
COOLSerdash
você pode criar uma variável que corresponda às interações de uma variável e usá-la para executar o mcp. Você faz assim. mydata $ gparank <- interação (mydata $ gpa, mydata $ rank)
Notquitesure
11
@Nova qual link você quer dizer? Aquele nos comentários? Aqui está o novo link para esse site.
precisa saber é o seguinte