Comparaison des niveaux de facteurs après un GLM dans R

25

Voici un petit historique de ma situation: mes données se réfèrent au nombre de proies mangées avec succès par un prédateur. Comme le nombre de proies est limité (25 disponibles) dans chaque essai, j'avais une colonne "Sample" représentant le nombre de proies disponibles (donc, 25 dans chaque essai), et une autre appelée "Count" qui était le nombre de succès ( combien de proies ont été mangées). J'ai basé mon analyse sur l'exemple du livre R sur les données de proportion (page 578). Les variables explicatives sont la température (4 niveaux, que j'ai traités comme facteur) et le sexe du prédateur (évidemment, homme ou femme). Je me retrouve donc avec ce modèle:

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

Après avoir obtenu le tableau Analyse de la déviance, il s'avère que la température et le sexe (mais pas l'interaction) ont un effet significatif sur la consommation de proies. Maintenant, mon problème: j'ai besoin de savoir quelles températures diffèrent, c'est-à-dire que je dois comparer les 4 températures entre elles. Si j'avais un modèle linéaire, j'utiliserais la fonction TukeyHSD, mais comme j'utilise un GLM, je ne peux pas. J'ai parcouru le package MASS et essayé de configurer une matrice de contraste, mais pour une raison quelconque, cela ne fonctionne pas. Des suggestions ou des références?

Voici le résumé que j'obtiens de mon modèle, si cela aide à le rendre plus clair ...

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
la source
2
glhtmultcompglht(my.glm, mcp(Temperature="Tukey"))model<-glm(y ~ Temperature*Sex data=predator, family=quasibinomial)
Salut, merci pour votre réponse rapide! Cependant, je dois faire quelque chose de mal car je ne reçois qu'un message d'erreur ... Je suppose que my.glm est le glm que j'ai effectué plus tôt (donc, "modèle" dans le cas). À quoi se réfère mcp? Je reçois un message d'erreur indiquant que les dimensions des coefficients et la matrice de covariance ne correspondent pas ...?
Anne
Il serait utile de modifier votre question et d'inclure la sortie du modèle.
COOLSerdash
3
Pourquoi avez-vous modélisé Temperaturecomme facteur? Vous n'avez pas les valeurs numériques réelles? Je les utiliserais comme une variable continue, puis tout ce problème est théorique.
gung - Rétablir Monica
3
Il est parfaitement raisonnable de vouloir savoir comment faire cela en général; votre question demeure. Cependant, en ce qui concerne votre situation spécifique, j'utiliserais temp comme variable continue même si vous l'aviez initialement considéré comme un facteur. Mis à part les problèmes avec des comparaisons multiples, la modélisation de la température en tant que facteur est une utilisation inefficace des informations dont vous disposez.
gung - Rétablir Monica

Réponses:

15

Anne, je vais brièvement expliquer comment faire de telles comparaisons multiples en général. Pourquoi cela ne fonctionne pas dans votre cas spécifique, je ne sais pas; Je suis désolé.

Mais normalement, vous pouvez le faire avec le multcomppackage et la fonction glht. Voici un exemple:

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 vous vouliez calculer les comparaisons par paires entre l' rankutilisation du HSD de Tukey, vous pouvez le faire de cette manière:

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

Remarque: Comme @gung l'a indiqué dans les commentaires, vous devez - dans la mesure du possible - inclure la température en tant que variable continue plutôt que catégorielle. Concernant l'interaction: vous pouvez effectuer un test de rapport de vraisemblance pour vérifier si le terme d'interaction améliore significativement l'ajustement du modèle. Dans votre cas, le code ressemblerait à quelque chose comme ça:

# 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 ce test n'est pas significatif, vous pouvez supprimer l'interaction de votre modèle. Peut glht- être que ça fonctionnera alors?

COOLSerdash
la source
1
Oh mon Dieu, merci beaucoup !! J'ai pu écrire la commande correctement cette fois et cela a fonctionné! Merci encore !
Anne
1
Question supplémentaire: existe-t-il un moyen d'obtenir plusieurs comparaisons sur l'interaction? J'ai des données similaires, où l'interaction (de la question initiale, qui serait Température * Sexe) est significative, et je me demandais s'il était possible de les comparer ensemble ...
Anne
1
Voulez-vous dire une comparaison multiple pour chaque niveau de l'interaction? Si oui, vous pourriez trouver ce site intéressant (le dernier paragraphe montre comment tester toutes les combinaisons possibles par paires).
COOLSerdash
vous pouvez créer une variable qui correspond aux interactions pour une variable et utiliser cette variable pour effectuer le mcp. Vous le faites comme ça. mydata $ gparank <- interaction (mydata $ gpa, mydata $ rank)
Notquitesure
1
@Nova de quel lien parlez-vous? Celui dans les commentaires? Voici le nouveau lien vers ce site.
COOLSerdash