Ordered categorical data analysis with a study from experimental economics, on the topic of ``storable votes''. See Chapter 15 in Regression and Other Stories.


Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstan")
rstan_options(auto_write = TRUE)
library("rstanarm")
invlogit<-plogis

Load data

data_2player <- read.csv(root("Storable/data","2playergames.csv"))
data_3player <- read.csv(root("Storable/data","3playergames.csv"))
data_6player <- read.csv(root("Storable/data","6playergames.csv"))
data_all <- rbind(data_2player, data_3player, data_6player)
data_all$factor_vote <- factor(data_all$vote, levels = c(1, 2, 3), labels = c("1", "2", "3"), ordered=TRUE)
head(data_all)
  school person round proposal value vote cutoff.12 cutoff.23 sd.logit
1      1    301     1        1    61    2    40.145    71.465   3.1555
2      1    301     2        1    43    2    40.145    71.465   3.1555
3      1    301     3        1     6    1    40.145    71.465   3.1555
4      1    301     4        1    11    1    40.145    71.465   3.1555
5      1    301     5        1    99    3    40.145    71.465   3.1555
6      1    301     6        1    13    1    40.145    71.465   3.1555
  factor_vote
1           2
2           2
3           1
4           1
5           3
6           1

Simple analysis using data from just one person

data_401 <- subset(data_2player, person == 401, select = c("vote", "value"))
data_401$factor_vote <- factor(data_401$vote, levels = c(1, 2, 3), labels = c("1", "2", "3"), ordered=TRUE)
head(data_401)
    vote value factor_vote
300    2    40           2
301    1     4           1
302    3    44           3
303    1    23           1
304    3    99           3
305    2    82           2
fit_1 <- stan_polr(factor_vote ~ value, data = data_401,
                   prior = R2(0.3, "mean"), refresh = 0)
print(fit_1, digits=2)
stan_polr
 family:       ordered [logistic]
 formula:      factor_vote ~ value
 observations: 20
------
      Median MAD_SD
value 0.09   0.03  

Cutpoints:
    Median MAD_SD
1|2 2.70   1.43  
2|3 5.84   2.19  

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

6 people

plotted <- c(101, 303, 409, 405, 504, 112)
story <- c("Perfectly monotonic",
           "One fuzzy and one sharp cutpoint",
           "Monotonic with one outlier",
           "Only 1's and 3's",
           "Almost only 3's",
           "Erratic")
n_plotted <- length(plotted)
data <- as.list(rep(NA, n_plotted))
fit <- as.list(rep(NA, n_plotted))
for (i in 1:n_plotted){
  #ok <- data_all[,"person"]==plotted[i]
  data[[i]] <- subset(data_all, person == plotted[i], 
                      select = c("vote", "factor_vote", "value"))
    fit[[i]] <- stan_polr(factor_vote ~ value, data=data[[i]],
                          prior=R2(0.3, "mean"), refresh = 0,
                          cores = 1, open_progress = FALSE,
                          adapt_delta = 0.9999)
}

Graph

par(mfrow=c(2,3), mgp=c(1.5,.5,0), tck=-.01)
for (i in 1:n_plotted){
  sims <- as.matrix(fit[[i]])
  n_cutpoints <- 2
  cutpoints <- rep(NA, n_cutpoints)
  for (i_cut in 1:n_cutpoints){
    cutpoints[i_cut] <- median(sims[,i_cut+1]/sims[,1])
  }
  s <- median(1/sims[,1])
  plot(data[[i]][,"value"], data[[i]][,"vote"], xlim=c(0,100), ylim=c(1,3),
        xlab="Value", ylab="Vote", main=story[i], yaxt="n")
  axis (2, 1:(n_cutpoints+1))
  temp <- seq(0, 100, 0.1)
  prob <- array(NA, c(length(temp), n_cutpoints+1))
  expected <- rep(NA, length(temp))
  prob[,1] <- 1 - invlogit((temp-cutpoints[1])/s)
  expected <- 1*prob[,1]
  for (i_cut in 2:n_cutpoints){
    prob[,i_cut] <- invlogit((temp-cutpoints[i_cut-1])/s) -
      invlogit((temp-cutpoints[i_cut])/s)
    expected <- expected + i_cut*prob[,i_cut]
  }
  prob[,n_cutpoints+1] <- invlogit((temp-cutpoints[n_cutpoints])/s)
  expected <- expected + (n_cutpoints+1)*prob[,n_cutpoints+1]
  lines (temp, expected, lwd=.5)
  for (i_cut in 1:n_cutpoints){
    lines(rep(cutpoints[i_cut],2), i_cut+c(0,1), lwd=.5)
  }
}

LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBhbmQgT3RoZXIgU3RvcmllczogU3RvcmFibGUiCmF1dGhvcjogIkFuZHJldyBHZWxtYW4sIEplbm5pZmVyIEhpbGwsIEFraSBWZWh0YXJpIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIHRvYzogdHJ1ZQogICAgdG9jX2RlcHRoOiAyCiAgICB0b2NfZmxvYXQ6IHRydWUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCk9yZGVyZWQgY2F0ZWdvcmljYWwgZGF0YSBhbmFseXNpcyB3aXRoIGEgc3R1ZHkgZnJvbSBleHBlcmltZW50YWwKZWNvbm9taWNzLCBvbiB0aGUgdG9waWMgb2YgYGBzdG9yYWJsZSB2b3RlcycnLiBTZWUgQ2hhcHRlciAxNSBpbgpSZWdyZXNzaW9uIGFuZCBPdGhlciBTdG9yaWVzLgoKLS0tLS0tLS0tLS0tLQoKCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQobWVzc2FnZT1GQUxTRSwgZXJyb3I9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGNvbW1lbnQ9TkEpCiMgc3dpdGNoIHRoaXMgdG8gVFJVRSB0byBzYXZlIGZpZ3VyZXMgaW4gc2VwYXJhdGUgZmlsZXMKc2F2ZWZpZ3MgPC0gRkFMU0UKYGBgCgojIyMjIExvYWQgcGFja2FnZXMKCmBgYHtyIH0KbGlicmFyeSgicnByb2pyb290IikKcm9vdDwtaGFzX2ZpbGUoIi5ST1MtRXhhbXBsZXMtcm9vdCIpJG1ha2VfZml4X2ZpbGUoKQpsaWJyYXJ5KCJyc3RhbiIpCnJzdGFuX29wdGlvbnMoYXV0b193cml0ZSA9IFRSVUUpCmxpYnJhcnkoInJzdGFuYXJtIikKaW52bG9naXQ8LXBsb2dpcwpgYGAKCiMjIyMgTG9hZCBkYXRhCgpgYGB7ciB9CmRhdGFfMnBsYXllciA8LSByZWFkLmNzdihyb290KCJTdG9yYWJsZS9kYXRhIiwiMnBsYXllcmdhbWVzLmNzdiIpKQpkYXRhXzNwbGF5ZXIgPC0gcmVhZC5jc3Yocm9vdCgiU3RvcmFibGUvZGF0YSIsIjNwbGF5ZXJnYW1lcy5jc3YiKSkKZGF0YV82cGxheWVyIDwtIHJlYWQuY3N2KHJvb3QoIlN0b3JhYmxlL2RhdGEiLCI2cGxheWVyZ2FtZXMuY3N2IikpCmRhdGFfYWxsIDwtIHJiaW5kKGRhdGFfMnBsYXllciwgZGF0YV8zcGxheWVyLCBkYXRhXzZwbGF5ZXIpCmRhdGFfYWxsJGZhY3Rvcl92b3RlIDwtIGZhY3RvcihkYXRhX2FsbCR2b3RlLCBsZXZlbHMgPSBjKDEsIDIsIDMpLCBsYWJlbHMgPSBjKCIxIiwgIjIiLCAiMyIpLCBvcmRlcmVkPVRSVUUpCmhlYWQoZGF0YV9hbGwpCmBgYAoKIyMjIyBTaW1wbGUgYW5hbHlzaXMgdXNpbmcgZGF0YSBmcm9tIGp1c3Qgb25lIHBlcnNvbgoKYGBge3IgfQpkYXRhXzQwMSA8LSBzdWJzZXQoZGF0YV8ycGxheWVyLCBwZXJzb24gPT0gNDAxLCBzZWxlY3QgPSBjKCJ2b3RlIiwgInZhbHVlIikpCmRhdGFfNDAxJGZhY3Rvcl92b3RlIDwtIGZhY3RvcihkYXRhXzQwMSR2b3RlLCBsZXZlbHMgPSBjKDEsIDIsIDMpLCBsYWJlbHMgPSBjKCIxIiwgIjIiLCAiMyIpLCBvcmRlcmVkPVRSVUUpCmhlYWQoZGF0YV80MDEpCmZpdF8xIDwtIHN0YW5fcG9scihmYWN0b3Jfdm90ZSB+IHZhbHVlLCBkYXRhID0gZGF0YV80MDEsCiAgICAgICAgICAgICAgICAgICBwcmlvciA9IFIyKDAuMywgIm1lYW4iKSwgcmVmcmVzaCA9IDApCnByaW50KGZpdF8xLCBkaWdpdHM9MikKYGBgCgojIyMjIDYgcGVvcGxlCgpgYGB7ciB9CnBsb3R0ZWQgPC0gYygxMDEsIDMwMywgNDA5LCA0MDUsIDUwNCwgMTEyKQpzdG9yeSA8LSBjKCJQZXJmZWN0bHkgbW9ub3RvbmljIiwKICAgICAgICAgICAiT25lIGZ1enp5IGFuZCBvbmUgc2hhcnAgY3V0cG9pbnQiLAogICAgICAgICAgICJNb25vdG9uaWMgd2l0aCBvbmUgb3V0bGllciIsCiAgICAgICAgICAgIk9ubHkgMSdzIGFuZCAzJ3MiLAogICAgICAgICAgICJBbG1vc3Qgb25seSAzJ3MiLAogICAgICAgICAgICJFcnJhdGljIikKbl9wbG90dGVkIDwtIGxlbmd0aChwbG90dGVkKQpkYXRhIDwtIGFzLmxpc3QocmVwKE5BLCBuX3Bsb3R0ZWQpKQpmaXQgPC0gYXMubGlzdChyZXAoTkEsIG5fcGxvdHRlZCkpCmZvciAoaSBpbiAxOm5fcGxvdHRlZCl7CiAgI29rIDwtIGRhdGFfYWxsWywicGVyc29uIl09PXBsb3R0ZWRbaV0KICBkYXRhW1tpXV0gPC0gc3Vic2V0KGRhdGFfYWxsLCBwZXJzb24gPT0gcGxvdHRlZFtpXSwgCiAgICAgICAgICAgICAgICAgICAgICBzZWxlY3QgPSBjKCJ2b3RlIiwgImZhY3Rvcl92b3RlIiwgInZhbHVlIikpCiAgICBmaXRbW2ldXSA8LSBzdGFuX3BvbHIoZmFjdG9yX3ZvdGUgfiB2YWx1ZSwgZGF0YT1kYXRhW1tpXV0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgcHJpb3I9UjIoMC4zLCAibWVhbiIpLCByZWZyZXNoID0gMCwKICAgICAgICAgICAgICAgICAgICAgICAgICBjb3JlcyA9IDEsIG9wZW5fcHJvZ3Jlc3MgPSBGQUxTRSwKICAgICAgICAgICAgICAgICAgICAgICAgICBhZGFwdF9kZWx0YSA9IDAuOTk5OSkKfQpgYGAKCiMjIyMgR3JhcGgKCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgcGRmKHJvb3QoIlN0b3JhYmxlL2ZpZ3MiLCJzYW1wbGVkYXRhNC5wZGYiKSwgaGVpZ2h0PTUsIHdpZHRoPTgpCmBgYApgYGB7ciB9CnBhcihtZnJvdz1jKDIsMyksIG1ncD1jKDEuNSwuNSwwKSwgdGNrPS0uMDEpCmZvciAoaSBpbiAxOm5fcGxvdHRlZCl7CiAgc2ltcyA8LSBhcy5tYXRyaXgoZml0W1tpXV0pCiAgbl9jdXRwb2ludHMgPC0gMgogIGN1dHBvaW50cyA8LSByZXAoTkEsIG5fY3V0cG9pbnRzKQogIGZvciAoaV9jdXQgaW4gMTpuX2N1dHBvaW50cyl7CiAgICBjdXRwb2ludHNbaV9jdXRdIDwtIG1lZGlhbihzaW1zWyxpX2N1dCsxXS9zaW1zWywxXSkKICB9CiAgcyA8LSBtZWRpYW4oMS9zaW1zWywxXSkKICBwbG90KGRhdGFbW2ldXVssInZhbHVlIl0sIGRhdGFbW2ldXVssInZvdGUiXSwgeGxpbT1jKDAsMTAwKSwgeWxpbT1jKDEsMyksCiAgICAgICAgeGxhYj0iVmFsdWUiLCB5bGFiPSJWb3RlIiwgbWFpbj1zdG9yeVtpXSwgeWF4dD0ibiIpCiAgYXhpcyAoMiwgMToobl9jdXRwb2ludHMrMSkpCiAgdGVtcCA8LSBzZXEoMCwgMTAwLCAwLjEpCiAgcHJvYiA8LSBhcnJheShOQSwgYyhsZW5ndGgodGVtcCksIG5fY3V0cG9pbnRzKzEpKQogIGV4cGVjdGVkIDwtIHJlcChOQSwgbGVuZ3RoKHRlbXApKQogIHByb2JbLDFdIDwtIDEgLSBpbnZsb2dpdCgodGVtcC1jdXRwb2ludHNbMV0pL3MpCiAgZXhwZWN0ZWQgPC0gMSpwcm9iWywxXQogIGZvciAoaV9jdXQgaW4gMjpuX2N1dHBvaW50cyl7CiAgICBwcm9iWyxpX2N1dF0gPC0gaW52bG9naXQoKHRlbXAtY3V0cG9pbnRzW2lfY3V0LTFdKS9zKSAtCiAgICAgIGludmxvZ2l0KCh0ZW1wLWN1dHBvaW50c1tpX2N1dF0pL3MpCiAgICBleHBlY3RlZCA8LSBleHBlY3RlZCArIGlfY3V0KnByb2JbLGlfY3V0XQogIH0KICBwcm9iWyxuX2N1dHBvaW50cysxXSA8LSBpbnZsb2dpdCgodGVtcC1jdXRwb2ludHNbbl9jdXRwb2ludHNdKS9zKQogIGV4cGVjdGVkIDwtIGV4cGVjdGVkICsgKG5fY3V0cG9pbnRzKzEpKnByb2JbLG5fY3V0cG9pbnRzKzFdCiAgbGluZXMgKHRlbXAsIGV4cGVjdGVkLCBsd2Q9LjUpCiAgZm9yIChpX2N1dCBpbiAxOm5fY3V0cG9pbnRzKXsKICAgIGxpbmVzKHJlcChjdXRwb2ludHNbaV9jdXRdLDIpLCBpX2N1dCtjKDAsMSksIGx3ZD0uNSkKICB9Cn0KYGBgCmBgYHtyIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CmlmIChzYXZlZmlncykgZGV2Lm9mZigpCmBgYAoK