More ordinal data display

March 30, 2013
By

(This article was first published on Wiekvoet, and kindly contributed to R-bloggers)

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.

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 

 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)

vc <- vcov(Res.clm)
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)
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)
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 on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

Mango solutions



RStudio homepage



Zero Inflated Models and Generalized Linear Mixed Models with R

Quantide: statistical consulting and training

datasociety

http://www.eoda.de





ODSC

ODSC

CRC R books series





Six Sigma Online Training









Contact us if you wish to help support R-bloggers, and place your banner here.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)