[This article was first published on Wiekvoet, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
The past two weeks I made a post regarding analyzing ordinal data with R and JAGS. The calculations in the second part made me realize I could actually get top two box intervals out of R. This demonstrated here. For that I needed the inverse logistic transformation. This I pulled out the ARM package, and I now realize the latter also contains a function to analyse ordinal data, which is the second part of this post.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Data
The data is as before the cheese data. I will pick up at calculating the model and printing the model summary. The subsequent calculation is actually fairly simple.
The first parameters in the summary refer to the the difference between cheese A and cheeses B, C and D. The second part are the thresholds between categories for cheese A. So, when I take the threshold of category 7|8 for cheese A (1.54) and add the differences, I get the 7|8 thresholds for the other three cheeses. Parameter values and their variances can be readily obtained. Associated standard deviations follow from the vcov matrix. Take the inverse logit and the desired result is there. It is almost too simple.
Res.clm <- clm(FResponse ~Cheese,data=cheese2)
summary(Res.clm)
formula: FResponse ~ Cheese
data: cheese2
link threshold nobs logLik AIC niter max.grad cond.H
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 208 -355.67 733.35 6(0) 3.14e-11 8.8e+01
Coefficients:
Estimate Std. Error z value Pr(>|z|)
CheeseB -3.3518 0.4287 -7.819 5.34e-15 ***
CheeseC -1.7099 0.3715 -4.603 4.16e-06 ***
CheeseD 1.6128 0.3805 4.238 2.25e-05 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Threshold coefficients:
Estimate Std. Error z value
1|2 -5.46738 0.52363 -10.441
2|3 -4.41219 0.42776 -10.315
3|4 -3.31262 0.37004 -8.952
4|5 -2.24401 0.32674 -6.868
5|6 -0.90776 0.28335 -3.204
6|7 0.04425 0.26457 0.167
7|8 1.54592 0.30168 5.124
8|9 3.10577 0.40573 7.655
co <- coef(Res.clm)[c language=”('7|8','CheeseB','CheeseC','CheeseD')”][/c]
vc <- vcov(Res.clm)[c 1=”style="background-color:” 2=”#f3f3f3;” 3=”-family:” 4=”Courier” 5=”New,” 6=”Courier,” 7=”monospace;” 8=”-size:” 9=”x-small;"> ” 10=” ” 11=” ” 12=” ” 13=” ” 14=” ” 15=” ” 16=” ” 17=” ” 18=” ” 19=”c('7|8','CheeseB','CheeseC','CheeseD')” language=”('7|8','CheeseB','CheeseC','CheeseD'),</span></div><div><span”][/c]
names(co) <- levels(cheese2$Cheese)
sd <- sqrt(c(vc[1,1],diag(vc)[-1]+vc[1,1]-2*vc[1,-1]))
data.frame(
`top 2 box`=arm::invlogit(c(-co[1],-co[1]+co[-1])),
`2.5% limit`=arm::invlogit(c(-co[1],-co[1]+co[-1])+qnorm(.025)*sd),
`97.5% limit`=arm::invlogit(c(-co[1],-co[1]+co[-1])+qnorm(.975)*sd),
check.names=FALSE
)
top 2 box 2.5% limit 97.5% limit
A 0.175676950 0.105533878 0.27795336
B 0.007407959 0.003235069 0.01687231
C 0.037118760 0.018617885 0.07264338
D 0.516712572 0.384663489 0.64646825
ARM
ARM (Data Analysis Using Regression and Multilevel/Hierarchical Models) is the package associated with the similar titled book (Gelman & Hill) which is certainly value for money. Having said that, I did not remember this being in the book. The package also contains a simulation function, which I used to get some posterior results for further processing.
library(arm)
library(arm)
Res.arm <- bayespolr(FResponse ~Cheese,data=cheese2)
Res.arm
bayespolr(formula = FResponse ~ Cheese, data = cheese2)
coef.est coef.se
CheeseB -3.25 0.42
CheeseC -1.63 0.37
CheeseD 1.59 0.38
1|2 -5.36 0.52
2|3 -4.32 0.42
3|4 -3.23 0.36
4|5 -2.18 0.32
5|6 -0.86 0.28
6|7 0.07 0.26
7|8 1.55 0.30
8|9 3.09 0.40
—
n = 208, k = 11 (including 8 intercepts)
residual deviance = 727.1, null deviance is not computed by polr
sims <- sim(Res.arm,n.sims=1000)
cosims <- coef(sims)
mycoef <- cbind(-cosims[,’7|8′],
-cosims[,’7|8′]+cosims[,grep(‘Cheese’,colnames(cosims),value=TRUE)])
mycoef <- invlogit(mycoef)
mycoef <- cbind(-cosims[,’7|8′],
-cosims[,’7|8′]+cosims[,grep(‘Cheese’,colnames(cosims),value=TRUE)])
mycoef <- invlogit(mycoef)
coefplot(rev(apply(mycoef,2,mean)),rev(apply(mycoef,2,sd)),
varnames=rev(levels(cheese2$Cheese)),
main=’Estimated Proportion Top 2 Box’)
To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.