(This article was first published on

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.**Wiekvoet**, and kindly contributed to R-bloggers)### 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="font-family:" 4="Courier" 5="New," 6="Courier," 7="monospace;" 8="font-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 his blog:**Wiekvoet**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...