2015-10-24 2 views
15

I équipés d'un modèle logit cumulatif logit sur les données ordinales en utilisant la fonction de polrMASS en utilisant (dans ce cas sur des données donnant la préférence pour les différents types de fromage):R: prédictions tracé de modèle ordinal MASSE POLR

data=read.csv("https://www.dropbox.com/s/psj74dx8ohnrdlp/cheese.csv?dl=1") 
data$response=factor(data$response, ordered=T) # make response into ordered factor 
head(data) 
    cheese response count 
1  A  1  0 
2  A  2  0 
3  A  3  1 
4  A  4  7 
5  A  5  8 
6  A  6  8 
library(MASS) 
fit=polr(response ~ cheese, weights=count, data=data, Hess=TRUE, method="logistic") 

Pour tracer les prédictions du modèle que je fait un tracé d'effet à l'aide

library(effects) 
library(colorRamps) 
plot(allEffects(fit),ylab="Response",type="probability",style="stacked",colors=colorRampPalette(c("white","red"))(9)) 

enter image description here

Je me demandais si, d'après les moyennes prédites rapportées par le paquet effects, on pourrait aussi tracer quelque chose comme la préférence moyenne pour chaque type de fromage avec les intervalles de conf 95% sur cela?

EDIT: à l'origine j'ai aussi posé des questions sur la façon d'obtenir des tests Tukey de POSTHOC, mais en attendant, je trouve que ceux-ci peuvent être obtenus en utilisant

library(multcomp) 
summary(glht(fit, mcp(cheese = "Tukey"))) 

ou en utilisant package lsmeans comme

summary(lsmeans(fit, pairwise ~ cheese, adjust="tukey", mode = "linear.predictor"),type="response") 
+0

question intéressante. Je suppose (comme vous le faites) que le problème se pose parce que vous prenez les moyens * après * vous avez créé les probabilités prédites. Voir [ici] (http://stackoverflow.com/questions/14390080/linear-predictor-ordered-probit-ordinal-clm) et [ici] (http://stats.stackexchange.com/questions/41006/predicting- ordered-logit-in-r/41025 # 41025) pour plus d'informations à ce sujet sur SE. En outre, avec 9 catégories, je choisirais simplement une OLS sur la variable de réponse qui produit presque exactement les mêmes estimations de points pour les catégories moyennes, ensemble avec des erreurs standard sensibles. Mais c'est une question intéressante. – Felix

+0

Oui, je pense que cela a à voir avec la moyenne sur l'échelle logit cumulative par rapport à l'échelle finale rétrotransformée. Donc, fondamentalement, je voudrais savoir comment faire la moyenne sur l'échelle des liens et ensuite rétrocalculer à l'échelle ordinale originale. Je sais que pour 9 catégories, je pourrais aussi faire des OLS, mais je voudrais aussi une solution générale pour moins de catégories, par ex. 3 ou 4. –

+0

Les graphiques de dynamite (ces graphiques en barres) sont simplement de mauvaises statistiques. Vous ne gagnez pas plus d'informations que vous ne le faites à partir du tableau 'wmeans' des statistiques récapitulatives. et en raison du fait que ceci n'est qu'une parcelle de statistiques sommaires, vous perdez toutes les données qui ont été utilisées pour le faire. les parcelles devraient montrer des données et non des statistiques sommaires. Je pense que cela résout votre problème puisque vous ne devriez pas le faire en premier lieu – rawr

Répondre

2

Russ La dixième m'a fait remarquer que la préférence moyenne et les intervalles de confiance de 95% peuvent être obtenus simplement avec lsmeans avec option mode="mean" (?models) (il en est de même pour un clm ou clmm ajustement du modèle en utilisant package ordinal):

df=summary(lsmeans(fit, pairwise ~ cheese, mode = "mean"),type="response")$lsmeans 
cheese mean.class  SE df asymp.LCL asymp.UCL 
A  6.272828 0.1963144 NA 5.888058 6.657597 
B  3.494899 0.2116926 NA 3.079989 3.909809 
C  4.879440 0.2006915 NA 4.486091 5.272788 
D  7.422159 0.1654718 NA 7.097840 7.746478 

qui me donne l'intrigue que je cherchais:

library(ggplot2) 
library(ggthemes) 
ggplot(df, aes(cheese, mean.class)) + geom_bar(stat="identity", position="dodge", fill="steelblue", width=0.6) + 
    geom_errorbar(aes(ymin=asymp.LCL, ymax=asymp.UCL), width=0.15, position=position_dodge(width=0.9)) + 
    theme_few(base_size=18) + xlab("Type of cheese") + ylab("Mean preference") + 
    coord_cartesian(ylim = c(1, 9)) + scale_y_continuous(breaks=1:9) 

enter image description here