Better discrete choice modeling through the rank ordered logit

Jim Savage calls the MaxDiff model of discrete choice a “mathematically incorrect model of a psychologically incoherent concept”¹.

Despite this lovely dunk, and some wonderful notes explaining why MaxDiff’s not great, the model remains frequently used, most prominently in market research, as implemented in Sawtooth. Why is this?

Beyond the most obvious explanations like inertia, it’s also surprisingly hard to find a fully fleshed out explanation of alternatives online. In this post, I aim to remedy that, with sections for:

1. An introduction to MaxDiff and its limitations
2. An illustration of how these issues impact final model quality
3. Alternatives: rank-ordered logit with connected choice graphs, and rank-ordered logit with unknown middle options

Introducing MaxDiff

If you’re reading this post, you likely have some familiarity with the basics of discrete choice models, but here’s a brief refresher².

We want to build a model of how people make decisions among discrete options. This could be:

• A pollster asking which candidates or political parties each respondent favors
• A marketing firm studying preferences amongst a variety of chocolate bars
• A political scientist asking which message voters find most convincing

One reasonable model³ for this is to say that individual i making choices among j options can be inferred to have some underlying utility \mu from the available choices, and that while they usually choose their most preferred option according to their utility function, there is some degree of randomness. The basic model is then

U_{ij} = \mu_{ij} + \epsilon_{ij}

where the \mu_{ij} can have rich provenance (demographics and other individual traits, the context in which the decision is made, the other options available…) but is fixed, and there’s some \epsilon_{ij} of randomness involved. To make this model easier to estimate, we assume that \epsilon_{ij} has a Gumbel⁴ distribution.

Under this model, the probability that individual i chooses alternative j from a choice set C_i is:

P_{ij} = \frac{\exp(\mu_{ij})}{\sum_{k \in C_i} \exp(\mu_{ik})}

This is the standard multinomial logit. We could further decompose this \mu_{ij} into a linear function of explanatory variables, but this is not central to my point, so I won’t build up notation for these.

Now, suppose we have a dataset where each individual i has made a choice y_i from their choice set C_i. The likelihood of observing this dataset under the multinomial logit model is:

L = \prod_i P_{iy_i} = \prod_i \frac{\exp(\mu_{iy_i})}{\sum_{k \in C_i} \exp(\mu_{ik})}

The likelihood is the product of the choice probabilities for each individual’s observed choice y_i. We can estimate the parameters of the utility function \mu_{ij} by maximizing this likelihood function (or, more commonly, the log-likelihood). This is how the single best choice data is naturally incorporated into the likelihood function. Each observation contributes a term to the likelihood based on the probability of the chosen “best” alternative under the model.

While we can (and will) push the basic logit model of this further to include respondent and choice level covariates, multilevel components, and other improvements, let’s think about the pressures of gathering data here for a moment, since that’ll motivate the desire for something like MaxDiff.

To estimate this, we recruit participants and ask them to choose their favorite option amongst a given choice set. As a way to control the difficulty of making a choice while still gathering enough data, we can limit the size of the set (choose the best of 10 items, instead of 20), and repeat the choice task.

Getting respondents and getting them to stick through a bunch of choice tasks is hard though, and it’s only natural to wonder: can we extract more with each choice set? One option would be to ask the respondents to rank ALL the options at once, but if you have a large choice set that sounds exhausting.

What if we asked people to choose their best and worst choices each time? After all, people might not have strong preferences amongst the middling 18 chocolate bars, but the best and worst seem more memorable, and that’s only 1 more choice.

And now, what if we don’t want to entirely change up our likelihood to handle the new rich source of data? Instead, what if we just treat the worsts as the opposite of the bests, which makes a sort of sense, and simplifies the likelihood one hell of a lot:

U_{ij}^W = -U_{ij}^B = -(\mu_{ij} + \epsilon_{ij})

Under this assumption, the probability of individual i choosing alternative j as the worst is: P_{ij}^W = \frac{\exp(-\mu_{ij})}{\sum_{k \in C_i} \exp(-\mu_{ik})}

… Now we’ve stepped in it, and Jimmy is mad at us.

Psychologically Incoherent

What I’ve introduced is the core of the MaxDiff formulation of discrete choice, before bells and whistles are introduced. This has some deep problems though; let’s start with the “psychological” ones. Human decision-making is an incredibly complex, not always logical process, and we’ll always be losing significant fidelity in boiling it down into a model. Here though, I’ll focus on explaining a handful of breakdowns in the relationship between reality and the world of MaxDiff models that are particularly harmful.

First, let’s talk about symmetry. With the MaxDiff likelihood above, we’re treating the worst as equal and opposite to the best. For example, though, I don’t like Joe Biden as much as I dislike Trump⁵. I really really like the best pizza in Brookyln, but New York pizza is all New York pizza, so it only gets so bad. This might be a reasonable simplification in some cases, but it’s hard to argue that baking this into our model faithfully mirrors reality.

Also, we’re not only asking them to be symmetric, we’re sort of conjoining the best and worst choice, asking them to share the same utility scale. In other words, the factors that make an alternative more attractive for the best choice are assumed to make it equally less attractive for the worst choice. An easy example is something like health risks- “that sushi place gave me food poisoning” is very relevant to my choice of worst restaurant, but the moment I have to think about food safety, a restaurant isn’t really anywhere relevant on the “best” side of the spectrum for me. Again, you can probably think of a case where this is a fine approximation, but in an ideal world, we won’t force ourselves to weld our notions of best to our notions of worst.

Finally, the error variances being treated as the same should feel pretty strange. I’m much, much more consistent in my selection of “bests” than “worsts”- why would I spend a bunch of time deciding which opinion of 20 is the absolute worst and which is just 19th worst? People tend to be much less consistent, and frankly much less engaged, with their worst choices. Why would we bake this into our model?

Putting this all together, the simplification in estimation that comes with MaxDiff also leaks into how the model “sees” the decision maker, and it meaningfully distorts the map in a way that does not necessarily reflect the territory.

Mathematically Incorrect

Beyond calling the model psychologically incoherent, Savage also says the model is mathematically incorrect. This feels like a slightly stronger insult, and indeed it is, in the sense that the model is failing on its own terms.

How? Well, we’ve explicitly laid out a model with Gumby errors:

Wait, no sorry sorry⁶. we’ve laid out a model with Gumbel errors \epsilon_{ij}. As a reminder, the Gumbel distribution is chosen for its mathematical convenience, as it leads to a closed-form expression for the choice probabilities in the logit model. the PDF is plotted below:

Probability density function of the Gumbel distribution. Source: Wikimedia Commons, https://en.wikipedia.org/wiki/Gumbel_distribution

Here’s where the trouble starts: like we discussed above, the MaxDiff model forcibly inserts an element of symmetry into the error distributions of the choices. However, as you probably already noticed, the Gumby Gumbel distribution is not a symmetric distribution! This is fine in the case where we’re just reasoning about the best choice, but even on its own terms the MaxDiff formulation doesn’t quite make sense.

The defects are not just theoretical or cosmetic: an illustration

So I’ve shown some ways I claim that the MaxDiff model of discrete choice falls short, but how much do they really matter? As someone with an appreciation for the agnostic perspective on statistical modeling, I think it’s important to not only show that there are theoretical senses in which a model might be flawed, but further prove that these technical blemishes harm model performance on metrics we care about.

To do this, let’s generate some synthetic data, and fit both the best-choice and MaxDiff models to it. For this section, I’ll be starting from Jim Savage’s best choice model and simulated data here, and then add in the “worst choice as negative best choice” logic afterwards.

Let’s start with the synthetic data. If the exact DGP isn’t very exciting to you, feel free to skim or skip this section at first. I’ll give an overview after the block.

Data Generation for Simulation Study

library(tidyverse)

Warning: package 'tidyverse' was built under R version 4.2.3

Warning: package 'ggplot2' was built under R version 4.2.3

Warning: package 'tibble' was built under R version 4.2.3

Warning: package 'tidyr' was built under R version 4.2.3

Warning: package 'readr' was built under R version 4.2.3

Warning: package 'purrr' was built under R version 4.2.3

Warning: package 'dplyr' was built under R version 4.2.3

Warning: package 'stringr' was built under R version 4.2.3

Warning: package 'forcats' was built under R version 4.2.3

Warning: package 'lubridate' was built under R version 4.2.3

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.3     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.4.3     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(rstan)

Warning: package 'rstan' was built under R version 4.2.1

Loading required package: StanHeaders

Warning: package 'StanHeaders' was built under R version 4.2.1

rstan version 2.26.13 (Stan version 2.26.1)

For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)

Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file

Attaching package: 'rstan'

The following object is masked from 'package:tidyr':

    extract

library(ibd)

Warning: package 'ibd' was built under R version 4.2.3

library(crossdes)

Warning: package 'crossdes' was built under R version 4.2.3

Loading required package: AlgDesign
Loading required package: gtools

library(combinat)

Attaching package: 'combinat'

The following object is masked from 'package:utils':

    combn

options(mc.cores = parallel::detectCores())

set.seed(604)

# Again, note that this code is just reproducing https://khakieconomics.github.io/2018/12/27/Ranked-random-coefficients-logit.html,
# although with a more real-world feeling number of respondents and
# and total choices. If this is taking uncomfortably long, feel free to scale
# down I, the results should reproduce just fine outside irresponsible I.

# Number of individuals
I <- 800
# Number of tasks per individual
Tasks <- 10
# Number of choices per task
J <- 5
# Dimension of covariate matrix
P <- 5
# Dimension of demographic matrix
P2 <- 6

# demographic matrix
W <- matrix(rnorm(I*P2), I, P2)
# Loading matrix
Gamma <- matrix(rnorm(P*P2), P, P2)

# Show W * t(Gamma) to make sure it looks right
W %*% t(Gamma)

               [,1]          [,2]          [,3]         [,4]         [,5]
  [1,] -1.147325812   1.453389560 -0.1201172125  0.964840281 -4.537077496
  [2,]  0.274237431   2.755783152 -3.0698355631  1.205484626 -0.334220382
  [3,] -1.776487023  -3.676627003 -4.1725965643  2.232580311  0.030518687
  [4,]  1.303242055   2.796902290 -1.0818975319  2.213676892  0.460143743
  [5,] -0.038077610   1.013494581  1.5886972343 -0.357147508  0.735694540
  [6,] -0.609004672  -2.399734277 -2.8631005810 -1.385657832  0.701776800
  [7,] -0.003016024  -2.605408894 -0.0723969226 -0.315353540  2.563850493
  [8,]  2.091206040  -0.275019611  0.7615789884  0.979190201  2.017939901
  [9,] -1.601434736  -1.203864883 -1.8161995764  0.075095339 -0.086070281
 [10,]  2.181393983   0.010819374  3.0136121632  0.444287730  3.280880660
 [11,] -0.923930456   5.505541925  0.2091929622 -0.312569250 -3.015452664
 [12,]  2.076792020  -4.464111386 -0.7351549450 -1.968268039  4.015238464
 [13,] -3.104754729  -0.703419923 -2.9470505419  0.066057565 -3.054629868
 [14,] -0.880319086  -3.420345971 -0.0277780987 -2.667776798  1.830987262
 [15,]  3.226864181   0.529723362  1.5324612847 -1.871606991  2.964084007
 [16,] -0.930974513  -3.108637587 -1.3524866917 -2.109222757  2.365901280
 [17,] -1.310019841  -4.169604411 -2.2627166822 -1.266567496  1.043406974
 [18,]  3.161574540  -2.611419524  3.0503506088  0.227149338  5.338110678
 [19,] -4.312627339   2.514695771 -3.2054050106  1.419674971 -5.850313526
 [20,] -0.369833548  -0.850701972  0.0437976453 -0.814807451 -0.217777666
 [21,] -0.093495250  -2.056789230 -1.9104855848 -0.660730465  1.329483831
 [22,] -0.067175662   1.879783406  0.4375834364  1.740624948 -2.363759883
 [23,] -0.856613421  -3.888184893  0.6467183551 -1.229485027  1.794434607
 [24,]  0.687448612  -2.803472954  1.6679059947 -0.254895844  1.529721953
 [25,]  1.720909217  -6.504963652  2.0340556967 -0.116379708  4.536006412
 [26,]  0.208506170   0.882608635 -0.4012337183 -0.465793586 -0.879646390
 [27,]  0.112184349   1.525029004 -1.7386757284 -0.331554096 -0.164149516
 [28,] -2.119012895   2.134066676 -0.3907224457  0.074333553 -3.719055698
 [29,]  2.480228434   2.103264918  2.5432286975 -1.296866193  1.938367170
 [30,]  2.138337683  -0.926221470  2.0946608957  0.926451734  3.379799523
 [31,] -0.914561807  -6.540433106 -0.4179541821 -1.086444475  3.049195665
 [32,] -1.062292815  -1.300877435 -0.4361758507 -0.799892721 -0.335557689
 [33,]  0.245720595   3.073001625  0.6991710796 -0.074686140 -2.225227033
 [34,] -0.377272607  -0.353870328 -1.3458442750 -0.910899967 -0.565737827
 [35,] -0.295052509   0.091813224  2.0721980064  0.882913806  0.498928760
 [36,] -2.521616146   0.977427103 -4.3256955047  1.887656083 -4.298542059
 [37,]  0.080044948  -3.371503258 -0.2716704826  1.157519771  2.152985108
 [38,] -0.008685601  -0.656543690 -0.3087439318  0.093581253 -0.702426038
 [39,] -0.452865455   2.817498436  0.8123390925 -1.872157949 -2.724494438
 [40,] -2.485016731   4.433176832 -2.3819198326 -2.055012388 -3.777695095
 [41,]  0.919314701   1.311984501 -0.0900807513 -2.035894542  0.948855357
 [42,]  1.187100142   1.360097874  0.6672974682  0.964509720  1.670374710
 [43,] -1.581338997   1.866155832 -0.3658008112  0.905249003 -2.945159079
 [44,] -1.889823182  -3.081847141  1.0495166855 -0.385461308 -0.853659352
 [45,]  1.436584668  -0.518680614  0.5533374897 -2.208769780 -1.039659093
 [46,] -1.491045520   2.237414633 -1.3836291477 -4.248664708 -1.247598550
 [47,]  1.341838094   1.264542976 -0.4475810722 -2.563594967  1.510455298
 [48,]  2.743537922  -5.020213997  0.8264976930 -0.344415550  5.410768104
 [49,]  2.603513418   3.346058668  2.6413899936  3.800534752  0.640711970
 [50,] -1.139672402   7.052317591 -0.1562487781 -0.163534252 -7.507161898
 [51,] -1.893776574   2.040924415  0.3700712218  1.168897436 -3.361268354
 [52,] -2.766244964  -3.350457840 -2.2300027111  0.466177598 -1.176316863
 [53,]  1.051265095  -0.367123446  1.0517269291  0.408515674  1.643631211
 [54,] -1.236930061  -1.185028345 -1.3981317392  0.324433518 -0.908758321
 [55,]  0.340330445   1.845788495  1.8664202308  0.252825786 -1.005874480
 [56,]  3.918194976  -4.413687519  1.4843814220 -0.458038880  6.665507986
 [57,]  1.261571660   5.578593399  2.1264783330 -1.211710550 -2.051998156
 [58,]  0.132447017  -0.183630094 -1.3015783488 -3.174628410  0.761193274
 [59,]  0.177767767   0.642973124  0.5653449537  2.824666462 -0.775694646
 [60,] -3.016442920  -1.595111319 -2.9277032233 -1.767701996 -1.510719034
 [61,] -0.480377639   0.060127459 -0.0077221630  0.795443333 -1.187731088
 [62,] -3.543415601   0.429381240 -2.6984588659 -1.471647171 -4.774057565
 [63,]  0.286764561  -0.579039329 -0.4605559673  0.854811811  0.751781396
 [64,] -2.034091897  -0.263222388 -1.7213749996  1.194130069 -2.888655258
 [65,]  1.512810608   2.608794043 -0.0938323589  3.544357756 -0.277532923
 [66,]  3.255482897  -4.435563394  3.6355028985  3.623299587  6.258966684
 [67,] -0.412404419  -1.215755986  1.4495470425  0.076860781 -0.999827874
 [68,]  2.367336905  -0.781368179  3.1853606621  0.775734078  1.193184720
 [69,] -0.902277060  -2.145380764 -0.1414340913  0.229382106  0.120942905
 [70,]  0.584096490  -0.238561426 -1.8005325973 -0.193155113  0.556455753
 [71,]  2.283679546  -3.287338668  0.8685959457 -3.182096474  3.882607935
 [72,]  4.626043086  -3.939590808  4.9731979922 -0.596820342  6.225318515
 [73,]  1.951199285   1.150419025  1.5934541582  0.059301868  2.308793191
 [74,]  1.496241974  -6.335195985  1.2603894395  2.571418344  3.929775887
 [75,] -2.538895792  -2.766574835 -2.6999004387 -2.014100295 -0.403374526
 [76,]  1.778949680   2.650396428  2.0791506946 -0.037130488  0.534264533
 [77,] -1.212100728   2.442953752  0.1950404092 -0.610212638 -2.689368809
 [78,]  0.887260911  -6.483892589 -0.3126799591  0.856856722  5.702372152
 [79,] -0.627714623   2.483061948  2.6897637129  0.003376006 -2.631401234
 [80,]  0.402922728   0.441148139  0.6169147116  0.421889156  0.577734315
 [81,] -1.419340454   0.025800748 -0.0351987058 -1.712028649 -0.367138194
 [82,] -0.818013748  -0.418637206  0.0372427652  0.851869216 -1.507908724
 [83,] -2.220465012   3.144092238 -0.6640153584  0.748190825 -5.330762918
 [84,]  2.365394529  -0.192615206  2.4255379176  0.402529096  3.131663053
 [85,]  0.352101145  -1.892590125  1.5660426646  0.069608994  2.247204152
 [86,] -3.826769157  -0.266559135 -1.5387074416 -2.220052962 -3.165995840
 [87,] -0.572952327   1.218360954 -0.3820010251  1.427246541 -0.835182239
 [88,] -1.375655142   3.362327590 -1.0818254409 -0.046499489 -3.869823850
 [89,]  0.697157410   3.139312336  1.5479172301  1.177933414 -2.870252935
 [90,]  0.349982893   4.853684518  2.2167216829 -0.751651949 -3.435385031
 [91,]  2.689693021  -2.154542802  3.9248033242 -0.632343115  5.006995771
 [92,]  2.147461559  -3.566419990  0.2148003851  0.554312904  4.093823466
 [93,]  1.434951093  -1.436192351  3.5370133246 -0.821415867  0.049581963
 [94,]  1.776298329  -5.767818211  0.2352658101  3.618976060  2.999211025
 [95,] -0.223138997  -2.122979775 -0.7034167588 -0.924787202  1.518849012
 [96,]  1.023487867   4.395375208 -0.0827556764  2.367699021 -2.796291231
 [97,] -3.063696373  -1.642236719 -3.7699827007 -0.474634240 -3.936073520
 [98,] -0.629739909  -0.776520856  0.7733825328 -2.600767739  0.289085852
 [99,]  1.997778954  -6.260621815  2.1883804789 -0.807008086  4.849920887
[100,]  3.537479830   3.467595377  4.3234933259 -1.769656581 -0.132934189
[101,] -0.836774957  -0.349898988 -1.1197247128  0.662003198 -2.374277301
[102,] -1.714479293  -0.606878227 -2.0338281892  1.460501154 -1.049031311
[103,]  3.469213483  -3.382349105  4.1260182008 -1.662649743  6.293060752
[104,]  3.201526813   6.244561529  5.6044489592  1.553683445 -0.433072498
[105,] -2.484505799  -2.132311127 -4.6687197790  1.251026152 -1.060343818
[106,]  2.959983547  -0.381444310  2.9977317458  0.249857660  2.808447698
[107,]  0.666856614   0.444735394 -0.7317306334 -0.737839466  0.313371564
[108,] -1.384671399   0.982298534 -3.2105810866 -2.952758507 -1.680539934
[109,] -0.222921826  -0.988743910  1.6294818082 -1.481295193  0.173106656
[110,]  1.590595327   5.404892300 -0.7737593719  1.153015865 -1.599618363
[111,]  2.818008297   0.300213470  1.9323230304 -0.505861613  1.526842614
[112,] -1.150525193  -0.739456189 -1.2902837050  0.825635918 -1.473932172
[113,] -3.160246512  -1.609875092 -1.8634933032  0.507429702 -1.377127118
[114,]  1.853285765   0.838386270  0.4442075115  3.588161468  2.129641137
[115,] -0.038168723   4.152736587  0.8552937100  2.358162578 -2.091705045
[116,] -3.041824445   1.222596263 -2.1187971634 -0.632459855 -2.940498959
[117,]  0.005854236  -2.220440695  1.8771692699  1.855087160  1.698327886
[118,] -1.745411590   1.756340914  0.8794652293  0.049838402 -2.658977182
[119,]  2.558215172  -0.057773934  1.6972909350  1.980287155  1.595839232
[120,] -1.346182514  -4.261041304 -0.7180215449  0.124573680  1.589423937
[121,] -1.167097720   1.327505281 -0.6871473661  1.294948108 -2.497022633
[122,] -1.580601817  -3.915179576 -0.5468114295  1.327796569 -0.395895937
[123,]  3.543305534  -2.749472513  2.7333434678  0.913761437  5.013309927
[124,] -4.862137216   1.821961104 -6.6153447188 -2.772410920 -4.926524351
[125,] -2.565802516   0.048958460 -1.0522613918 -0.560930688 -2.272633542
[126,] -0.406992020   1.438129749  0.5620736673  1.970858739 -2.349174500
[127,] -0.728408129  -4.606261717 -1.9114059999  0.237351971  2.245823980
[128,]  0.954292400   0.240677679  2.0974612476  0.850151298 -1.883533392
[129,]  0.747923124   0.095019392 -1.2441459810 -0.510758282  1.608079453
[130,] -0.914430334   2.124500772  0.4885946238  2.999426585 -0.429998181
[131,] -0.359819540  -0.253552587 -0.1132966961 -0.200815617 -1.648253615
[132,] -0.139328268   7.034489159  0.6250800426 -0.170535689 -3.810217970
[133,]  0.953522385   1.884010369  1.1547927881 -0.312382764  1.379458936
[134,] -1.857727648  -1.146619917 -2.7326924813 -0.857607617 -0.595583835
[135,] -2.476550853   5.792733745 -0.3408842392  1.078955340 -5.819511378
[136,]  0.541670970  -3.180398748 -0.1610077968  2.305518553  2.575943006
[137,]  0.309640572   9.714652930  0.0938091071 -1.479894037 -3.092546972
[138,] -1.397122976   0.136784980 -3.3757715659 -1.426894182 -1.653748468
[139,] -2.612448104  -3.585558658 -2.6868253571 -0.171930483 -1.114902143
[140,]  0.352546035  -4.240933874 -1.0407189143  0.208588905  3.371469349
[141,] -3.668231027  -1.069134966 -3.8915683166 -1.000472226 -2.664830096
[142,]  2.553308174 -11.112539431  0.1167527872 -0.856026250  8.402730480
[143,] -0.424918290   1.624886199 -1.9578544984 -2.600295532 -0.600916472
[144,] -1.287314256   1.276821566 -0.9678671244 -1.268447260 -0.969383679
[145,] -0.632165862  -2.997690519  0.6401426790 -0.885560634  1.909204796
[146,]  0.532830858   1.332600922  2.0603421220 -0.416613225 -0.901738534
[147,]  1.540560940   3.561378335  1.8650690538 -0.818773162 -0.804894245
[148,]  0.629580556   2.245842790  2.5548795267  2.812633031 -1.467932558
[149,]  0.689517253  -4.412227181  1.1885964108 -1.431090347  3.518173141
[150,]  1.469975030  -3.505353628  2.2210507752 -1.623588521  4.614175476
[151,] -1.879405439   2.492590933 -0.8209443349  1.358453023 -1.932861038
[152,]  1.284802681  -1.603351763  2.5509952267  0.103186345  1.590995771
[153,] -0.518357312   1.520496894 -0.1049530572 -1.360796854 -1.002540873
[154,] -2.035583786   1.540903107 -1.6809522726 -0.798315062 -1.372671547
[155,]  1.365243210   1.856570645  3.0469864077  0.535801057  0.767015498
[156,]  1.700472147  -0.816657576  2.2121371353  0.271139362  1.990247881
[157,] -0.423109570   0.971688184  0.7216906450 -0.083891209  0.118230267
[158,] -0.841208862  -0.442723728  0.4199048893 -1.778977801 -0.867300281
[159,]  3.468085784   0.317782620  2.2703820574  1.988241109  1.442909750
[160,] -1.011212013  -2.540982305 -1.2813905508  0.781692721  0.124039233
[161,]  0.392894704  -5.317134588  0.7811303497  1.429310491  2.490071184
[162,] -1.931843597   4.907334086 -2.1455729645  2.647480778 -6.466075552
[163,] -3.781403832   4.017454293 -4.3496406091 -1.628341787 -6.628727306
[164,]  0.980492061   1.558050655 -1.1222253855 -1.024031550 -0.238917623
[165,] -3.620149394  -3.180970580 -4.6171890206 -1.784218351 -1.582808941
[166,] -1.789439979   0.421005361 -2.3126583691 -2.545659252 -1.124149187
[167,] -3.037035863   1.482902779 -4.0808936706  0.229985434 -3.715256993
[168,] -0.451051708  -3.739637759 -2.1926029239  0.569611062  1.474048658
[169,]  2.558105304   3.578045141  3.4994877056 -2.598877095  0.521487280
[170,] -2.668865263  -2.937307995 -3.6000880426  1.006330533 -2.308283910
[171,] -0.629253656   2.316638291  0.1959717326 -0.843487591 -0.458876760
[172,] -2.526830485   5.580101549 -2.1835854767  2.864041538 -6.282063257
[173,] -2.896383094   1.613015994 -2.1882028789  0.953678064 -2.812582555
[174,] -2.126136747  -1.408927322 -1.4901156366 -0.108056146 -0.810012782
[175,]  0.032191073   0.286540791  3.5306208729  0.461497680 -1.039204827
[176,]  2.879985756   0.305246777  1.5081647033 -0.688096990  3.467366126
[177,]  1.742609063  -0.591923248  0.5355973194 -2.221030509  1.982229934
[178,]  0.867383883  -1.223075770 -0.7517734651 -1.367410397  1.913453886
[179,] -0.104577213  -1.608611576  0.5709293862  0.332918428  2.004652769
[180,]  0.818823197  -3.553940675  2.3051189573  0.400908785  2.388338243
[181,] -1.963965677   0.810941534 -0.6998286687 -0.130454708 -2.936106038
[182,]  0.866753563   1.103379573  1.4701392916 -0.384914150 -0.719300450
[183,]  1.075088091   1.542711228  0.9604435703 -0.447268544  1.376065136
[184,] -2.717763276   3.616257801 -2.8156324972 -2.434903485 -5.123506644
[185,]  1.669532671   2.460479907  1.6832929387 -0.923080548  0.999847914
[186,] -4.648863959   2.461572524 -2.8060141910 -2.058874784 -6.390285650
[187,]  0.746116196  -5.709606201 -0.4536416437  1.684924215  2.310522161
[188,] -1.688720688   0.098756922 -1.7480798713 -1.557039792 -1.497607268
[189,] -1.363626484  -0.492397127 -1.8774761491  2.936029454 -1.015145119
[190,] -1.189492345  -3.536179736 -1.6973115124 -2.850060289 -0.150286136
[191,] -0.920743392   1.702088131 -1.2613004554 -0.128682949 -1.797273963
[192,]  3.029813673  -1.604138219  0.6448555084  0.614426166  3.527531324
[193,] -0.534149467  -6.264726120 -0.3782466940 -0.723594171  2.101432577
[194,] -3.824601917   2.815285338 -3.8322990385 -1.218507638 -5.082096098
[195,]  1.739799546  -2.621197389  1.2704575938 -0.771554377  3.023381451
[196,] -0.139306618   0.021946800 -1.0953285386 -0.273040643 -0.615705881
[197,]  1.242138598  -0.827290958  2.2803561329  2.090801097  0.281611843
[198,]  2.547365427  -0.718370934  1.7379140599 -1.554420533  3.651011434
[199,]  1.690492059   2.615566494 -0.3004369557  2.216552860 -0.028995586
[200,] -1.476849894  -0.210325358 -0.5894008516  0.917000417 -1.319248572
[201,] -0.449010798   0.731538415 -3.2950548238  0.585151802  0.499530338
[202,]  1.996812641  -2.538206587  1.4092494807 -0.851289275  2.972899004
[203,] -2.273057123   0.443327068 -1.5600480614 -1.187380822 -1.187891044
[204,] -3.039237390   1.231504931 -2.4029885110 -2.747038189 -4.058419031
[205,]  0.816166718   0.063374317  1.5774084078  4.547430572 -0.822935052
[206,] -0.814368598  -3.236086029  0.8560657174  1.310962088  0.753733122
[207,] -1.310726016   2.899045668 -1.1873692087  2.292840427 -4.459857007
[208,]  1.897888537  -3.360258420  0.9642859479  0.174503671  3.264650382
[209,] -1.377524064   5.147961000 -0.9819411254 -1.225044732 -4.161151863
[210,]  0.327902215  -0.698974873  0.7186816108 -0.694784244  0.234130780
[211,] -1.448902136   6.259048533 -1.2164213904 -1.411471142 -5.400895580
[212,] -0.522003051  -1.739029394 -0.8670067012 -2.770375227 -0.561036385
[213,]  2.673411988  -1.016495853 -0.6747520188  0.591249191  2.784038960
[214,] -1.215892977  -0.415808045 -1.0895045534  0.820895793 -2.250666618
[215,] -0.754155402   1.813781340 -1.0290878069 -0.609781967 -1.834208817
[216,] -0.807334091   7.068133310 -0.0821869415 -0.073573285 -2.995571770
[217,]  0.812653201   1.779177053 -0.2594524909  1.323654156 -1.095825872
[218,] -0.896751958  -0.361832540 -1.9362862156  1.795237499 -2.103083333
[219,]  0.209757854   3.583665799  0.2326097286 -1.566542797 -1.950016750
[220,] -0.607822241  -3.538951600 -2.7917321052 -0.442502962  0.283975904
[221,]  0.052146486   0.346397912  0.3416016254  1.463618957 -0.675893045
[222,] -5.310786249   5.767981616 -3.5079415103  0.916968336 -8.040401468
[223,] -0.691953186  -8.173689661 -0.3346364606 -0.215522285  3.140273580
[224,]  1.991069616  -0.596374018  0.5708604947  2.699472595  2.309430888
[225,]  0.078614759  -2.648703771  0.5062579158 -0.529760198  2.363780176
[226,]  2.175891410  -2.126548526  3.1736763034  1.294763122  2.463010180
[227,]  1.812951905   6.012398846  2.9967276141 -1.067353000 -3.450174183
[228,]  0.453049361  -0.232522850  1.4115287786 -1.736627297  0.507827544
[229,] -0.789007351   6.816253416  0.4974357963  1.120724203 -4.784708633
[230,] -1.575380204   1.172590311 -1.8607763401 -0.481842623 -3.691723827
[231,]  4.536526172  -3.156827344  1.9730953283  1.536107112  7.224333166
[232,]  0.304963349   0.036860067  0.8908635300 -1.785552557 -0.084320804
[233,]  0.551523434  -2.776311271  0.9790498414 -4.106160804  1.939161011
[234,]  0.703926169  -2.633192133  2.1855389175  0.668578013  2.884839927
[235,]  0.538887098  -2.163024200 -1.4539048566 -0.202130035  2.363363480
[236,]  1.367698124  -8.225190890  1.9844771734 -2.731895856  5.789233515
[237,]  2.824749194   4.321248513  1.8156925135 -0.427485763 -0.400294495
[238,]  1.393585533   1.480917362  0.7577597501 -1.357877401  1.070671309
[239,]  0.163552603  -0.578976550 -0.4883511487 -0.124428028  1.394408287
[240,] -0.920475884  -2.461520224 -1.0119835104  0.572654228  0.774039116
[241,]  1.982356203  -2.104508622  1.9643822426  1.842204669  3.949965471
[242,]  0.896462357  -1.484021354  0.5511371570 -1.295333668  2.809166041
[243,] -0.511007131  -0.747369394 -1.4982529998 -0.720837309 -0.049399521
[244,]  0.917626129  -3.291233076  1.3945606344  0.921575372  1.561545491
[245,]  0.839547974   1.650484012  0.1382287949  0.052416342 -0.932493792
[246,] -0.885494723   2.243508299 -0.5336227929  0.441883190 -1.775383435
[247,] -0.561973385   3.022758903 -0.7192179181 -0.457258703 -3.423247083
[248,] -1.305460710   5.253572824 -0.7021705409 -0.595639530 -3.786391834
[249,]  1.780132470   2.952712483  2.4939830867  1.475215234 -1.109003757
[250,] -1.356617662  -3.798825450 -2.3723579240 -2.649162422  2.065106628
[251,]  0.620675684  -1.259884049 -0.4700035428 -0.532646813  1.757537363
[252,] -0.058443036   3.106592500  2.7416651365  0.732194809  0.394657339
[253,] -1.159981538  -2.798949470 -0.8226539238 -2.607236880  0.988752042
[254,]  0.422641461  -0.864548204 -0.5650287922  1.409186249  0.432565093
[255,] -1.846547856   2.987349625 -2.4850219076  0.048490928 -2.421457945
[256,] -0.686255788   9.036875671  1.0756765000 -2.743353993 -4.496535647
[257,]  3.477634078   4.277376064  4.0902262869  0.433673400  1.815721651
[258,]  0.366882229  -6.055565213 -1.8476864156 -1.358202013  3.696759737
[259,] -1.692373702   5.386758317 -2.1524104913  1.672686577 -4.697762049
[260,]  1.943820084  -3.332214806  2.8575517713  1.140286571  3.554252967
[261,]  2.361391236  -2.146309725  2.4414752533  0.620583802  2.471858473
[262,]  2.451122721  -2.642804980 -0.3167921963  0.555782212  2.910926101
[263,] -2.562877130  -1.957868683 -3.2157912395  0.212888819 -0.522834488
[264,] -0.314433381   2.246267515 -0.7856006930 -0.394418069 -0.526038409
[265,]  0.288729938   0.222169307  1.6599138499 -1.376314073  0.862982408
[266,]  1.136783569   0.656235064  2.4824861239  1.964967849 -0.784760141
[267,] -2.824170801  -3.235570874 -2.1175974093 -0.094849037 -0.964841934
[268,] -3.277865848   1.760657537 -2.4134164484 -0.354930577 -5.195509392
[269,] -0.197985224   3.168330103 -0.1155721572  0.417135719 -1.760878324
[270,] -0.084495207  -0.421624898  1.3572956658 -1.773251038  1.509732778
[271,]  0.491427706  -3.095326065 -0.2069028117  0.995565149  2.247850372
[272,]  3.242559098   1.414395475  5.1585054446  3.069712127  2.476792126
[273,]  0.789003817   4.118414844 -0.2748804343  1.827961816 -4.126480903
[274,] -1.564744016  -0.701007429 -1.4063802825  0.597984003 -1.728985266
[275,] -0.746608048   1.346641439 -1.1332276864 -1.950216250  1.409067739
[276,] -0.250864970   2.249833478  0.2146301686 -1.452409000 -1.422875383
[277,]  0.999964906   1.542858584 -0.6411493364  1.202993774  0.896258261
[278,] -1.050198925  -3.792564575 -1.8759986065 -1.160351914  1.109446525
[279,] -1.060945275  -1.031792292 -0.3823303054  2.158079796 -1.483892572
[280,] -1.378020472   0.368215094 -3.2661952196 -1.638018888 -2.688433257
[281,]  0.153248293  -4.850542821 -0.5899299915 -0.725820418  2.630699167
[282,] -0.536571645   1.582055167  0.3709058639 -0.972421650 -1.913293489
[283,] -1.820060812   4.659002929 -1.1978200396  0.895666673 -4.961505430
[284,] -0.427099207   0.181793032 -1.7785365421  0.819952456 -1.050559958
[285,] -0.260361517  -1.423161826 -0.4064427527  2.121456521  0.169050890
[286,] -0.003748518   3.312654742 -1.1604807667 -0.417769596 -3.149985094
[287,]  2.812989352   1.619439565  4.0393513777 -0.788674282  2.080161055
[288,] -1.780721380   3.779178413 -2.1598868414 -0.536542168 -3.864738722
[289,] -2.735446228  -3.074115281 -1.5958466348 -0.747521239 -0.818903078
[290,]  0.337393852   1.543932843  0.6416575157  3.156805570 -1.568128255
[291,]  1.795486252   1.628566198 -0.4484571547 -3.143091163  2.399445758
[292,] -0.664075645  -2.415838300 -1.8822646136 -0.938905604 -0.272235763
[293,] -1.644891708   1.788543611 -0.1099643461 -1.799749507 -1.641883587
[294,] -1.835911684  -1.925047206 -3.4328655375 -2.462787270  0.475285171
[295,] -0.990776279  -1.765580053 -3.2955435297 -1.626794182  1.036678577
[296,] -0.245271457  -4.831363195 -1.5226459261 -0.043770371  3.666030816
[297,]  0.569967287   1.737752394  0.1355817010  1.510964361 -1.187605275
[298,] -0.299584538   3.360761771  1.1210215973 -1.395869183 -1.784654725
[299,]  0.477194454   2.443255325  0.5395754328  2.862167343 -0.420041862
[300,] -0.824156639  -4.132407143 -0.0219441275 -0.010716763  1.661706649
[301,] -1.260393289  -3.482049069 -0.8700705113 -0.682520467  2.323427297
[302,] -0.201296664  -1.116414363  1.7560019280 -1.767778904  0.200013786
[303,] -2.820823397   1.383841289 -1.5204703754 -3.530233489 -3.094894822
[304,]  2.219972811   3.649120217  2.9749625320 -0.055624671  0.307163272
[305,]  0.216853882  -0.665801975 -0.6269717748 -1.450179831  1.564789726
[306,] -4.800468856   5.478749919 -2.7897751169 -3.920517503 -6.904228425
[307,] -2.515283377  -2.592977459 -0.4688634456  0.244979591 -2.202451171
[308,]  1.197537457   2.379341310 -0.4595128053  2.286482608 -0.413465784
[309,] -2.012585043   2.838666030 -1.0447548720  0.552168827 -4.086716182
[310,]  2.358000477  -4.054649686  0.2669417685  1.543611910  3.295832475
[311,] -4.047714357  -4.408506285 -4.0891241465 -0.377377080 -0.866602945
[312,]  2.517577875   1.060988622  0.9718986985  1.153093903  0.714208479
[313,] -1.387150611  -0.871804702 -3.0397679104  2.696390974 -0.970477644
[314,] -0.631537950  -0.706169481 -2.6871867133 -1.471914053 -0.771659698
[315,]  2.297999881  -0.985575526  2.2902151774  1.684186562  2.943480911
[316,] -0.801397750  -5.527458296 -0.8356829415 -1.002822106  2.359299639
[317,]  3.081757791  -2.597842048  2.5077062967 -0.646799136  5.309396624
[318,]  1.526928857   1.827403311  1.5081841049 -0.087289039  0.497086648
[319,]  1.746192869  -3.156304504  2.1630106889 -2.775199500  2.833922316
[320,] -0.640844853  -1.707424612  2.6734105992 -1.223312014  1.593243152
[321,]  0.685432307  -2.761920328 -1.1438698229 -2.287135191  2.855302436
[322,] -3.580130836  -2.040506000 -5.2081651059 -2.635526086 -2.317890390
[323,] -0.008273563   2.950932294 -2.0929988845  0.024072140 -1.022442235
[324,] -2.341921052  -2.925617168 -0.8502426989  0.419269148 -2.282209548
[325,] -2.553772574  -6.352513160 -2.8675647885 -1.605145143  1.260437929
[326,] -2.058585249  -0.556104637 -1.7641074708  0.697550857 -0.925183763
[327,] -1.259480118   1.250441985 -2.8342113507  0.342850416 -0.568671546
[328,] -2.356427708   2.764884150  0.3517233016 -1.452263995 -5.334750436
[329,]  4.496090997   3.498766405  5.0683134336 -0.273712834  1.821998093
[330,]  1.035975859   1.903736935  0.0007916891  0.204415935  0.538509571
[331,] -2.451373154   5.259521405  0.2135237829 -1.624064203 -4.354078979
[332,]  0.043891910  -2.361961609 -0.7901052975  1.204027900  0.248944407
[333,] -2.532197966  -2.950513593 -1.8170632831  0.452354600  0.103487354
[334,] -1.690819506   2.680252951 -2.1882567675 -0.958951785 -4.087785035
[335,] -0.745221727   0.172301580 -0.7114413375 -1.295037090 -1.448698679
[336,] -0.828644750  -0.010770368  0.3800549781 -1.763020531 -0.307968315
[337,] -1.832777604   2.836506111 -1.4246544753 -1.184711778 -3.053256210
[338,]  1.519566763  -1.713892455  1.2092995136  2.039331146  2.340921972
[339,]  1.027249569   0.965102261  0.8566585044  0.424076590  3.098111147
[340,]  1.278320836   0.002497314  2.0690965640  0.766337606  1.453103988
[341,] -1.513215123  -5.131717998 -0.3821831180 -0.180453070  1.660833921
[342,]  0.359356037  -0.438468435  0.0917964465  0.163040652  0.562550194
[343,]  0.644615494   4.137694671  1.2838575261  0.147163915 -1.598879814
[344,]  0.213643335  -2.015199933  0.5159102127  1.464150195  1.377162392
[345,] -2.432923548   0.480195587 -1.8595872864  3.293856484 -3.005345228
[346,] -0.332333899  -0.482991292 -1.3185459178 -0.311970278 -1.687868685
[347,]  0.026409277   2.985440517  0.0160709439 -1.032682466 -0.632932543
[348,]  1.357084452   7.043885938  3.9079073029  1.740237359 -3.098172669
[349,]  0.456870932  -1.729319444 -0.4657653266 -0.414305785  0.725963334
[350,] -0.137666418  -2.158134292 -0.0741607676  0.144782912  2.753867844
[351,] -0.420192610  -7.374829689 -3.6695585011  0.860330078  4.311390698
[352,] -2.394323308  -5.389673624 -2.9626618241 -2.736825312  0.440184037
[353,] -0.068616531   3.292823390 -0.3383370646 -2.814942876 -1.256326086
[354,] -3.309667902   7.894702862 -1.7805857253 -0.084267066 -7.122762320
[355,]  1.314323070   3.444542770  1.1813573183  2.505931081 -0.398728562
[356,]  0.722143552  -4.790203071  0.4450831850  0.056212071  3.333812400
[357,] -0.771819515  -2.380335732 -1.2786225831 -2.341698298  0.332323311
[358,] -1.498312176  -5.273867306 -2.1380934069 -2.084250623  0.378591654
[359,]  0.299824084   1.449786014 -1.6976929883  3.306379370 -0.399401363
[360,]  2.329647757   1.038322514  1.6267597690 -0.684361602  1.293079293
[361,]  0.471356448   2.148515598  1.1252024542  0.546380605 -2.051023094
[362,] -1.250048799   0.678557307 -0.4973077146  1.576890530 -2.236894892
[363,] -1.514186757  -0.182337139 -0.7278413142 -3.240579937 -0.784419498
[364,] -2.996545580  -1.083828349 -2.5488364158  0.422588728 -1.995390088
[365,] -1.940556019   0.902307028 -1.3976860485 -1.367354015 -2.880138172
[366,] -2.791523173   4.076062975 -1.4399774454  1.573466276 -5.392101755
[367,] -2.408091034   1.299720336 -2.2852973584  1.881722586 -3.269797226
[368,]  3.820875226  -0.761538723  4.7306349561  0.674699068  3.557105530
[369,] -2.798760966   3.528654380 -3.0476683086  0.001778689 -4.604119147
[370,]  2.817791607   5.889432340  2.0295493274  0.719133960 -0.273473710
[371,] -0.126572242  -4.522054257  0.0357855299  0.008632871  2.232267623
[372,]  5.995531560  -2.265335206  4.8308268359 -0.629143304  7.143317426
[373,]  2.156928700   1.465282866  3.2409417105 -0.465831267  1.245576284
[374,]  1.030409310   2.948730452  1.9366168633  1.961540885  1.033447316
[375,]  0.432928558  -3.779201779 -0.7141907566 -0.221686496  1.742527060
[376,]  0.916526972   7.242238756  0.5544108669  0.105733108 -3.838443828
[377,] -3.487569245  -2.474051906 -2.4946733400  0.057406443 -2.208738315
[378,] -0.893835500   0.446854119 -0.9283223531  0.474138154 -1.091448085
[379,] -2.150241296   0.185050472 -0.6913233328  0.036837967 -2.644252901
[380,] -0.237253795   1.194272970 -0.4908948597  0.404721306 -1.693469179
[381,] -0.869537906   3.319607423 -2.2584913405  0.723497507 -2.060267399
[382,] -0.177616238  -2.339784868 -0.5337778997 -0.777606057  2.205508032
[383,] -0.565171052   0.666878155 -0.0340599129 -0.331408156 -1.545994058
[384,] -3.300993520  -0.380403127 -4.9918180665 -2.525357465 -0.335073758
[385,]  0.399802123  -1.475915971  0.0028480702 -4.508386865  2.821703958
[386,]  0.129954113   1.731369756 -0.0361766382  0.185206127 -1.513388943
[387,] -0.843443516  -0.637235310  0.5209086907  0.169651516 -0.741539078
[388,]  3.426952024  -4.500311904  1.0603295279  1.529183634  6.581931024
[389,]  3.565194231  -1.542060007  2.9750779269 -0.739484119  2.510948861
[390,] -3.280779964   4.062607431 -0.9053832876 -2.259670290 -3.966705459
[391,]  2.010960586   2.509440292  2.4796209672 -0.758585114  0.501385509
[392,] -0.346670385   0.293827742  0.8869619617  0.571785447  0.025925780
[393,] -0.176510929  -1.123325371  0.1672903444  0.387087930 -0.364143040
[394,] -0.836094906  -1.822918783 -0.1372034019 -1.030647288 -1.067025193
[395,] -2.959071753   6.044727186 -1.4066800674 -0.669100636 -5.948108958
[396,] -1.353431160  -2.943140168 -2.6669603898 -0.688079948 -1.448911292
[397,]  0.117004261   0.975015560  0.9054251362 -2.539219140  0.265730487
[398,]  0.356975915   0.029392264  0.3241238195 -0.785106534 -0.571822651
[399,] -2.023034288   3.893035331  1.9660923016 -1.645486441 -2.623749863
[400,] -0.500986882   2.321564694  1.2618729804 -1.308438505 -1.874042761
[401,]  0.404141821  -0.523101936 -0.6077818270  0.772714239  1.715825845
[402,] -0.640543501  -8.157576238 -2.8219361766  1.487057294  4.151240261
[403,] -1.800693181   4.116729169 -0.0851213698  1.566888798 -4.875282466
[404,]  3.185196501   0.520111507  1.0518299445  3.129502081  1.609747442
[405,] -1.020123600   3.049150439  1.1117514443 -1.340415412 -2.780888197
[406,]  4.309932927   0.386020583  3.4561104360 -0.445672046  5.746869838
[407,] -3.121702156   2.123491114 -2.5150757556 -0.568799765 -2.979108358
[408,] -0.015241285  -2.446383919 -1.1097423524 -0.041111957  0.976616685
[409,] -0.196127354   6.966414072 -0.0334626406  2.048691608 -2.865138005
[410,]  0.373525316  -1.171129072 -0.7893948354  2.384201274  0.382008279
[411,]  0.108900657  -0.810103386 -1.7045007460  1.370580639 -1.180496674
[412,]  0.831090826  -2.565610263  0.2406260185  0.591343194  0.566718511
[413,]  2.906847610  -3.508634991  1.2029936018  2.422803927  3.967051094
[414,] -0.596200885   0.544154904 -1.8685619219 -0.925784670 -1.460321718
[415,] -1.305128236  -0.900257464  0.4670259623  1.524989221 -2.103353381
[416,]  1.470058174   3.342565974  0.4359323038 -0.771130850 -0.205355886
[417,]  0.677411999   1.376181142  1.6843918604  2.231978526  1.214850686
[418,] -0.548564821  -2.625386782  1.6284479887 -0.261429189  1.460294863
[419,] -3.327350762   0.022108449 -2.3087874343 -0.999521916 -2.639165668
[420,]  1.519445791  -7.362959652  0.5153801542  0.170843358  5.700242905
[421,]  0.439385011   0.683238548 -0.1081730710  2.912977284  0.198953393
[422,]  3.173652638   1.334905516  3.0801652688  0.863362947  3.830287037
[423,]  1.554826821  -1.583934812  1.2402222574  2.529253528  2.020625776
[424,]  1.226966468  -3.622143929  1.0725571281 -0.727733812  3.196872654
[425,] -1.054943051  -1.483962006  0.2559551852 -3.665160226  0.457707432
[426,]  1.483368985  -2.455070045 -0.2739190190  0.600194768  0.555487717
[427,] -1.282941291   2.701874717 -0.4618195805  2.423807255 -2.688645367
[428,] -0.739740224  -2.300149467 -1.8634552325 -1.608391641  1.233255634
[429,]  0.879792731   1.779915334  1.3628628616  0.381397913  1.153251878
[430,]  1.206926607   0.351834514  0.5993600183 -0.530403360  0.770329841
[431,]  0.471988696  -2.079862164 -1.4782468514 -1.199637402  1.894200406
[432,]  1.816111523  -1.609570706  1.3713828636 -1.304475868  3.860225529
[433,]  2.032058115   3.007732914  3.1372808396 -0.096510729  0.831496344
[434,]  0.003572781  -2.103794423  0.9935247659  0.314942465  0.954295882
[435,]  0.304649305  -0.457526603  1.9407027989 -1.916293232 -0.401132173
[436,]  0.936449471   2.557489163 -0.3025254689  2.000114128 -1.721016640
[437,]  4.055775414  -3.852973083  2.1989533460 -1.639153845  5.480100733
[438,]  0.906832715   0.085631007  1.0367607106 -1.704632841  2.735549575
[439,] -0.538818434  -3.165449278 -1.5903444479 -1.561475929  1.469334707
[440,] -1.545381677   0.372687643 -2.0679884336 -1.482697067 -1.458572298
[441,]  2.075986634   3.432644046  0.3447383549 -0.926921160  1.047362044
[442,]  0.265103514   2.569231213  1.4780904026 -2.475381770 -1.099134945
[443,]  0.200342652  -1.177979470 -1.2539377069  0.586925064  0.337304403
[444,] -1.785079286  -2.712932332 -1.5051830322 -1.524511528  0.152819178
[445,] -1.035847660   3.769964122  1.1392004918  1.269082479 -3.010631109
[446,]  1.295976024   4.243907552 -0.3754817460 -1.254022767 -1.637106570
[447,]  1.017534385   1.791657397  2.0984124712 -0.960023146  0.434061837
[448,] -0.219033262   1.215302523  0.0187702129  0.588743204  0.831678318
[449,] -1.120190295  -0.104342017 -2.1823418697 -1.044332409 -1.851926295
[450,] -0.351316982  -0.148857987 -1.6688488575  2.332145883 -2.517597784
[451,] -0.749160343  -1.179728954  0.3399263399 -3.390386702  0.364363719
[452,]  3.632046943  -3.613266085  3.7538131006  1.194468037  4.472975588
[453,]  1.903897989   0.592007701  0.8835773391 -1.013941847  1.120218602
[454,]  1.223445726  -1.741374705  1.6617553337 -1.252918647  2.753773661
[455,] -0.941614169  -4.533435627 -1.8942575906  0.474935477  2.815084098
[456,]  1.637943688  -2.023704266  0.3759644672  0.409263248  1.726142271
[457,] -0.623481531   3.118016903  0.7347811503 -1.176778259 -0.106689134
[458,]  0.520036396  -1.158691092  0.4863931263 -1.568240430  0.960339971
[459,]  0.433968296  -0.446494166 -1.0810135302 -0.546000796  0.457577502
[460,]  1.413356244   5.442970300  0.4921085584  0.120897118 -1.811488549
[461,]  2.580662827   4.956397145  3.4975858006  1.850520862 -0.448516132
[462,] -2.084154746   2.581951723 -0.3829811570  0.994110721 -4.330909923
[463,] -3.758311554   0.006369473 -3.4385333584 -0.004529990 -2.090805455
[464,]  2.233372190  -1.747556571  2.3374873153  0.762338646  3.822333227
[465,] -2.421740707   3.513676229 -0.8337129443 -0.557430892 -3.593505048
[466,] -0.336660140   0.412802904 -2.3359097923 -0.160441130 -0.371951828
[467,] -1.465698569  -0.939923326 -3.0827039944 -0.380541565  0.125236603
[468,] -1.165735184   0.132146186 -1.6296004033  1.538743196 -1.302208920
[469,] -3.277372425  -0.035820041 -1.1058153964 -1.397956531 -3.358248544
[470,]  1.429328848   3.411462270  3.4115940841 -0.604048457  0.025718416
[471,]  2.453854564   4.084050610  1.8285919156  1.995518384 -1.156147694
[472,] -1.682814907   0.296375646 -0.4175517168  0.666486417 -2.056190672
[473,] -3.065298815   2.985699679 -2.7155226499  0.149130260 -4.254352864
[474,] -2.075750967  -1.188369057 -2.6244369968  0.980214350 -2.344367145
[475,]  5.755523575  -1.230789495  7.6875553303 -1.043490612  6.492020675
[476,]  0.654089682   2.155228795  0.4045567371 -0.941674482 -0.983062790
[477,]  0.693530007   2.202452402  0.5916435631  1.952856801 -2.337994314
[478,]  0.376746643   6.258486032  1.3007160622  2.038211211 -3.054021186
[479,]  2.312255171   3.246888455  3.4013795648  2.314332178  0.913937176
[480,] -2.659874345  -0.920508087 -2.6066703437 -1.455427321 -2.088144988
[481,]  2.442221589   0.165276924  2.7865796958  1.812285138  0.269527340
[482,] -0.143398850   2.358948143 -0.7263290326  0.074017618 -0.045361179
[483,]  0.096539487   1.065388304 -0.4350798345  1.028069222  0.649652264
[484,]  4.396060817   0.568427438  6.5438081039 -1.608908298  5.307983948
[485,]  1.034929454  -5.034442027  2.6886018106 -0.119318917  3.067543222
[486,]  1.261213579   2.140328931  0.4130152357 -1.252742441  0.431791809
[487,]  3.777735905   5.738368319  4.2797913614  0.336295263  0.900497936
[488,]  1.742453751  -0.509073903  0.8046036945  1.391903854  2.337276256
[489,]  0.598862123   5.602560176 -0.4892512706 -2.240561353 -1.647455187
[490,] -1.808966335   0.950544799 -0.4402703968 -1.114836792 -1.958751752
[491,]  1.147639951   3.045089771  0.6658656411 -0.468945869 -1.181528375
[492,] -0.302052105  -0.899231441 -0.9298090014  0.511887115 -0.149915393
[493,] -2.787062704   2.190348120 -0.8378460550  2.523854051 -5.293623579
[494,] -2.384163980   1.978987035 -1.5274853613 -0.027067227 -2.984903112
[495,] -2.339427452  -1.394737558 -2.9129292257 -0.592136919 -1.020572785
[496,]  2.591168272  -2.024848621  2.3191618940 -1.425546396  4.379194752
[497,] -2.315606888  -6.659263510 -0.7931744013  0.023099316  0.823265807
[498,] -1.479801643   0.199177540 -0.9951182020  2.950851013 -1.434877535
[499,] -2.598236434  -2.438550326 -2.3071029643 -1.614365594 -2.396126880
[500,]  4.684731619  -1.124806653  3.8999456391 -0.340925229  3.834046033
[501,]  1.718395624  -1.424069028 -1.1337578872 -1.438901085  2.011743809
[502,] -0.350617794  -0.650804759  0.4089715200  0.317384267 -0.475468421
[503,] -2.618672379   0.249212142 -2.7262094428 -3.266938767 -1.595143932
[504,]  1.160879361  -1.501700885 -0.3159113539  0.413675509  2.482799724
[505,] -1.058723481   1.772339794 -1.1449453640  0.528593060 -2.323751198
[506,]  1.461827652  -5.008462143  3.1742598518  0.293775466  4.729719631
[507,] -2.016149234  -0.963373341 -1.3853518956 -1.479719122 -0.591673555
[508,] -1.447306227   2.943607283 -2.5193851982  1.081540993 -4.028585710
[509,]  0.777620812   4.038731506 -1.2196686492 -0.948977677 -0.984959441
[510,] -2.274215743  -1.357999032  0.0348886473  0.920368579 -1.107955211
[511,] -3.147334257  -3.709682406 -2.2336661545 -1.154749624 -0.524728928
[512,]  2.178205062   0.118253465  2.5280841472  0.869509988  2.593017427
[513,]  0.416310461   3.263269275  1.0279936725  1.043265833 -1.982523273
[514,]  0.606487359  -1.736180222 -0.6890234684 -0.355851859  2.726844003
[515,] -0.504930141   4.796381794 -0.6793226240 -2.054837392 -3.945210391
[516,] -0.727309342   1.140961595  0.5812838120 -2.064489028 -1.053020530
[517,]  1.015979172   0.596699059  1.0168315394  0.095138489 -0.634886808
[518,]  2.177206805  -3.371566537  2.8686656322 -0.304577606  2.989458196
[519,] -5.452032484   4.376416265 -3.5403543505 -3.833492725 -6.587118634
[520,]  0.162546595   2.523109007  1.0413874246 -1.180434554  0.138107171
[521,]  1.001972200   2.148634614  0.3872183001  0.581784041 -1.770168648
[522,]  1.702542415   3.484408050  2.9648111772  3.069611588 -0.728037465
[523,]  1.402555873  -2.383322954 -1.4062534572  3.066911443  1.970978299
[524,] -2.307829558   2.975389109 -1.8844070926  0.445689800 -2.270805124
[525,]  0.668466457   4.649333212 -1.1470772328  0.349380962 -4.030893923
[526,] -1.927554531   1.176985649 -1.2100627011 -1.437482026 -2.242505659
[527,] -1.995456874   3.253252831 -0.9084933631 -0.242762166 -1.477475091
[528,] -3.251116635  -3.007832462 -3.3878926286 -2.243406869  0.226053791
[529,]  1.877642266   1.413479051  2.6512272462 -1.857131020  1.569857678
[530,] -1.049972421   0.020855932  0.2574370298 -0.343728751 -1.620806296
[531,]  2.834461246   5.016248333  3.5964275759  1.145922272 -0.453784061
[532,]  0.167344829  -3.509939864 -0.9309146546 -0.761886328  3.706009262
[533,] -3.049823154   2.236110310 -1.3448897266 -1.856894815 -3.485971384
[534,]  0.526536665   1.687006831  2.2215332469  2.672287465 -1.788753501
[535,]  1.633273340  -0.957197665  1.1389166145 -2.201998378  1.705893280
[536,] -1.139992739  -6.855312183 -1.2195266328 -0.200085706  2.581023693
[537,]  0.301942751  -0.029612354  0.1945668202  1.949447208  1.151611701
[538,]  1.290475345  -0.505408272  1.2702275229  0.092288252  1.431451260
[539,] -0.903461686  -1.702447760 -1.3842013519 -0.590125734  0.862797469
[540,] -0.172331315   5.331763579 -0.0909878019 -2.226112440 -2.545325375
[541,] -0.846052891  -2.613269793 -0.5303659582  1.921666687  0.906156217
[542,]  2.301965393   1.934757780  1.9001921097 -0.578429322  0.952667518
[543,] -1.002025349   0.588354242  0.8412176287 -0.513837550  0.637652298
[544,] -2.179999749  -3.750957639 -3.0202304397 -1.625441400 -0.663515299
[545,] -4.375896478  -0.333128177 -7.4313975399 -4.971850132 -4.514082469
[546,] -2.042611143  -0.607990125 -4.0177168229 -2.165745297 -1.679500922
[547,] -1.482734907  -1.986418551 -0.7001209317 -0.176834302  0.083835112
[548,] -1.638337586   3.242482956 -0.1670055422  0.740646580 -4.931943173
[549,] -0.833042627  -3.045476059 -1.6075459718  0.576270315  1.700108596
[550,] -0.273594983   1.778294975  0.3105959629  0.463975627 -0.711636790
[551,] -1.516525050  -2.271379529 -2.2436502280 -1.505621126  0.380675730
[552,] -1.638628048   2.175281219 -0.9591803903  0.014619932 -4.362384150
[553,] -1.170113159  -5.110500837  0.0323073667 -3.375816540  2.411445866
[554,]  1.673547179  -4.682467997  0.9679521071 -0.439571948  5.898395560
[555,] -0.689177059   5.406589849 -0.7006791297  0.009517968 -4.398153718
[556,]  1.709926311  -1.731562071  2.8129877774  1.263098166  2.317410080
[557,] -0.045110443  -2.094785463 -1.0889114387 -2.309107602  1.475913089
[558,]  3.846896184   0.442567230  3.9678876903  3.036868241  3.514141829
[559,] -0.690657480  -2.630552256 -0.8036164400  0.967667758  1.143878452
[560,]  2.480899813  -1.667020830  2.7925722470  1.663729333  1.489271092
[561,]  2.816367686  -1.821895999  3.7001080740  2.194194189  3.786378457
[562,]  1.527166114   3.330686868  1.1913812338 -1.102823796 -1.635001920
[563,] -0.993053446   1.719644698 -1.6143464201 -2.333644030 -0.484120492
[564,] -0.529606899   1.881044865 -1.1785461855  2.369653360 -3.317146095
[565,]  2.390283371   0.145878063  0.6613034887 -0.138057843  2.469828457
[566,]  1.129896823   2.585141888  2.7704905597  1.592002101  1.220857443
[567,]  1.785546031  -4.862551967  1.2474729044 -3.032824828  6.235628315
[568,] -0.359694320   1.086528394  2.0351272609 -1.946065522 -1.179936196
[569,] -1.460835217  -4.781142448  0.4166102951 -0.095297087  1.184300782
[570,]  1.180775426  -4.302735910  0.7851905185 -0.120296399  4.253901732
[571,]  0.824557804   1.079477161 -0.1133365775 -1.660113364  0.846366667
[572,] -0.439714482  -1.803001446 -0.0283155285  2.189670321  0.226888826
[573,] -0.396525959   2.966841234  1.6890565322 -1.492667968 -1.476711061
[574,]  0.640134316   2.062480839  0.6155279504  0.385553537  0.458440356
[575,]  0.435915340  -3.867260452 -0.3556649542  0.437707854  3.238045410
[576,] -1.586046425  -5.562491172 -2.3182650616  1.757738944  3.550290232
[577,]  3.734796245  -1.939343130  4.3159695238  0.923435123  4.448483743
[578,]  0.107885766  -4.096395609 -0.4058648615 -3.178922727  2.959180706
[579,]  0.292150660   6.650027940  1.7147024926  3.220320859 -2.993357343
[580,] -1.163388542  -0.490282694  1.7005644700 -2.278423975 -0.497236815
[581,] -0.254436813  -4.771251632 -0.1272173523 -2.161902706  3.328925880
[582,]  2.286309408  -4.989246450  3.7939482052 -0.939825615  4.845695279
[583,]  5.300738589   3.675560080  4.1071046641  1.124470977  3.068312370
[584,] -1.144907560  -4.627227342 -1.3567687296 -0.636048634  1.082324308
[585,]  1.435481821  -6.273809796 -0.4600575320  0.316032451  5.853268555
[586,] -0.044150599  -5.779164097 -2.4050619238  1.097357495  3.807099962
[587,] -0.884418342   1.609573560 -0.1043158725 -0.906490847 -0.420072700
[588,]  1.871024220   3.038820318  1.8631497672  0.083433549 -0.962048310
[589,]  1.597289267   2.152264536  1.3053984303  0.683902708 -0.215311968
[590,] -1.684722729   0.767191432 -2.2342991074 -0.587728084 -1.316599312
[591,]  0.901063871  -1.361110323  0.7593274019 -0.895492892  1.415411214
[592,]  0.060863767   2.326388765 -0.0247762495  1.103084471 -1.151651297
[593,] -1.253942822   2.222581918 -0.3786485312 -0.485735841 -1.549899438
[594,]  2.551567369   3.321406275  1.8981325762 -1.717846394  1.397453980
[595,] -1.666703405  -1.450441056 -1.6259547328 -1.121364885 -0.421868125
[596,]  0.985996824   1.598265899 -0.7873917128  1.118490957 -1.670996235
[597,] -2.054492328   3.100898453 -1.9596489969  1.536128709 -3.205598502
[598,]  0.614683422   9.865447526 -1.4644369565  2.428347437 -5.693914819
[599,] -3.969505747  -0.185765213 -5.1041285887 -1.198211593 -4.289183849
[600,] -1.249656772  -3.620582850 -3.0971775281 -2.169322322  0.332091834
[601,]  2.807177819  -7.311841229  2.4925955012  1.431159760  5.842353201
[602,]  1.697350714   4.211742396  3.3581603092  1.660672852 -1.513427914
[603,]  1.474316567  -0.172098358  1.5186415113 -1.564820007  2.200312026
[604,]  2.515899542   2.666596800  2.7525894498 -0.064088397  2.057706557
[605,] -0.686808814  -6.591696447 -3.9454376130  0.949991578  1.792121113
[606,]  2.443595233  -0.837717552  3.0713582766 -1.799166146  4.127095734
[607,] -0.062316217   4.331168021 -0.6765890736  0.389575183 -2.952638179
[608,]  2.392137234  -1.822984291  3.0146553005  2.315240582  1.412687147
[609,]  1.405128553   1.413225465  3.1081923843 -0.673972513  0.682742213
[610,]  1.598438524  -0.598902471  1.1318735989 -0.392652065  1.358472836
[611,]  2.945996741   0.219922255  3.3446295689  0.855022546  1.575001311
[612,]  2.441540333   1.432535916  1.6889227906 -0.041138774  0.517575220
[613,]  0.662791737  -4.206601401  1.7852459379 -2.960369668  2.335310826
[614,]  1.409769982   3.246361620  2.3872913793  0.879130992 -1.186169693
[615,] -3.713201200   1.554657929 -1.7953499868  2.890796982 -3.910438316
[616,] -0.367605547  -1.811389413 -1.4817665532 -1.057084934  0.358328835
[617,] -1.814956955   2.529764189 -2.6889061331 -0.867660171 -3.302865613
[618,] -2.250927305  -3.630320097 -3.3203844063  2.265118517 -0.567337065
[619,]  3.185804668   0.076784835  2.2794335947  0.598405700  3.012051340
[620,]  0.612485037  -3.639865941 -0.9873471322 -2.434328512  2.013078548
[621,] -2.666346316   3.408414503 -1.0137783546 -1.128261358 -4.107276833
[622,] -0.277152361  -1.555063943  1.5449620638 -2.096654772  1.628745241
[623,] -0.564162497  -0.794000344  0.9531111685  2.767378955 -0.320098658
[624,] -1.558464987  -0.680002369 -0.7668893262  2.199817044 -2.443047178
[625,]  0.962794783   4.184693190  1.4998903575  0.718249904 -1.584003578
[626,]  0.920494353   3.025654055  1.5866716514 -0.071368049 -0.537474981
[627,] -0.792606159   0.781473720 -1.2125693788 -1.549729192 -2.240582619
[628,] -1.272960629   2.143286753 -1.4465666848  1.017271451 -3.279670121
[629,]  0.300454258   7.698656926  0.3007005469  0.610739357 -3.616025233
[630,] -1.405415753  -3.151886648 -1.9268935361 -0.454877748  0.980769113
[631,]  4.912622179  -5.162062540  3.3426931963  0.003814061  8.589078840
[632,]  0.480629541  -2.536661228 -0.6980510994 -1.926167438  2.228495218
[633,] -0.590130968  -3.352279868 -0.0612257228 -2.107486348  1.111808882
[634,] -1.054843156  -4.083860206  0.3854095851 -4.061332503  1.328373369
[635,]  0.217405937  -2.777682942 -0.4668315493  0.663442913  2.787685290
[636,]  2.243994751   0.768431800  4.1809990679 -1.948713096  1.937764956
[637,] -0.273548320  -4.550628956  0.4570303741 -0.952786547  2.156153931
[638,]  1.576767296   1.948546702  1.6215404919  0.144540619  1.176143143
[639,]  0.120669640  -0.871217197  0.9920319593  1.927565469  0.609587247
[640,]  0.906207717  -2.128306809 -0.7146282865  0.304020170  0.781921269
[641,] -3.456747208   5.961674485 -2.7179913632 -1.618433491 -3.999603301
[642,] -0.325269685   0.867991765 -1.4770639820 -0.781287706  1.017122843
[643,]  0.853920969   3.209140303  1.2683764805  1.909054788 -1.809231827
[644,]  1.558159195  -0.014461332 -0.7652724412  0.363188258  1.694745774
[645,]  1.919514333  -4.360346040 -0.0733399533  2.738784661  4.016452586
[646,] -1.158946360   6.516551674 -1.2714695126  1.135282199 -3.741180003
[647,] -2.347792498   3.281993171 -2.3536491530 -0.013153782 -4.550046070
[648,]  0.153184558   2.657409893 -0.7834774795  0.586570719  0.264892327
[649,] -0.075421816  -0.775408755  0.9653533443  1.924853136  0.554536603
[650,] -1.954970366   1.060069478 -2.9805575173 -0.252732875 -2.188717240
[651,] -0.986370464  -0.191289347 -2.5261327038 -1.063989770 -2.287301636
[652,]  1.989891273   4.020737668  0.8597697458 -3.512436279  0.652772052
[653,] -0.655660850   1.482542245  0.0233291917 -0.131102960 -2.798189534
[654,]  2.269737641   4.741250449  3.3258732971  1.984517982 -2.110310246
[655,] -3.738487292   3.763709588 -0.5226914389 -0.845292758 -4.760786178
[656,]  2.892186440   7.797695345  2.0448110862  0.215256879 -1.971526044
[657,] -4.041574525   4.242700420 -4.2601662544 -0.228402007 -6.104531143
[658,] -1.271420268  -3.965343661 -2.1320814150  0.598475261  1.500802383
[659,]  1.480175536  -0.305564529  1.2536996029 -0.961733524  1.814771745
[660,] -0.817766506  -0.573934864 -1.0385565888  2.462949563  0.549626275
[661,] -1.731724241   2.143396807 -0.9526306498  0.945961229 -1.830294342
[662,]  0.585149091   1.461671951  0.7586155195  1.900448974  0.492164620
[663,] -0.495709609   0.795564790  1.4584312385  1.649826683 -1.043117067
[664,]  0.470755908   1.296015036 -0.3172529956  0.918461943 -0.546944680
[665,]  3.457516906   2.890464758  3.4268584058  2.073116618  0.514225756
[666,] -0.984814535   5.190223662 -0.3571826800 -1.048280531 -2.473581727
[667,] -1.985051719   3.423216273 -0.1113935773 -0.557473528 -2.640056157
[668,]  1.832628587  -0.117776083  0.0356676424 -0.790898757  1.020076957
[669,]  0.601267502  -4.358435023  0.7933610446 -0.396714879  3.265145308
[670,]  1.713910830  -4.031945959  0.1570854388  1.882177065  3.383697621
[671,] -1.307008249   0.038764740 -0.2964032984  0.191008923  1.078037819
[672,] -1.358350068   9.230426208 -1.9648817632 -0.483711310 -6.263034644
[673,]  0.018399244  -6.301507331  0.0794882448 -0.178218998  2.943565818
[674,]  1.915898072   3.436674845  2.0127391342  0.843292647  0.556229629
[675,] -0.432599053   2.920379926  0.7948193036  0.204437662 -1.689926517
[676,] -1.992099908  -4.525137984 -1.4350708342 -0.494388714  0.711151351
[677,]  1.242120300  -2.573578856  0.5977245908  1.282401111  2.942843917
[678,] -1.145095487  -1.678504543  1.3899078712  0.656994549  0.332112827
[679,] -0.015690354  -4.393341040 -0.4095047587 -1.014950847  1.609557688
[680,] -0.056740498   7.246872133  0.9277882301  0.417821650 -2.964592327
[681,]  2.863001528   5.815099590  4.2050426473  0.034795117  0.722651793
[682,] -0.633851638   0.325610917 -2.9581543536  1.453603804 -1.561752756
[683,]  0.559220504   3.795815984 -0.3743257713  1.907908674 -0.810082219
[684,]  4.290333916   3.038887492  5.4016954071  1.574921554  3.273358137
[685,] -0.186794961  -0.405063155 -1.2739476718 -1.339539254  0.516570883
[686,] -3.106455899   1.705549934 -0.1634416855  0.142874203 -2.439664505
[687,]  1.130368008  -0.403120147  2.7368152317 -0.872782742  2.257992772
[688,]  0.395712921  -0.467215972 -0.0383191199  1.596814717 -1.392959400
[689,] -0.442939192   3.166005875  0.8747197733 -1.350961189 -0.973520549
[690,]  0.931682377  -3.752253291 -0.4845000549 -0.939445655  4.636790750
[691,]  3.290438502  -1.017672977  3.6825304889  0.337282254  4.678899145
[692,] -0.129510032  -0.505534205 -1.0205582578 -0.851913840  0.400727476
[693,]  0.514182192  -0.034610868  0.4721612257  1.033964641  0.949413118
[694,] -3.789008591   5.568587405 -2.3942592927 -0.491288004 -6.019993440
[695,]  2.271818612  -3.803097061  2.2211626074  0.703865127  4.922844826
[696,] -1.802301919  -0.794704670 -2.6342686374  0.641338030 -0.967936776
[697,] -2.195866141   2.025110636 -4.3008245008  3.356712750 -3.704175123
[698,]  0.214279931   1.036613828  0.6056408562 -1.481286239 -0.424443406
[699,] -4.365172549   0.187278757 -3.6618565528 -1.111036203 -3.903103259
[700,] -2.323497343  -0.320817663 -3.4068608234  1.607645459 -2.465016854
[701,]  1.044393633  -0.428398260  2.9262018472  1.567778654  0.999212100
[702,]  3.013580982  -0.796684430  2.7426273978  0.782634764  3.935755587
[703,] -0.070239091  -0.912825248 -0.1166317544 -0.470468335  1.254705196
[704,] -0.059963352   1.973892108  2.5170468154 -0.517458996  0.944609703
[705,]  2.590800760   5.333389430  2.5248245784  1.098300755 -3.161284126
[706,] -1.937571459  -3.545243310 -1.5855075438  0.340579276  0.118113564
[707,] -0.887914411  -3.483447899  0.1522108967 -1.231990784  3.269438095
[708,] -1.993934244  -0.756633361 -0.9638599280 -0.295013342 -0.321252029
[709,]  1.360879924  -0.424875101  0.7919573152 -1.024809312 -0.740646124
[710,]  0.676872148  -7.387174185  0.3715712062 -0.274664031  5.728938905
[711,] -0.449568822   2.127826640 -0.6750041341  0.432909702 -1.223875281
[712,] -0.353136069  -1.243011625  1.2526274956 -0.723871149  2.453420748
[713,]  0.271822297  -2.757667594 -0.0268512381 -1.993731689  2.818838556
[714,]  1.744307788  -1.702890373  0.1769017192 -0.879694049  4.175886178
[715,]  0.845854851   3.529937516  1.7911764224  1.682764062 -2.805446422
[716,]  0.911507841   5.394251355  1.7834682476  0.909599944 -2.869235165
[717,] -1.589666866   1.066421762 -1.6788048738  0.233739233 -1.984749517
[718,]  0.072676669   0.620257578  0.9230352605  0.367181712  0.875466143
[719,]  2.733833857  -1.627313953  1.3777588314  1.796585327  3.816847293
[720,] -1.567208031  -5.512627795 -2.3589619366 -0.007309992  0.875208551
[721,] -2.198977014  -7.540483938 -2.2145260112 -4.026983014  3.076707518
[722,]  0.211824342  -2.961546107  0.3984818951 -0.127663410  0.949887745
[723,] -0.242084838  -5.963216409 -1.9744710270 -1.140494379  4.379596896
[724,]  1.450507863   3.824337152  1.5752727914 -1.464988140  0.230288676
[725,] -0.330011897   7.797341487  0.2243535543  0.339937203 -5.222482640
[726,]  0.967186902   1.027704996 -0.9749764273  0.679918304  1.173466260
[727,]  0.310260769  -3.555189209  1.6184844055  1.303855692  2.857119038
[728,]  1.499787743  -6.809780038  0.5587978299  0.110214460  5.836555775
[729,] -1.554679937   3.593510241  0.0369922555 -0.640463796 -3.605683613
[730,] -2.003823483  -2.581327508  0.6285289255  0.028555318 -2.207076533
[731,] -0.227649281   1.886149333 -2.0845628002 -0.684688475 -1.788037190
[732,] -1.582911001   1.720968334 -1.0879294626 -4.951521331 -2.197972639
[733,] -0.476707192   4.674850588  0.8798566315 -2.640304340 -1.952383371
[734,] -0.402301478   2.004681409  0.4213662970  2.062354522 -0.761396983
[735,]  2.376127983  -3.072144128  2.6275935762  0.357626250  4.638671921
[736,]  0.881755693   1.304887493  1.6105670523  0.101049133  1.091238517
[737,] -2.240316474  -1.251649511 -1.6553466360  0.551624995 -0.427749905
[738,] -0.248697575  -2.384506310 -1.2135470015 -0.653784618  0.314643520
[739,]  0.434364025  -5.812317943  0.0620217349  2.204676077  3.246828405
[740,] -1.511312333   5.426706244  0.0307704846 -1.667977429 -2.912716168
[741,]  0.378134117  -2.227811270  0.8673000798  1.490642741  1.959776491
[742,] -0.885010370   1.584634889 -0.7366072092 -1.557049495 -0.234626007
[743,]  1.037919328  -2.457155326  1.5768421021  2.019570436  3.975374884
[744,]  1.690845925  -3.764061109  3.1130911387  1.951504874  2.618796393
[745,]  1.176998570   2.569612740  2.4960547322 -0.097386131 -0.160344469
[746,]  2.659626184   3.976778884  1.7595454730  0.823788730 -0.133357933
[747,] -2.269270463  -0.427690104 -1.3528046972 -0.698679557 -1.751949714
[748,]  2.178331713   1.134928081  1.3621819551  0.274208634  1.500396212
[749,]  1.197797085   0.133440675 -0.0803064531 -0.351949701  0.619688874
[750,] -1.104484541  -1.084351911  0.7691604146 -0.247406131 -1.219955039
[751,] -0.329593596   3.939781632  0.9424645150  1.206475499 -2.867678260
[752,]  0.010858794   0.028398559 -0.2107277207 -1.885100717  0.001396571
[753,] -1.479616958   2.731647357 -2.4164217671 -1.014898113 -2.085992305
[754,]  0.195745915  -3.202733305  0.8415752331 -0.079580354  2.632328641
[755,]  0.191399766   5.355333889  0.8862173860  0.076566862 -2.965684458
[756,]  2.305663578   1.407666944  2.3718010880 -1.793925219  0.311560361
[757,]  0.439864986   2.351877113 -0.2832473031  1.614253890 -1.792755643
[758,]  0.240579365   1.352740126 -0.8273565995  1.926354533 -1.276892509
[759,]  2.008900463  -2.843462241  1.0854392540 -0.132149892  2.756883438
[760,]  1.707347249  -3.547018578  1.3319123670 -0.313639326  3.289679363
[761,] -2.306343178   2.134064878  0.7435479682  0.259133250 -3.731965137
[762,]  3.005897520   1.962585636  3.9217290971  2.106037662  2.003094511
[763,] -2.624116364  -3.637138243 -2.5080971034 -0.413621980 -1.376724426
[764,]  0.779825304  -1.660827361  0.2163319832 -1.504609014  2.115022659
[765,]  0.909121331   0.193773300  1.4828518438  1.124248351  0.604931725
[766,]  1.253796757  -7.152543640  0.9778368409  0.789846536  4.200889320
[767,] -0.411597350  -6.303594948  0.7468999138  4.766816000  2.984946174
[768,] -0.698006025   4.472618910  0.0041807399  2.023273101 -2.230073237
[769,]  0.610022588  -0.196092967  2.7154887360  1.335621975  1.704899592
[770,] -1.209949272  -1.565391353 -0.0571311044 -3.515267148  1.283125706
[771,]  1.750078702  -4.218703627  1.4554044927  1.299355291  3.078215404
[772,]  0.507448900   1.283153025  1.3666862658 -1.219305136 -0.996568105
[773,]  1.540318056  -0.633239937  1.1932147476  1.851560294  3.300948840
[774,]  1.529795364   1.285885393  1.2113249076  0.131788086  0.441395772
[775,]  3.260148533   4.869505397  3.9457963782  0.518121243  0.815770297
[776,]  0.186580617  -1.818970504  1.9457621155 -1.509474471  1.806481046
[777,] -0.686455425   0.133385108 -0.9887047819 -0.933247930 -1.654272574
[778,] -1.568897702   1.290683715 -0.8589525665  0.568819816 -1.983675957
[779,] -0.395123967   4.772802593 -0.4663853558 -1.336303626 -2.201785096
[780,]  2.819021073   3.495155571  5.2626436958  1.529406919 -1.054558292
[781,]  1.080017687  -1.055016411  0.0911826557  1.769497663  1.431831269
[782,] -1.041286385   2.029462474 -2.7268325315 -2.684900070 -1.238199439
[783,] -2.359415274   0.692818593 -2.8617585620  0.028006792 -3.094883315
[784,] -0.434176986   1.416043164 -1.3438242357 -3.381264451 -0.584957352
[785,]  1.502935017  -9.895691883 -0.2830694337  2.035493488  6.868280944
[786,]  2.271699898   2.253528334  1.4955602255 -0.351720562  0.230735265
[787,] -2.162031903  -0.934477157 -0.6172285853 -2.283534170 -0.848020956
[788,] -4.117280669   0.953345090 -2.0811057689 -3.282284031 -4.499898820
[789,]  0.520669030  -2.724844099 -0.5389557936 -0.408134873  2.321375000
[790,]  3.373968880   0.417418920  1.3588290420  0.726449764  1.927558995
[791,] -2.362530279   2.652603444 -2.6537653916 -1.873838603 -2.942545393
[792,]  1.325366584  -1.417318168  0.0337180707  0.503419175  2.077593583
[793,] -0.846683806  -3.028072687 -1.2302037198  2.615584927  0.481792688
[794,]  0.150650384   0.063978018  0.6704218708  0.803128273 -0.225200563
[795,] -2.085416032   2.423736711 -3.5296847657 -0.477660067 -3.757569896
[796,] -3.696173237   0.996584657 -2.5646429761 -0.329361570 -4.232942081
[797,] -0.851170496   0.165851719 -0.2157931645 -0.149200861 -0.345746086
[798,]  4.278500744   1.609834226  3.9205400300 -1.277134949  3.291661044
[799,] -1.747572181   4.054734706 -1.3492809247 -1.534892400 -3.305015726
[800,]  0.498296810   0.506639979  0.6839018920 -1.723128274  0.169399823

# Correlation of decisionmaker random slopes
Omega <- cor(matrix(rnorm(P*(P+2)), P+2, P))

# Scale of decisionmaker random slopes
tau <- abs(rnorm(P, 0, .5))

# Covariance matrix of decisionmaker random slopes
Sigma <- diag(tau) %*% Omega %*% diag(tau)

# Centers of random slopes
beta <- rnorm(P)

# Individual slopes
beta_i <- MASS::mvrnorm(I, beta, Sigma) + W %*% t(Gamma)

# Again, quick plot to sanity check
plot(as.data.frame(beta_i))

# Create X -- let's make this a dummy matrix
X <- matrix(sample(0:1, Tasks*I*J*P, replace = T), Tasks*I*J, P)
# Each of the rows in this matrix correspond to a choice presented to a given individual
# in a given task

indexes <- crossing(individual = 1:I, task = 1:Tasks, option = 1:J) %>% 
  mutate(row = 1:n())

# Write a Gumbel random number generator using inverse CDF trick
rgumbel <- function(n, mu = 0, beta = 1) mu - beta * log(-log(runif(n)))
mean(rgumbel(1e6))

[1] 0.576678

# Ok, now we need to simulate choices. Each person in each task compares each 
# choice according to X*beta_i + epsilon, where epsilon is gumbel distributed. 
# They return their rankings. 

ranked_options <- indexes %>% 
  group_by(individual, task) %>% 
  mutate(fixed_utility = as.numeric(X[row,] %*% as.numeric(beta_i[first(individual),])),
         plus_gumbel_error = fixed_utility + rgumbel(n()),
         rank = rank(-plus_gumbel_error),
         # We're going to use the order rather than the rank in the Stan part of the model
         order = order(rank),
         # And here we create a dummy vector for the best choice
         best_choice = as.numeric(1:n() == which.max(plus_gumbel_error)),
         worst_choice = as.numeric(1:n() == which.min(plus_gumbel_error))
  )


tt <- ranked_options %>% 
  group_by(individual, task) %>%
  summarise(start = min(row), 
            end = max(row)) %>% 
  ungroup %>%
  mutate(task_number = 1:n())

`summarise()` has grouped output by 'individual'. You can override using the
`.groups` argument.

While this simulation isn’t exactly simple, let me pause to point out something I’ve not done. I have not made any choices here that are designed to make MaxDiff look poor. Since we want to use a model of realistic complexity, this has simulated structure such that it’ll benefit from the random coefficients formulation I’ll use in a moment. Similarly, this has load-bearing covariates at both the individual and choice level, since we would almost always have those in the real world. But the core utility function we’ve created for respondents boils down to simulating data under the assumptions of the multinomial choice model with Gumbel errors we’ve been discussing all along.

Let’s specify the best choice model in Stan⁷:

Best choice logit model

best <- "// Again, note that this code is just reproducing https://khakieconomics.github.io/2018/12/27/Ranked-random-coefficients-logit.html
data {
  int N; // number of rows
  int T; // number of inidvidual-choice sets/task combinations
  int I; // number of Individuals
  int P; // number of covariates that vary by choice
  int P2; // number of covariates that vary by individual
  int K; // number of choices
  
  vector<lower = 0, upper = 1>[N] choice; // binary indicator for choice
  matrix[N, P] X; // choice attributes
  matrix[I, P2] X2; // individual attributes
  
  int task[T]; // index for tasks
  int task_individual[T]; // index for individual
  int start[T]; // the starting observation for each task
  int end[T]; // the ending observation for each task
}
parameters {
  vector[P] beta; // hypermeans of the part-worths
  matrix[P, P2] Gamma; // coefficient matrix on individual attributes
  vector<lower = 0>[P] tau; // diagonal of the part-worth covariance matrix
  matrix[I, P] z; // individual random effects (unscaled)
  cholesky_factor_corr[P] L_Omega; // the cholesky factor of the correlation matrix of tastes/part-worths
}
transformed parameters {
  // here we use the reparameterization discussed on slide 30
  matrix[I, P] beta_individual = rep_matrix(beta', I) + X2 * Gamma' + z*diag_pre_multiply(tau, L_Omega);
}
model {
  // create a temporary holding vector
  vector[N] log_prob;
  
  // priors on the parameters
  tau ~ normal(0, .5);
  beta ~ normal(0, .5);
  to_vector(z) ~ normal(0, 1);
  L_Omega ~ lkj_corr_cholesky(4);
  to_vector(Gamma) ~ normal(0, 1);
  
  // log probabilities of each choice in the dataset
  for(t in 1:T) {
    vector[K] utilities; // tmp vector holding the utilities for the task/individual combination
    // add utility from product attributes with individual part-worths/marginal utilities
    utilities = X[start[t]:end[t]]*beta_individual[task_individual[t]]';
    
    log_prob[start[t]:end[t]] = log_softmax(utilities);
  }
  
  target += log_prob' * choice;
}"

Much of this logic handles the multilevel component and covariates, but the real key is the last few lines where the likelihood is built. That’s where the problem will arise when we extend the above model into MaxDiff:

MaxDiff model

best_worst <- "// Again, note that this code is starting from https://khakieconomics.github.io/2018/12/27/Ranked-random-coefficients-logit.html
// I've just adapted to to include the worst choice
data {
  int N; // number of rows
  int T; // number of inidvidual-choice sets/task combinations
  int I; // number of Individuals
  int P; // number of covariates that vary by choice
  int P2; // number of covariates that vary by individual
  int K; // number of choices
  
  vector<lower = 0, upper = 1>[N] choice; // binary indicator for choice
  vector<lower = 0, upper = 1>[N] worst_choice; // binary indicator for worst choice
  matrix[N, P] X; // choice attributes
  matrix[I, P2] X2; // individual attributes
  
  int task[T]; // index for tasks
  int task_individual[T]; // index for individual
  int start[T]; // the starting observation for each task
  int end[T]; // the ending observation for each task
}
parameters {
  vector[P] beta; // hypermeans of the part-worths
  matrix[P, P2] Gamma; // coefficient matrix on individual attributes
  vector<lower = 0>[P] tau; // diagonal of the part-worth covariance matrix
  matrix[I, P] z; // individual random effects (unscaled)
  cholesky_factor_corr[P] L_Omega; // the cholesky factor of the correlation matrix of tastes/part-worths
}
transformed parameters {
  // here we use the reparameterization discussed on slide 30
  matrix[I, P] beta_individual = rep_matrix(beta', I) + X2 * Gamma' + z*diag_pre_multiply(tau, L_Omega);
}
model {
  // create a temporary holding vector
  vector[N] log_prob;
  vector[N] log_prob_worst;
  
  // priors on the parameters
  tau ~ normal(0, .5);
  beta ~ normal(0, .5);
  to_vector(z) ~ normal(0, 1);
  L_Omega ~ lkj_corr_cholesky(4);
  to_vector(Gamma) ~ normal(0, 1);
  
  // log probabilities of each choice in the dataset
  for(t in 1:T) {
    vector[K] utilities; // tmp vector holding the utilities for the task/individual combination
    // add utility from product attributes with individual part-worths/marginal utilities
    utilities = X[start[t]:end[t]]*beta_individual[task_individual[t]]';
    
    log_prob[start[t]:end[t]] = log_softmax(utilities);
    log_prob_worst[start[t]:end[t]] = log_softmax(-utilities); // Here's the problem!
  }
  
  target += log_prob' * choice;
  target += log_prob_worst' * worst_choice;
}"

Comparing the two, you can definitely see the implementation simplicity MaxDiff provides; we’re essentially just building out log_prob_worst and pulling it into the likelihood. Boom, done.

So how does this perform? Well, let’s fit the two and find out. The next code block has some plumbing to fit the models and extract results, so feel free to jump to the results plot below:

Fit the two models, and gather predictions to compare

# Again, note that this code is largely just reproducing https://khakieconomics.github.io/2018/12/27/Ranked-random-coefficients-logit.html,
# and then adding in details need to compare it to MaxDiff

data_list_best_choice <-list(N = nrow(X),
                             T = nrow(tt),
                             I = I, 
                             P = P, 
                             P2 = P2, 
                             K = J, 
                             # NOTE!! This is the tricky bit -- we use the order of the ranks (within task)
                             # Not the raw rank orderings. This is how we get the likelihood evaluation to be pretty quick
                             choice = ranked_options$best_choice,
                             X = X, 
                             X2 = W, 
                             task = tt$task_number, 
                             task_individual = tt$individual,
                             start = tt$start, 
                             end = tt$end) 

data_list_best_worst <-list(N = nrow(X),
                             T = nrow(tt),
                             I = I, 
                             P = P, 
                             P2 = P2, 
                             K = J,
                             choice = ranked_options$best_choice,
                             worst_choice = ranked_options$worst_choice,
                             X = X, 
                             X2 = W, 
                             task = tt$task_number, 
                             task_individual = tt$individual,
                             start = tt$start, 
                             end = tt$end) 

compiled_best_choice_model <- stan_model(model_code = best)

Warning in readLines(file.path(makevar_files)): incomplete final line found on
'C:/Users/timma/Documents/.R/Makevars.win'

compiled_bw_choice_model <- stan_model(model_code = best_worst)

Warning in readLines(file.path(makevar_files)): incomplete final line found on
'C:/Users/timma/Documents/.R/Makevars.win'

best_choice_fit <- sampling(compiled_best_choice_model, 
                            data = data_list_best_choice, 
                            iter = 800)

Warning: There were 15 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

best__worst_choice_fit <- sampling(compiled_bw_choice_model, 
                            data = data_list_best_worst, 
                            iter = 800)

Warning: There were 2 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

# Now make some predictions to compare
best_choice <- as.data.frame(best_choice_fit, pars = "beta_individual") %>%
  gather(Parameter, Value) %>%
  group_by(Parameter) %>%
  summarise(median = median(Value),
            lower = quantile(Value, .05),
            upper = quantile(Value, .95)) %>%
  mutate(individual = str_extract(Parameter, "[0-9]+(?=,)") %>% parse_number,
         column = str_extract(Parameter, ",[0-9]{1,2}") %>% parse_number) %>%
  arrange(individual, column) %>%
  mutate(`True value` = as.numeric(t(beta_i)),
         Dataset = "Best Choice")

best_worst_choice <- as.data.frame(best__worst_choice_fit, pars = "beta_individual") %>%
  gather(Parameter, Value) %>%
  group_by(Parameter) %>%
  summarise(median = median(Value),
            lower = quantile(Value, .05),
            upper = quantile(Value, .95)) %>%
  mutate(individual = str_extract(Parameter, "[0-9]+(?=,)") %>% parse_number,
         column = str_extract(Parameter, ",[0-9]{1,2}") %>% parse_number) %>%
  arrange(individual, column) %>%
  mutate(`True value` = as.numeric(t(beta_i)),
         Dataset = "MaxDiff")

# Combine the datasets
combined_data <- bind_rows(best_choice, best_worst_choice)

# Plot results
ggplot(combined_data, aes(x = `True value`, y = median, color = Dataset)) +
  geom_linerange(aes(ymin = lower, ymax = upper), alpha = 0.3) +
  geom_point(aes(y = median), alpha = 0.5) +
  geom_abline(intercept = 0, slope = 1) +
  labs(y = "Utility Estimates",
       title = "Utility Estimates from the Two Models",
       subtitle = "With 90% CIs") +
  scale_color_manual(values = c("Best Choice" = "blue", "MaxDiff" = "red")) +
  facet_wrap(~Dataset) +
  theme(legend.position="none")

To be clear, what’s plotted here are respondent level utilities for each choice. Concerningly, it doesn’t even seem clear that the best-worst choice model performs better here despite the extra data. The fact that we have to do this at all doesn’t bode well, but let’s compare the RMSE’s of the predictions versus the (simulated and thus known) truth to confirm:

# Is the MaxDiff model actually better?
combined_data %>%
  group_by(Dataset) %>%
  summarize(RMSE = sqrt(mean((`True value` - median)^2)))

# A tibble: 2 × 2
  Dataset      RMSE
  <chr>       <dbl>
1 Best Choice 0.218
2 MaxDiff     0.333

To say this more emphatically, MaxDiff gets twice as much data from each respondent, and it’s not even consistently the case that it is better⁸! And again, it’s not like I chose some obscure DGP where MaxDiff is uniquely poor, we’re working with data that aligns precisely with the model assumptions.

So to summarize, I’ve described some psychological and mathematical issues with MaxDiff. Then, via a comparison to a best choice only model, I showed that these differences are more than theoretical naggles: for pretty vanilla data, the MaxDiff model can actually perform worse than the best choice formulation.

It’s worth pointing out that we only get this type of precise comparison in a simulation study. On real world survey data with a similar DGP, you’d just see the models give fairly similar answers, with no reason to suspect the slight differences might actually be bad tweaks on MaxDiff’s part. With simulation, however, we can look closer, and the flaws I’ve discussed become increasingly apparent.

Alternative #1: rank-ordered logits with connected graphs of choices

Ok, so you get it. MaxDiff isn’t great. But respondents are still expensive, and you still have a choice modeling task to complete at work, and you really, really want to do better than just using the best choices. Maybe your boss thinks Sawtooth is fine since it’s been used for a long time. What’s the alternative?

In this section, I’ll lay out what we’ll call the rank-ordered logit with a connected graph of choices. If this sounds complicated, don’t worry; the model isn’t actually going to get much more complicated.

Instead, the plan is mostly to be cleverer about how we build our choice sets in the first place.

Explaining the “connected graphs”

What we really want, I’d argue, is full choice ranking data, for each respondent. But of course, respondents rapidly become unhappy and start satisficing when they are forced to rank 20 different largely similar options at the same time. How can we get the simplicity from the respondent POV of MaxDiff, but the rich ranking data we want?

The clever trick here is that by asking the right 3 item choice sets, and asking the survey taker to pick best and worst, we can collect every choice needed to construct a full ranking. We’ll then be able to simply fit a rank-ordered logit model to these data, which makes far fewer simplifying assumptions about how the choices and their utilities interact.

There are some downsides I’ll be transparent about here: for K choices, you need each respondent to do K-1 choice tasks. If you have 20 options, this can get tiring fast, and you won’t really be gathering “extra” opinions per respondent anymore by doing the minimum K-1 rounds. Here’s a hand worked example of what this looks like for intuition:

Illustration: K-1 comparisons are needed

Let’s say we want to get someone’s ranking over 2024 presidential candidate choices, from the set {Biden,Trump,RFK,Stein}.

Choice Set 1:

Options: {Biden, Stein, RFK} Best choice: Biden Worst choice: RFK

Choice Set 2:

Options: {Stein, RFK, Trump} Best choice: Stein Worst choice: Trump

Choice Set 3:

Options: {Biden, Stein, Trump} Best choice: Biden Worst choice: Trump

Now, let’s reconstruct the full ranking from these choice sets:

From Choice Set 1, we know that Biden > Stein > RFK. From Choice Set 2, we know that Stein > RFK > Trump. From Choice Set 3, we know that Biden > Stein > Trump.

Combining these partial rankings, we can infer the full ranking: Biden > Stein > RFK > Trump

You’ll notice that you can’t get the full ranking for e.g. 4 options in any less than 3 sets, but you well could add more sets to gain some statistical efficiency if you so choose. Of course, in larger choice sets, you’ll probably end up closer to K-1 sets given respondent fatigue, but for smaller sets it’s a very helpful option to know exists.

Most of the innovation here was in the choice set design, but there’s some in the model too, so let’s turn to that. The likelihood for the Rank-Ordered Logit model just extends the intuition from the best-choice model to account for the full ranking.

For a given ranking r of J alternatives, the likelihood takes the form:

P(r) = \prod_{j=1}^{J-1} \frac{\exp(\mu_{ij})}{\sum_{k=j}^J \exp(\mu_{ik})}

where \mu_{ij} is the utility of the alternative ranked j-th for individual i. Here that is in code; the key bit is the rank_logit_lpmf function.

Model 3: Rank-Ordered Logit for Connected Graphs

ranked <- "// Again, note that this code is just reproducing https://khakieconomics.github.io/2018/12/27/Ranked-random-coefficients-logit.html
functions {
  real rank_logit_lpmf(int[] rank_order, vector delta) {
    // We reorder the raw utilities so that the first rank is first, second rank second... 
    vector[rows(delta)] tmp = delta[rank_order];
    real out;
    // ... and sequentially take the log of the first element of the softmax applied to the remaining
    // unranked elements.
    for(i in 1:(rows(tmp) - 1)) {
      if(i == 1) {
        out = tmp[1] - log_sum_exp(tmp);
      } else {
        out += tmp[i] - log_sum_exp(tmp[i:]);
      }
    }
    // And return the log likelihood of observing that ranking
    return(out);
  }
}
data {
  int N; // number of rows
  int T; // number of inidvidual-choice sets/task combinations
  int I; // number of Individuals
  int P; // number of covariates that vary by choice
  int P2; // number of covariates that vary by individual
  int K; // number of choices
  
  int rank_order[N]; // The vector describing the index (within each task) of the first, second, third, ... choices. 
  // In R, this is order(-utility) within each task
  matrix[N, P] X; // choice attributes
  matrix[I, P2] X2; // individual attributes
  
  int task[T]; // index for tasks
  int task_individual[T]; // index for individual
  int start[T]; // the starting observation for each task
  int end[T]; // the ending observation for each task
}
parameters {
  vector[P] beta; // hypermeans of the part-worths
  matrix[P, P2] Gamma; // coefficient matrix on individual attributes
  vector<lower = 0>[P] tau; // diagonal of the part-worth covariance matrix
  matrix[I, P] z; // individual random effects (unscaled)
  cholesky_factor_corr[P] L_Omega; // the cholesky factor of the correlation matrix of tastes/part-worths
}
transformed parameters {
  // here we use the reparameterization discussed on slide 30
  matrix[I, P] beta_individual = rep_matrix(beta', I) + X2 * Gamma' + z * diag_pre_multiply(tau, L_Omega);
}
model {
  // priors on the parameters
  tau ~ normal(0, .5);
  beta ~ normal(0, 1);
  to_vector(z) ~ normal(0, 1);
  L_Omega ~ lkj_corr_cholesky(4);
  to_vector(Gamma) ~ normal(0, 1);
  
  // log probabilities of each choice in the dataset
  for(t in 1:T) {
    vector[K] utilities; // tmp vector holding the utilities for the task/individual combination
    // add utility from product attributes with individual part-worths/marginal utilities
    utilities = X[start[t]:end[t]]*beta_individual[task_individual[t]]';
    rank_order[start[t]:end[t]] ~ rank_logit(utilities);
  }
}"

Take a second to reason through the above code; we’re able to use all of the information a full ranking provides to us, without having to force any interrelationship between best and worst to exist that may or may not exist. Also, notice that there’s nothing here specifically for the “connected graph of choices” part of model- that’s all in the choice set generation, which we’ll do next.

Building the choice sets

Let’s try making the choice sets the new way, and then compare models again.

To keep things simple, we’ll generate our graph of options as a balanced incomplete block design. This is a quick way to get the connectedness property we want in a reasonably statistically efficient design. It’s worth noting though that the literature on designing maximally efficient discrete choice experiments is deep, and has some beautiful connections to the broader experimental design literature; here’s a footnote if you’re the sort of person who finds that exciting⁹.

Updated data generation for connected graph illustration

# Number of individuals
I <- 400
# Number of tasks per individual
Tasks <- 10
# Number of choices per task
J <- 3
# Dimension of covariate matrix
P <- 5
# Dimension of demographic matrix
P2 <- 6

# demographic matrix
W <- matrix(rnorm(I*P2), I, P2)
# Loading matrix
Gamma <- matrix(rnorm(P*P2), P, P2)

# Show W * t(Gamma) to make sure it looks right
W %*% t(Gamma)

                [,1]         [,2]         [,3]          [,4]         [,5]
  [1,] -1.5412624799 -1.492174492 -5.116226558  -1.347807230  0.001574089
  [2,] -3.7257266234 -2.723139959 -2.822841903   0.101675159  0.763523848
  [3,]  0.4730049846  0.910864917  4.083111863  -1.471631667 -1.279559170
  [4,] -1.1289564758 -0.842181449 -2.937264693  -0.437928873  1.117150515
  [5,] -4.1905191154  3.318393965 -0.599327838  -3.310317907 -3.217285041
  [6,] -2.6492353829 -1.033830233 -2.811919648  -2.527213567  0.016064346
  [7,]  1.4591711093  4.353062595  4.524177128   3.791394053 -3.055228317
  [8,] -0.5661501940 -3.742740573 -5.053927670  -3.149040498  1.821531612
  [9,] -0.0569210566 -2.512887096  1.150352221  -4.158098423  2.544617312
 [10,]  0.8174140064 -3.257770914 -6.679072420   5.217851104  2.543820586
 [11,] -3.3644771332 -1.449694472 -5.814100769   0.034075107 -1.421669636
 [12,] -0.6646885105 -1.994607998  0.229234553  -2.930599125  2.334724287
 [13,] -0.4554012271  6.725575915  3.974224000   1.624482334 -2.525514532
 [14,] -0.4290414679 -2.914172693 -6.812881630   3.427975008  1.301192440
 [15,]  2.3451100762 -1.868093119  3.661664471   1.256928318  1.815057453
 [16,] -0.7424357061  1.314333858  2.326296328  -1.651657324 -1.242121310
 [17,] -1.1355535184 -1.554960793  3.500814111   2.138937785  0.399665190
 [18,] -2.4670221123 -2.350556927 -0.608124647   3.760988423  0.919302772
 [19,]  0.6336296989  1.098046588  2.083267335  -1.385890290 -0.900592817
 [20,] -0.1783054530 -1.782880847  1.301147942  -9.646402726  1.613337908
 [21,]  2.8588590755 -1.054297886  1.907116980  -7.060972023  1.038812638
 [22,] -0.2039825209 -1.247507832  2.470660505   1.239479888  1.106626373
 [23,] -3.4746506631 -3.508828446 -5.449731641  -7.648194070  2.333166463
 [24,]  2.5724355418 -0.160439092  4.824422518   4.478013435  0.809276982
 [25,] -0.9755464590  2.706429159  2.085153879  -7.002689457 -1.163424409
 [26,]  0.2742467555  2.783637461  0.487778649   1.538890117 -2.805469328
 [27,] -2.6327365951  0.647359172  1.006432011   0.690032677 -1.621858539
 [28,]  1.8763278673 -1.272990544  2.078108412  -5.243199305  2.082184979
 [29,]  2.2188968660  0.240736613 -0.815478862  -2.744060820 -0.048770516
 [30,] -3.9628733890 -0.113607772  0.224711894  -5.198084608 -1.114163049
 [31,] -4.2641345894  1.996915047  0.099624723   1.331299406 -1.516121432
 [32,] -0.5524263714  1.311462389 -1.214233177   2.584592840 -1.630199534
 [33,]  1.9944600720 -0.752170898  2.436853931  -5.344579832 -0.083796489
 [34,] -0.2059555289 -0.851039534  2.687940642  -6.656340585  1.886319532
 [35,] -1.7705938256  1.706629894  1.336601439  -5.630928568 -0.541385096
 [36,] -0.5273080110 -2.018241480 -0.700311999   3.621702144  0.279030488
 [37,] -2.7185644156  0.360535474 -0.169927022  -5.380951726  0.191448735
 [38,] -3.8600831761  3.791903631  0.447887021  -1.294971709 -1.569381142
 [39,]  5.2576883153 -4.031141475  0.981260876  -0.389974915  2.434063355
 [40,]  0.2264332470 -2.668998647 -2.151615699  -4.955313403  1.697963807
 [41,] -1.1162185151 -3.255401499 -3.544450808   1.238311086  2.634532420
 [42,] -1.4336601196 -2.100986861 -1.829304497  -6.780311417  2.468169302
 [43,] -3.2243263508  1.324538045 -0.541223217   3.753983082 -0.617814088
 [44,] -1.3758730861  1.959017116  5.831226366  -2.086105419 -0.270401213
 [45,] -1.9743016454  1.054005527  1.528013201  -4.542528278  1.105393765
 [46,]  3.0898340361  3.222170045  6.285250416  -0.762194443 -1.581173608
 [47,] -1.4055645007  2.087779424 -1.225676172   1.160761651 -0.050036984
 [48,]  0.4646083115 -0.981080358 -2.013141413  -1.061568791  0.095956075
 [49,]  0.6117466464  3.065568660 -1.186154449   1.919869983 -1.166239243
 [50,] -2.2382817344 -1.861037151 -1.984131835  -6.223047509  2.271733572
 [51,] -1.6750183655  0.314000156 -2.193940420   1.445680792  0.104323168
 [52,] -0.5424056953  4.635270387  5.493769693  -2.790746468 -0.647672683
 [53,] -1.0869464214  0.772149226 -0.457213500  -3.837057494  0.152893366
 [54,]  0.8314561763 -0.195856405 -2.802297421   9.629050659 -0.943798566
 [55,]  1.3873534192 -0.163905851  0.876675905   1.955158280  0.082625241
 [56,] -1.0028628371 -1.065792177 -0.387355992   2.934993534 -0.701583668
 [57,] -1.7069408645  1.963128354  2.492700571   0.324471935 -0.483430577
 [58,] -2.9952769643 -0.713819727 -7.860000192   5.633282644 -0.744285185
 [59,]  0.7948305572 -1.938533263 -1.219037953  -0.832381925  1.071975571
 [60,] -2.5497991144 -1.208014293 -4.972682949   1.382277683  1.349566702
 [61,]  1.2185172356  3.385639253  2.740549092  -0.303520107 -3.243991049
 [62,]  2.0346715861 -2.017907006 -1.661937529   1.712486110  0.269443568
 [63,]  1.1003817042 -0.069002060  4.935401334   0.300455176  2.408042531
 [64,] -0.6564698814 -0.826226916 -1.970112617  -6.712398580  0.958957912
 [65,]  0.3144406730  0.767761174  0.459991335   3.879362644  0.869803554
 [66,]  2.0929055159  0.854741777  5.230006519  -8.064751996 -0.145558505
 [67,]  0.7152519594  1.385720251 -2.925501025  -2.113992995 -1.654562613
 [68,] -1.0074155568 -0.430985382 -3.697355886   1.944809497  0.805377438
 [69,] -4.5877847402 -1.115382588 -6.601460733   6.511372347 -1.834392685
 [70,]  1.3149305093  1.375782792  2.276168002   4.572992536 -0.390324479
 [71,]  0.8831744784 -0.026936584  0.556239946   2.553588961  0.063788647
 [72,] -1.7266497964  0.112322933 -1.308367798  -5.016971355  0.921644440
 [73,]  1.7179999725 -1.072504568 -5.036695338   7.313133472 -0.605940628
 [74,] -2.5728198961 -0.860120812 -0.060575910  -0.999756268  1.768373949
 [75,] -0.5476467617 -1.159374665 -2.006579237  -2.268839513  0.187039974
 [76,]  0.1650960090 -3.708035880  0.102718618 -10.606112375  1.488022737
 [77,] -0.9721551093  2.158634003  3.189175876  -2.438725094 -1.106092489
 [78,]  0.4057151033  4.001283801  3.280511664  -6.940937119 -2.809459727
 [79,]  2.9478267878  1.135935148  2.142583188   0.664822137 -0.144571229
 [80,]  0.2755206984  0.472526312 -2.351192378  -1.849787822 -0.947808873
 [81,] -0.0248319352  0.814213514 -1.233446682  -5.156781210  1.299755378
 [82,]  0.1432628546  0.606484309  1.641764869  -2.981491335 -0.035143959
 [83,] -0.1959375904  4.032697198  4.292382205  -6.784462132 -1.637216987
 [84,] -1.0461641405  2.159511058  0.445951483   1.447644970 -1.506760901
 [85,]  0.3722287881  1.408425046  0.797829118  -3.724560916 -0.010446977
 [86,] -1.5351583649  2.038973232  1.638698528  -5.710825425 -0.304028004
 [87,]  1.7345228902  0.773871492 -2.567213936   0.384591360  0.002782755
 [88,] -1.2409675750  3.304432311  1.213745581  -4.811381380 -2.814723220
 [89,]  0.4792349258 -4.812218155 -5.104655158   0.764346408  1.267067517
 [90,] -2.8027120518 -1.429076516 -4.645677056  -0.247897534  0.077751243
 [91,]  4.8231762505  3.742086153 10.088221297  -2.527660071 -0.175154337
 [92,] -1.3797727159 -1.073837477 -2.111056183  -2.062980915  1.448070356
 [93,]  0.7100897071 -2.246621479 -1.580015393   4.568521557  1.621422544
 [94,]  4.8942011463  2.309683041  8.396684919  -4.541263246 -0.124288457
 [95,]  1.2013912261  2.895327974  2.410685637   1.870889346 -0.326469559
 [96,] -0.8505517239 -0.625527476 -2.527230492   7.491415972  1.899747380
 [97,] -1.5709793562 -1.079210748 -1.881011222   3.249166424 -0.313077176
 [98,] -3.5116383882  1.217197562 -5.823922097   1.055318646 -1.699281316
 [99,] -0.1567694031  1.792952395  1.768021069  -1.417357333  0.685322649
[100,]  0.9652400055 -3.179355675 -2.179006231  -3.020443444  1.655407424
[101,]  0.4823556858 -1.682962326 -1.131194413   1.254709438  1.020400258
[102,] -2.3044516832  4.097745443  2.977431366  -0.150843752 -1.990983859
[103,] -1.0175255260 -0.427489190  0.712764242  -2.452938001  1.702683003
[104,] -0.2698443257 -1.716970057  0.521238247   4.764269458  0.585221974
[105,]  2.5899921548  1.154353679  0.592179144   0.480689307 -0.342398771
[106,] -2.6282040844  0.781529018 -3.151733589   0.404017379 -0.855167573
[107,] -3.7578394431  0.725436982 -2.082676660   1.446067621 -0.210110152
[108,] -3.6568247106 -1.318774812 -3.388348305   3.349438898  0.545780946
[109,]  1.4037370745 -0.904889966  1.023022088  -2.192989050  0.481397270
[110,] -1.0884221704  2.806843203 -0.847539160  -0.818310380 -1.482521731
[111,]  6.6603289114 -4.923114094  0.894412761  -1.518693251  3.008703122
[112,] -0.9104537828 -1.758171794 -3.431275531   5.141225445  0.647030477
[113,] -1.5535052837  1.161424565  0.760838076  -7.247328990  1.499210539
[114,]  1.0682567477 -1.963708505 -0.148556728  -7.212397905  0.426371232
[115,] -0.0880081758  0.091873626  2.853263490  -4.913293403  1.406393441
[116,] -0.2691936928 -3.725975958 -3.560845295   5.169427116  2.031946658
[117,]  2.3775023933  0.892673643  4.204472259  -2.035742715 -0.659865801
[118,]  4.6448924786 -2.527037309 -0.616799447  -5.128486182  3.076479420
[119,] -4.3226389583  4.480860926 -0.974356555   6.303175192 -1.996697530
[120,] -1.2542342150  5.354119929  1.799319635   6.999319448 -3.252113365
[121,] -0.9852929611  0.951391968  1.551200330  -1.176931049 -1.308358452
[122,]  2.2852650200 -1.820923818 -0.555728479  -7.215712948  3.217758249
[123,]  1.1961343876  2.528461493  2.538794751   3.999815765 -2.261056170
[124,]  3.9434028987  0.468410063  2.953459433  -1.344569521 -0.650084090
[125,]  3.7113597666 -4.288455188 -3.817740127  -3.594567122  3.179976415
[126,]  0.1099306967 -2.651007315 -3.608223708  -6.043665129  1.578558040
[127,]  1.0696418393 -0.729894560  2.463588344  -0.522842823  0.040370391
[128,] -0.6787163713 -1.801860564 -3.739704396   7.603190616  0.554145723
[129,] -1.2859450363 -3.953714915 -3.736954229  -4.708744127  3.096942958
[130,] -0.3756811068 -1.153890578 -2.946557245  -4.158037571  0.176978925
[131,] -1.2186553868  2.246970408  2.212810159   3.881813134 -1.738799509
[132,]  1.9836235785 -2.887767397 -0.721256825  -2.384730879  2.013760760
[133,]  2.0292746447  3.028565720  7.656624845  -0.305502572 -1.068319678
[134,] -0.9594714878 -1.493083785 -1.351516594  -2.302979749  1.525383755
[135,] -1.9792088825 -1.850492713 -1.214298773   3.982149692  0.244041239
[136,] -0.4742228524 -2.106332820 -2.540659370  -3.294491828  0.530532763
[137,]  0.0571285304  0.232013344  1.960116037  -0.397719999  0.734578272
[138,]  0.5310229010 -1.094103849 -0.264782256  -3.190251020  2.244983085
[139,] -2.1888499907 -2.635693944 -0.988012506   3.885343352  1.547730404
[140,]  0.6508365425 -4.204347775  0.832578124  -4.900228760  2.863470736
[141,]  0.7740411320  2.823505148  1.409851042  -7.475820472  0.151679308
[142,]  0.0268653240  3.115058538  3.038754355  -2.264443676 -2.216430120
[143,] -0.0058500487  0.729084642  2.290477993  -0.513549995 -0.740090261
[144,] -3.3244178992 -2.863492582 -2.361741730  -3.665838545  0.580048164
[145,] -0.7783620259  3.129133473  2.223952122   2.255136322 -1.426592834
[146,]  1.8875741423 -0.617481752 -0.746104767  -0.419698149  1.588773517
[147,] -0.6827788185 -1.523238910 -1.341994806  -4.136552330  2.362416106
[148,]  2.1507516389  3.030654883  5.159838750  -3.983748678 -0.550696980
[149,] -1.6748774124 -0.288494688 -0.714806839  -9.996319038  1.140863331
[150,]  5.1522519392  0.307630449  7.003830438  -0.909847210  2.197515606
[151,]  1.9385032415 -2.249163256 -1.144024683  -2.044815308  2.281543871
[152,] -2.7899597051 -0.088987475  1.351830693   0.260738582 -0.709741001
[153,] -2.3066595278 -1.056895459 -1.638617662   2.749698665  0.754547895
[154,]  0.9800091577 -2.544119957 -3.086461348   4.519232301  2.606182930
[155,] -3.8686894949 -1.620129695 -2.351941041   0.602648533  0.714014301
[156,]  2.5383634852 -5.379374126 -0.716504320  -3.974143744  3.927922450
[157,]  1.1053159646 -0.422707780  2.008146213  -3.850246343  0.338358739
[158,]  0.3264388388 -0.593516866 -3.075237461   0.719912097 -0.112876982
[159,]  1.6012790692  1.812533336  1.797755581   3.660860044 -0.423047308
[160,]  1.2922500314 -0.459464352  2.215192655   1.352833069  1.206678403
[161,] -0.7622101570 -2.522655980 -4.213362587   4.912426613  0.103644226
[162,]  4.5404230491  0.721701168  1.555109280   1.299873582  0.628581632
[163,] -0.9540353397 -0.195694636 -0.161224837  -4.122842818  1.180071013
[164,] -3.6417478377 -0.394231940 -6.264626108  -0.608377097 -1.390711784
[165,]  2.1091447064  1.353748920  4.461868678   0.273458003  0.120112304
[166,] -0.0290936061  1.144045508  4.092093006  -1.936552910 -1.197377036
[167,]  1.4860203983 -0.584293520 -0.004069713   5.621312976  0.620950890
[168,]  0.3290043984 -1.612743964 -0.544021563  -1.438185235  0.963252517
[169,] -5.1572797410  5.401016313 -1.559269433  -0.231835231 -3.984857171
[170,] -2.6157717958  1.433195442  0.527952453   0.181755243 -1.578262619
[171,] -1.8740753120 -0.248012267  0.212349396  -3.833235112  1.451500183
[172,]  1.8811258923 -2.366310147  2.547098882  -2.288771052  1.765463673
[173,] -4.1492921663  4.586286028  1.142542347  -0.665391466 -2.760908221
[174,]  2.4120628060 -0.917798806  0.975102447  -3.880958877  1.111384476
[175,]  1.9592358643 -0.167929048  6.404296348  -0.358888539  1.558850196
[176,] -2.0423960871  0.433842229 -0.578914293  -4.341287710 -0.347192922
[177,] -1.4543777885  1.146837670  1.464567863  -3.091032963 -0.515891125
[178,] -1.4380563715  0.019447794 -0.579070961  -1.792870353 -0.197552068
[179,] -0.1278726908 -1.292790720 -1.737233273   0.038747398  0.376501230
[180,] -0.7135982395 -0.889820042 -3.973535998   4.333053174  0.152650807
[181,] -2.8937211604  0.041303667 -5.145715862  -0.765775591 -0.347524612
[182,] -0.6430920314 -0.295679308  3.406405768  -2.663172863 -0.159770649
[183,]  6.0267348568 -2.463655549 -0.494349697  -3.139696659  1.895277663
[184,]  0.7688866439 -1.255548051  0.727220543   1.427414584  1.339018973
[185,]  2.0103475737 -2.522123475  1.720728104  -3.390667013  2.680437375
[186,]  4.4206082338 -2.940748918 -0.949502862  -4.621107230  3.461719747
[187,]  5.4624956592 -3.767211808  0.194013035  -5.080256795  3.944792150
[188,] -0.2706198245 -1.174388056  0.220000788  -3.580691316  1.836674882
[189,]  0.6574342287 -0.266791033 -2.338009015  -1.313941916  1.072824906
[190,]  1.7708181838 -0.807595730  0.448269453  -0.264439364  2.080144818
[191,]  2.9277332197 -2.290685185  1.130965769  -2.648946888  2.294051234
[192,]  2.3376289738 -4.810455537 -5.743077235   2.245823198  1.551264919
[193,] -0.0445550487  2.050968440 -0.788438235  -2.000760923  0.049066720
[194,] -0.2721325117  3.543700004  3.880504277   0.071514588 -0.860674416
[195,] -2.8633225359  3.517464994  4.208164397  -2.116180407 -3.263312850
[196,]  5.3891735367 -3.600419049  2.771758646  -1.447716125  3.950849184
[197,] -2.2153679890 -2.406417057 -5.966472124   1.801221269 -0.122900748
[198,] -1.7246958740  0.743113811  3.653829564   1.135371579  0.820246965
[199,] -2.8245208248 -2.772037686 -2.803603884  -1.638859318  1.880857446
[200,] -2.1115753585  1.117852869  1.207609266  -0.182881622 -1.077794168
[201,] -0.0771998669 -4.656055217 -3.720279083  -0.971893333  2.461171395
[202,] -1.6639499127  3.480818549  3.186987522  -8.754646778 -0.317926593
[203,] -0.9280184504  1.506599081  0.646954668   1.382803939 -1.938449705
[204,] -0.4752444829 -1.275334331 -1.406909561   0.130797029 -0.995125045
[205,] -2.5620703689 -1.459616772 -4.390515016   3.062913186 -0.510308217
[206,] -1.6964887349  2.837466505 -0.641594375   2.783606274 -2.779898541
[207,] -4.6243368303  4.618496157  3.620502553   0.183556444 -2.999291378
[208,] -1.4746760840 -2.057579226 -4.030329668  -3.336698020 -1.132639487
[209,] -5.9629354032 -0.191445220 -4.329583840  -2.596478948 -0.689677710
[210,] -0.4764465345  1.745071235  1.635180190   1.322559406 -1.551122539
[211,]  0.7858118471  2.375928591  6.473000633  -5.145734534 -0.705885438
[212,] -2.0538425360  0.091292768  1.367583540  -4.274493629  1.190022621
[213,] -1.3686152481  4.686531589 -1.122755018   3.681627519 -3.869159447
[214,]  3.2086560459 -1.258015992 -0.460510159  -9.164499174  2.790891399
[215,] -3.7943660819  3.984004385 -2.011043112   0.804491187 -3.493218115
[216,] -1.6678732222  5.675634450  8.287417610  -5.213123717 -2.593581338
[217,] -1.1411340005  1.793941542 -0.732179612  -1.131659801 -2.240653511
[218,]  1.5877035866  1.407216793 -0.227111096   4.064433118 -0.970644315
[219,]  0.0213908464 -0.731075672 -1.702992361  -0.143091758  0.555208459
[220,]  1.2886190926  0.057318999  1.159205088 -11.900004402  1.339579101
[221,] -1.8874029196 -0.678092548 -1.922415493   2.536637609  0.176763838
[222,] -1.5481535790 -2.233664448 -3.931024438  -0.637148938  0.832523101
[223,] -3.1411791971  0.814817837 -1.791176437   0.067955883  0.770568438
[224,] -1.9656430756  3.230815079 -0.304634323   6.015425431 -1.824901702
[225,] -0.8718787657 -0.818840179  1.162383931  -0.708362934  1.397718748
[226,]  1.3325762625  1.445431300  2.659789287  -3.852118743 -2.155659613
[227,]  3.1841480716  2.102119921 -0.529331143   7.033242831 -2.707838078
[228,]  0.5196335301 -1.898114180 -2.264547089  -0.282825994  1.360109767
[229,] -1.6695109399  2.172579151  5.105986681   1.632536986  0.062691160
[230,] -2.1632159146  1.310763498 -2.434157235   1.286787544 -3.209746369
[231,]  0.0752715667 -0.240473372  0.164362746   0.331146101  0.134191803
[232,] -1.5988099760  5.448585338  4.785373130   0.368892045 -2.651398295
[233,]  2.4213864542 -2.152410691  4.724474115  -2.370211890  0.819671695
[234,]  1.2137433694 -4.215403271 -3.103534709  -9.346053778  3.678107036
[235,] -2.8397972037  3.825611752  1.307163761  -0.075750946 -3.080975910
[236,] -0.3262216246 -1.095678089 -2.925265561  -5.651569659  1.357939949
[237,]  1.2345872009  0.999166920  3.338740210   3.001784768 -0.994845288
[238,]  3.3868373207 -1.073254036  1.160564541  -0.373667439  2.246072143
[239,] -0.0004368311  1.984666603  2.678258421   0.044847876  0.317588987
[240,]  1.3609634388  2.096965468  1.260048693  -5.888458446  0.389243378
[241,]  3.2411279233 -1.891805623 -4.214081801   4.416902066  2.602978454
[242,] -2.7945577172 -1.104266271 -3.124331307  -6.694818023  2.103459667
[243,]  3.4899558567 -3.158850165  2.025287374   1.190410438  2.354619764
[244,]  2.8685019199  0.147372899  0.052249406  -1.686380913  1.520371056
[245,]  0.4407005088 -0.309341024  1.968090398   1.754090422  0.840552749
[246,]  0.4770725161 -2.379991071 -2.457669900   5.610080881  1.245976959
[247,] -1.3774356510 -2.182405551  0.136256314  -5.308583575  3.050991568
[248,] -2.3317120002  2.972877054 -0.592511418   0.757276508 -2.365377220
[249,] -1.2109381551 -0.268419650 -1.964015959   4.592277453  1.153552090
[250,]  8.2995536627 -4.848705764  3.341974461  -6.494704793  5.911996483
[251,]  4.4519126307 -1.222155245  1.973744988  -1.795958206  1.832629276
[252,] -2.8358609484  4.035047944  3.197290382   5.889419669 -2.793148700
[253,]  0.1735477201 -2.890729833 -4.033903460   0.763938445  1.679779840
[254,]  0.9916637745 -2.363646665  0.308022000   9.217571147  0.552854863
[255,]  0.1286702704 -1.753017494  2.990927854  -3.916214290  2.252696842
[256,] -3.0122367133  2.495467168  1.927867480   0.142405776 -1.786866625
[257,] -0.9432058286  0.143417466 -2.563617567   6.178052217  0.271914242
[258,]  1.0882938541  1.351837340  4.098870234  -2.012676795 -0.226826839
[259,]  2.5860141711 -0.289256836 -0.914675391  -1.664340430  0.556265495
[260,] -2.4457014342  2.178972355  2.153378187  -0.223719649 -2.402919174
[261,]  0.1457771480 -0.164694728  2.735616397   1.408696933  1.257969428
[262,] -4.7358857314  0.781087755  0.120120181  -0.730628058 -0.119145333
[263,]  1.8619444055 -1.678991032  2.643384992  -9.590256589  3.013505861
[264,]  1.1678683277  1.383399284  1.107692624  -3.070776033 -0.200076369
[265,] -2.0882424395  0.894465085 -0.915395936  -0.570762035 -0.015445036
[266,] -2.1141622513 -1.710927293 -0.077442988  -1.729499986  0.451725307
[267,]  0.3536386083  0.010927679  0.046937123  -3.813459412  0.933169264
[268,] -4.0197262095  2.369429440  3.584370545  -1.847007882 -0.871994427
[269,]  1.0997328660  7.747151983  8.261930139   2.629590525 -7.128945152
[270,] -2.2497716033 -1.121489710  0.060012025   6.045267586  0.699975151
[271,] -2.4918971173  3.872210353 -5.399106105   3.716266010 -4.443892725
[272,]  2.2895407874  2.331136532  6.253221578  -1.765980613 -1.448486227
[273,]  0.3567716687  1.533750574  1.325156379  -0.687018270 -0.554505501
[274,]  4.1093340293  1.144033063  3.077562557  -1.010487732  0.047011582
[275,]  0.7199603249 -0.671760402 -1.643340138  -1.805501803 -0.725343889
[276,] -1.1241502576  0.305943149 -0.579438747  -0.327705509 -0.312556120
[277,] -2.6235020669 -0.753496612 -3.889567627  -0.111757330  2.089658643
[278,] -3.3197459793 -2.312400949 -4.879802256  -1.007380467  1.064915835
[279,]  1.2467640802 -0.243609493  3.831151783  -7.446608384  2.017421902
[280,] -0.4810261503  2.628295900  2.882606665   0.577265948 -0.747689286
[281,]  3.4355030503  0.418218687  4.551389076   0.032214044  1.030933930
[282,]  1.5587065220 -4.354283522 -3.047823906  -2.085625730  1.985805211
[283,] -2.5532753739 -0.420570692  0.950776776  -3.315261795  0.398531703
[284,] -1.9458551434 -0.516262307 -1.760888750   7.229026143 -0.977479876
[285,] -1.9554143876 -2.000126858 -1.380702044  -1.325813108 -0.164382971
[286,]  0.3226550343 -2.561782501 -1.133588760  -0.321319492  2.262324510
[287,]  0.0502264452  4.046519646  0.817080464   8.984927525 -2.180837831
[288,] -1.7553511792  1.636438468  1.100463731   8.768995160  0.516296192
[289,] -0.0977606392 -2.793190676 -3.503937226  -2.445680922  0.710205286
[290,]  0.8037679703 -1.893508195 -3.339620655   4.013768293  1.187895707
[291,] -0.1548909875  1.760769330  1.282157811   0.645293941 -1.227291443
[292,]  2.9487681541  1.273837592  4.234409483  -3.464707154  1.703767268
[293,]  0.9599507814  6.051816624  4.268043832   5.554173771 -2.756390139
[294,] -2.1215692029  2.088479660 -0.091410396  -1.373266919 -0.537585980
[295,] -4.7531221693 -1.676852985 -4.542555984   1.773380186  0.433548877
[296,]  0.5419129330 -0.129916339  1.995913942   7.123535386 -0.074064828
[297,]  3.6950218515 -0.870476496  3.243938886  -4.182923198  2.795772909
[298,] -2.0809597737  0.208690990  2.273319084   1.423882574 -1.687377516
[299,]  1.4221278863 -0.501367710  0.035961812   4.122425866 -1.483101475
[300,] -0.3417691134 -2.116609045 -0.737975669   0.939547334  2.361158839
[301,]  1.3086492138  2.829847392  5.129402267   1.980408398 -2.009114492
[302,]  2.7637510247 -5.614449850 -6.235308012  -2.274170290  3.479506567
[303,] -2.0067482475  0.506826436  0.994977761   1.176581400  0.407421671
[304,] -2.3426833419  1.204946497  1.976441339  -9.014257342 -1.165731571
[305,] -1.7754009433  0.264859973 -2.897006697   0.267096243  0.229904567
[306,]  1.7069107567  0.263833742  0.731837028  -0.821499862 -1.201779921
[307,]  4.3459838333 -2.767812279  3.651435502  -0.558623218  1.882426894
[308,] -1.2561599715 -3.399856498 -6.893266020   2.980103189  0.458123764
[309,] -2.5895465056  6.232426601  2.036836456  -0.152024117 -3.451530867
[310,]  2.4490355263 -0.299455921  0.873034206  -3.625072865  1.641336088
[311,] -1.0024483978 -3.089324051 -4.911194845  -0.208256858  1.722434207
[312,]  0.3277856250  6.619762615 10.027982385   0.547909089 -5.352563263
[313,] -0.4138086439 -2.909200437 -6.022041961   2.853755441  1.568581140
[314,]  3.8143620300 -1.138003383  1.060894096  -3.343694692  2.986428684
[315,]  1.0389468466 -4.139732399 -2.299931804  -1.936714581  1.759737773
[316,] -2.8010586427  0.855179408 -3.937889107   3.922723463 -2.689481377
[317,]  6.0068667200 -2.452052067  0.926575212   4.950704980  1.312297273
[318,] -0.1843325608  1.364126137 -1.814461417   2.806065020 -0.993909431
[319,] -1.3714045030  0.556409910 -1.678003448  -1.641227464 -0.150295198
[320,]  0.3734495482 -3.813660303 -4.297412575  -2.989620528  2.972251935
[321,] -3.8111874424 -1.103618011 -5.718179987   0.038854005 -1.518296751
[322,] -0.4503106706 -1.593864412 -2.800863481  -1.191181804  2.734001096
[323,] -1.1646377032  1.377078821  3.748657213  -1.252972228 -0.567072995
[324,] -1.9196671897  0.240295466  1.163759700   3.036373295 -0.589744496
[325,] -2.4786493693  2.440088528  0.722633034   0.007714635 -1.832410270
[326,]  0.2759757669  1.803648954  4.229651552  -2.757308911 -0.119949920
[327,] -1.4406946784  0.280358344 -0.668566278   3.817194167 -1.048930730
[328,] -3.7195307994  0.132371450  0.048227163   6.582063532 -2.608982649
[329,] -0.0198730560  0.403507312  2.032030301  -2.353155918 -0.661874066
[330,]  2.6522478012 -0.875752249 -2.118692911  -8.584447666  0.676784210
[331,]  1.7482016944  3.463659667  5.083062612   0.407479241 -0.477145735
[332,]  3.2057224200  1.805907691  5.097592442   0.313955928 -0.194772102
[333,]  1.0830141534 -1.219515926 -0.606352500  -4.938554433  0.662739471
[334,]  2.1939255089  1.378546836  5.532263744  -4.825718910 -0.512513147
[335,]  0.2895851813  2.221719876  2.759203968  -1.138951352 -1.466324285
[336,] -1.6680114621  7.255253549  6.021719464   4.366479080 -4.849974098
[337,]  2.9149340641  1.813955090  3.099250059   2.831824079 -1.697545040
[338,] -1.3055925344 -3.543215074 -2.028013686   6.229299571 -0.004534920
[339,]  1.1292342049 -0.458117182  2.101925922  -1.413211931  1.020912005
[340,]  3.6242052353 -2.936378240 -4.912682226  -4.466661455  3.117054237
[341,] -0.6032583870  1.083878524  1.386151328  -6.760306755  0.273400085
[342,] -5.8408875246  0.364345708 -1.606113345  -3.941886283  0.297201502
[343,]  2.3528814653 -3.199928006 -0.156254785  -4.858713803  1.871573879
[344,]  2.9689024190 -1.871368564  0.907082987   1.068395241  1.253374176
[345,] -0.1994619398 -3.723697419 -3.951286919  -4.694055033  3.340486308
[346,] -1.1824969542 -1.223344393 -1.199173883   3.133309888 -0.561366430
[347,]  2.1900337315 -0.135932924  0.934022455  -1.823322418  0.871435934
[348,]  0.4574408621 -2.406204410 -0.073214145  -3.376441453  2.788976112
[349,] -1.9836893780 -4.130532452 -4.658974990   0.449849836  2.733680364
[350,] -1.1000589540  1.046733721 -1.399941834   8.164718767  0.014890627
[351,] -2.6052711423  4.492324992  4.856245845   1.956843650 -2.807793189
[352,] -2.0544991100  2.604083874 -2.878668632   3.878712133 -1.720694879
[353,]  4.0210815172 -1.419503975  6.110612202  -3.640316119  1.849353830
[354,]  2.4241853441 -0.508739794  2.508038096  -1.007456478  1.885966852
[355,]  0.2533916759 -0.509847528 -1.244787862   0.726972540  0.399856659
[356,] -4.3325638910  3.046534506  1.852349735  -7.777479756 -1.169035608
[357,] -3.2703336649  0.098420032 -2.445248917  -2.368231253 -0.377835768
[358,] -1.0002713655  4.085988836  0.689416392   1.569156500 -2.370921762
[359,] -2.7068700043  3.215339037  0.827195295  -1.834827563 -2.057384924
[360,]  1.9138198953 -0.156235204  3.225552439  -1.257812321  1.036676835
[361,]  2.0649131951  1.603782589  5.928903080  -1.478982412  0.402097546
[362,]  1.5732949929 -0.578526750  2.543265185  -0.820911660  0.355650098
[363,]  5.5018308327 -1.581824249  2.939668643   0.528774125  1.813284599
[364,] -1.8178895237 -0.837440767 -3.272859663  -0.524138343  0.163560944
[365,]  0.6379840448  1.488900809  1.950844328   0.150114299 -1.592521509
[366,]  5.1579950812 -3.064644101 -0.890625098   1.718199334  2.732690106
[367,] -4.9957252204 -0.782929898 -1.679901358   1.702519270 -0.206548879
[368,] -0.9901613421 -1.058866068 -2.217511692  -4.797444928  1.079949938
[369,]  1.3527372504 -6.651716236 -3.924733157   5.478559947  4.597974468
[370,]  0.3439111873 -0.114700380  0.475479632   2.040013425  0.165362378
[371,] -2.0698741776 -0.002093949 -2.623280008   0.402484692 -1.260062369
[372,] -3.6942055001  2.230241853 -0.745881858   0.408072506 -1.511184334
[373,]  2.6792189684  2.024161668  0.311133695   2.837114741 -1.228428290
[374,] -4.1200290084 -3.954575226 -6.182715237   4.540783902  1.141889815
[375,] -3.0389910195  2.127803778  1.473365199   4.749093504 -2.376674140
[376,]  2.9276471689  0.266925494  0.904247090   1.941260499 -0.437891765
[377,] -0.0996330446  3.480036806  1.227569879   2.408211004 -2.432264291
[378,]  5.6891190090 -2.271058567  1.273417213   1.822326223  0.172380537
[379,]  1.4151926928  4.432804917  2.225055529   0.604446965 -1.661457302
[380,]  4.6351507482 -2.861036630  0.581493898  -2.762629132  2.931101124
[381,] -2.0845835826  0.437190129 -4.311451776   1.695410225 -3.460968084
[382,]  0.5260692962  1.988186521  1.693206474   0.095309877 -0.654517664
[383,]  1.6234055938  0.249275099  2.935306962  -2.413443091  0.441335925
[384,] -0.0125945077  0.328927190 -1.000813898  -2.438906545 -0.335985361
[385,] -2.2971449989 -0.665226781 -2.276579138  -2.175669722 -1.186653499
[386,]  0.6265146542  4.274682751  2.014068285   3.072674562 -3.599539185
[387,]  1.2231314642  1.026097772  1.410984317   0.100134558 -0.538417996
[388,] -1.4290822752 -0.049352590  1.606991285  -1.730864779  0.473596109
[389,] -0.9372528372  1.615507190  1.480250427   5.063850969 -0.793673728
[390,]  1.8775016231  2.664194118  4.333083432  -0.301266972 -1.478078120
[391,] -2.8079047800 -0.583140754 -3.070339763  -4.734645850 -0.018427698
[392,] -1.3133931201 -4.057036154 -4.347763891   1.440233916  1.382953551
[393,] -1.5239131507 -0.239784368 -4.013799937   1.777625263 -1.034984876
[394,]  2.5665471713  1.542321872  1.205232466   3.159327011 -2.264593305
[395,]  3.3713851129 -3.331020990  0.488637965  -0.333368183  2.975326211
[396,]  3.8354075022  0.392123374  0.404967955   0.357983871 -0.791072036
[397,] -1.0800485176 -2.179510817  1.018442906  -3.379460578  1.394947831
[398,]  3.2589579784 -1.477800490  3.221959360  -2.584628239  0.968719704
[399,] -1.7802306370  3.391903955 -0.259946656   0.034371577 -2.626135392
[400,] -2.0346276380  0.687385242  0.123381126  -4.826722847 -0.954933517

# Correlation of decisionmaker random slopes
Omega <- cor(matrix(rnorm(P*(P+2)), P+2, P))

# Scale of decisionmaker random slopes
tau <- abs(rnorm(P, 0, .5))

# Covariance matrix of decisionmaker random slopes
Sigma <- diag(tau) %*% Omega %*% diag(tau)

# Centers of random slopes
beta <- rnorm(P)

# Individual slopes
beta_i <- MASS::mvrnorm(I, beta, Sigma) + W %*% t(Gamma)

# Again, quick plot to sanity check
plot(as.data.frame(beta_i))

bibd <- find.BIB(trt = Tasks, b = I*Tasks, k = J)

# The below function is not guaranteed to produce a valid BIBD; it may produce
# a close but imperfect design, which is sufficient for our purposes. I will
# however confirm it's connected, since we'd start to lose quite a lot of efficiency
# if it were not. 
# Returns 1 if the design is connected
is.connected(bibd)

[1] 1

# Create a mapping from BIBD options to attribute levels; idea is to make
# identically shaped inputs to our previous test, but using the conditions
# from the BIBD instead of sampling randomly for X.
option_attributes <- matrix(sample(0:1, 10*P, replace = TRUE), nrow = 10, ncol = P)

# Function to convert BIBD task to attribute matrix
bibd_to_attributes <- function(bibd_row) {
  matrix(c(
    option_attributes[bibd_row[1], ],
    option_attributes[bibd_row[2], ],
    option_attributes[bibd_row[3], ]
  ), nrow = 3, byrow = TRUE)
}

# Generate X matrix based on BIBD
X <- do.call(rbind, lapply(1:nrow(bibd), function(i) bibd_to_attributes(bibd[i,])))

indexes <- crossing(individual = 1:I, task = 1:Tasks, option = 1:J) %>% 
  mutate(row = 1:n())

# Write a Gumbel random number generator using inverse CDF trick
rgumbel <- function(n, mu = 0, beta = 1) mu - beta * log(-log(runif(n)))
mean(rgumbel(1e6))

[1] 0.5771369

# Ok, now we need to simulate choices. Each person in each task compares each 
# choice according to X*beta_i + epsilon, where epsilon is gumbel distributed. 
# They return their rankings. 

ranked_options <- indexes %>% 
  group_by(individual, task) %>% 
  mutate(fixed_utility = as.numeric(X[row,] %*% as.numeric(beta_i[first(individual),])),
         plus_gumbel_error = fixed_utility + rgumbel(n()),
         rank = rank(-plus_gumbel_error),
         # We're going to use the order rather than the rank in the Stan part of the model
         order = order(rank),
         # And here we create a dummy vector for the best choice
         best_choice = as.numeric(1:n() == which.max(plus_gumbel_error)),
         worst_choice = as.numeric(1:n() == which.min(plus_gumbel_error))
  )


tt <- ranked_options %>% 
  group_by(individual, task) %>%
  summarise(start = min(row), 
            end = max(row)) %>% 
  ungroup %>%
  mutate(task_number = 1:n())

`summarise()` has grouped output by 'individual'. You can override using the
`.groups` argument.

Now we have simulated data that represents a connected graph design with 10 total choices, and 15 choice sets per respondent. Each choice set contains 3 options, and each pair of choices appears together in at least one set across all respondents.

Evaluating our new model

So now we’ve introduced the ranked choice logit model, and simulated data that fits the connected graph structure. How do the three models perform?

Again, some data plumbing in the drop down below, so feel free to jump to the results plots after it.

Comparing ROL to best choice, and MaxDiff models

data_list_best_choice <-list(N = nrow(X),
                             T = nrow(tt),
                             I = I, 
                             P = P, 
                             P2 = P2, 
                             K = J, 
                             # NOTE!! This is the tricky bit -- we use the order of the ranks (within task)
                             # Not the raw rank orderings. This is how we get the likelihood evaluation to be pretty quick
                             choice = ranked_options$best_choice,
                             X = X, 
                             X2 = W, 
                             task = tt$task_number, 
                             task_individual = tt$individual,
                             start = tt$start, 
                             end = tt$end)

data_list_best_worst <-list(N = nrow(X),
                            T = nrow(tt),
                            I = I, 
                            P = P, 
                            P2 = P2, 
                            K = J,
                            choice = ranked_options$best_choice,
                            worst_choice = ranked_options$worst_choice,
                            X = X, 
                            X2 = W, 
                            task = tt$task_number, 
                            task_individual = tt$individual,
                            start = tt$start, 
                            end = tt$end)

data_list_ranked_rcl <- list(N = nrow(X),
                             T = nrow(tt),
                             I = I, 
                             P = P, 
                             P2 = P2, 
                             K = J, 
                             # NOTE!! This is the tricky bit -- we use the order of the ranks (within task)
                             # Not the raw rank orderings. This is how we get the likelihood evaluation to be pretty quick
                             rank_order = ranked_options$order,
                             X = X, 
                             X2 = W, 
                             task = tt$task_number, 
                             task_individual = tt$individual,
                             start = tt$start, 
                             end = tt$end)

compiled_best_choice_model <- stan_model(model_code = best)
best_choice_fit <- sampling(compiled_best_choice_model, 
                            data = data_list_best_choice, 
                            iter = 800)

Warning: There were 2 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

# Compile and fit MaxDiff model
compiled_maxdiff_model <- stan_model(model_code = best_worst)
maxdiff_fit <- sampling(compiled_maxdiff_model,
                        data = data_list_best_worst,
                        iter = 800)

Warning: The largest R-hat is NA, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

# Compile and fit ROL model
compiled_rol_model <- stan_model(model_code = ranked)

Warning in readLines(file.path(makevar_files)): incomplete final line found on
'C:/Users/timma/Documents/.R/Makevars.win'

rol_fit <- sampling(compiled_rol_model,
                    data = data_list_ranked_rcl,
                    iter = 800)

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: The largest R-hat is NA, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

best_choice <- as.data.frame(best_choice_fit, pars = "beta_individual") %>%
  gather(Parameter, Value) %>%
  group_by(Parameter) %>%
  summarise(median = median(Value),
            lower = quantile(Value, .05),
            upper = quantile(Value, .95)) %>%
  mutate(individual = str_extract(Parameter, "[0-9]+(?=,)") %>% parse_number,
         column = str_extract(Parameter, ",[0-9]{1,2}") %>% parse_number) %>%
  arrange(individual, column) %>%
  mutate(`True value` = as.numeric(t(beta_i)),
         Dataset = "Best Choice")

maxdiff_choice <- as.data.frame(maxdiff_fit, pars = "beta_individual") %>%
  gather(Parameter, Value) %>%
  group_by(Parameter) %>%
  summarise(median = median(Value),
            lower = quantile(Value, .05),
            upper = quantile(Value, .95)) %>%
  mutate(individual = str_extract(Parameter, "[0-9]+(?=,)") %>% parse_number,
         column = str_extract(Parameter, ",[0-9]{1,2}") %>% parse_number) %>%
  arrange(individual, column) %>%
  mutate(`True value` = as.numeric(t(beta_i)),
         Dataset = "MaxDiff")

rol_choice <- as.data.frame(rol_fit, pars = "beta_individual") %>%
  gather(Parameter, Value) %>%
  group_by(Parameter) %>%
  summarise(median = median(Value),
            lower = quantile(Value, .05),
            upper = quantile(Value, .95)) %>%
  mutate(individual = str_extract(Parameter, "[0-9]+(?=,)") %>% parse_number,
         column = str_extract(Parameter, ",[0-9]{1,2}") %>% parse_number) %>%
  arrange(individual, column) %>%
  mutate(`True value` = as.numeric(t(beta_i)),
         Dataset = "Rank-Ordered Logit")

# Combine all results
all_results <- bind_rows(best_choice, maxdiff_choice, rol_choice)

ggplot(all_results, aes(x = `True value`, y = median, color = Dataset)) +
  geom_linerange(aes(ymin = lower, ymax = upper), alpha = 0.3) +
  geom_point(aes(y = median), alpha = 0.5) +
  geom_abline(intercept = 0, slope = 1) +
  labs(y = "Utility Estimates",
       title = "Utility Estimates from the Three Models",
       subtitle = "With 90% CIs") +
  scale_color_manual(values = c("Best Choice" = "blue", "MaxDiff" = "red", "Rank-Ordered Logit" = "green")) +
  facet_wrap(~Dataset) +
  theme(legend.position="bottom")

At a squint, rank-ordered logit > best choice > MaxDiff as I’ve argued. Notice though that unlike the MaxDiff model, both best choice and rank-ordered logit look unbiased; the ROL model simply seems more efficient. Let’s check the RMSE to confirm:

all_results %>%
  group_by(Dataset) %>%
  summarize(RMSE = sqrt(mean((`True value` - median)^2)))

# A tibble: 3 × 2
  Dataset             RMSE
  <chr>              <dbl>
1 Best Choice        0.483
2 MaxDiff            0.556
3 Rank-Ordered Logit 0.417

The ordering is as we’d expect from the plots.

If I have free choice over the design of a choice experiment, this sort of connected graph design, analyzed with the rank-ordered logit is what I’d generally prefer to use. Designing the choice graphs to make this work well isn’t particularly hard with modern software, and the rank ordered logit model is widely available, not just in a sophisticated Stan model like I cribbed from Jim Savage here, but also in Stata, and other commonly used languages.

Hopefully the theoretical discussion and simulation study so far are persuasive: you can and should aim to do better than the MaxDiff model of discrete choice.

Alternative #2: a backup plan- rank-ordered logits with unknown middle preferences

There’s still an interesting question left here: what to do if you already ran an experiment designed without connectedness in mind, or your survey platform won’t easily allow you to make a connected design?

It’s not my preferred choice if I have control over design, but in this scenario we can use a Rank-Ordered Logit with unknown middle preferences.

Allowing unknown options in rank-ordered logit models

Let’s reuse the example above where we ranked my preferences over 2024 presidential candidates. The complete ordering for me is: \textcolor{blue}{\text{Biden}} > \textcolor{green}{\text{Stein}} > \textcolor{yellow}{\text{RFK}} > \textcolor{red}{\text{Trump}} If we collected data from me though by only asking my best and worse choices, we’d have: \textcolor{blue}{\text{Biden}} > ??? > ??? > \textcolor{red}{\text{Trump}} To keep using a rank-ordered model here, we’d need to figure out how to reason about the middle two options being unknown. Allison and Christakis, 1994 discuss various ways to handle this, but the one that seems most promising and that we’ll discuss operates as if the true ranking exists, but is unknown to the modeler. That seems pretty reasonable, since we’re not saying the middle preferences have to be strong, just that people likely aren’t totally ambivalent between the two in general. To handle this, we’ll first enumerate the possible ways the middle options could pan out. Here’s there’s just two: \textcolor{blue}{\text{Biden}} > \textcolor{green}{\text{Stein}} > \textcolor{yellow}{\text{RFK}} > \textcolor{red}{\text{Trump}} \vert \textcolor{blue}{\text{Biden}} > \textcolor{yellow}{\text{RFK}} > \textcolor{green}{\text{Stein}} > \textcolor{red}{\text{Trump}} Then, to construct a likelihood, we’ll marginalize over all the possible middle options. Since the two ways the middle options could shake out are mutually exclusive, that implies: P(\textcolor{blue}{B}>\textcolor{green}{S}>\textcolor{yellow}{R}>\textcolor{red}{T} \text{ or } \textcolor{blue}{B}>\textcolor{yellow}{R}>\textcolor{green}{S}>\textcolor{red}{T}) = P(\textcolor{blue}{B}>\textcolor{green}{S}>\textcolor{yellow}{R}>\textcolor{red}{T}) + P(\textcolor{blue}{B}>\textcolor{yellow}{R}>\textcolor{green}{S}>\textcolor{red}{T})

Now, let’s express this using the notation from Allison and Christakis (1994), which will help us handle the subscripts more easily. For a given individual i, let \mu_{ij} be the deterministic utility of the alternative ranked j-th. The probability of observing a particular complete ranking r (where r is one of the possible rankings consistent with the observed partial ranking) is:

P(r) = \prod_{j=1}^{J-1} \frac{\exp(\mu_{ij})}{\sum_{k=j}^J \exp(\mu_{ik})}

where J is the total number of alternatives. For our example with unknown middle preferences, the likelihood becomes:

L(\textcolor{blue}{B}>\textcolor{green}{S},\textcolor{yellow}{R}>\textcolor{red}{T}) = \frac{\exp(\mu_{i\textcolor{blue}{B}})}{\sum_{j \in {\textcolor{blue}{B},\textcolor{green}{S},\textcolor{yellow}{R},\textcolor{red}{T}}} \exp(\mu_{ij})} \cdot \left[ \frac{\exp(\mu_{i\textcolor{green}{S}})}{\exp(\mu_{i\textcolor{green}{S}}) + \exp(\mu_{i\textcolor{yellow}{R}}) + \exp(\mu_{i\textcolor{red}{T}})} + \frac{\exp(\mu_{i\textcolor{yellow}{R}})}{\exp(\mu_{i\textcolor{green}{S}}) + \exp(\mu_{i\textcolor{yellow}{R}}) + \exp(\mu_{i\textcolor{red}{T}})} \right] \cdot \frac{\exp(\mu_{i\textcolor{red}{T}})}{\exp(\mu_{i\textcolor{red}{T}})}

The central terms are the key to the unknown ranking model¹⁰. This lets us fit a ranked model and avoid MaxDiff, even with unknown middle options.

Of course, this approach isn’t without issues. Obviously, with a bit more foresight on design we could’ve learned more from the middle options while still asking respondents to make two choices per task. Putting that aside, the number of unknown middle options grows as J_{unknown}!, which makes calculating all the likelihood terms a real, real pain in the ass. If we re-do our initial example with 10 total alternatives, the middle 8 will be unknown, and there will be 40,320 permutations.

Building the ROL with unknown middle preferences in Stan

I’ll present two versions of the ROL with unknown ties model. First, the full marginal likelihood model described above, and then the Efron approximation you can use when the number of unknown middle options gets unwieldy. We’ll build towards comparing MaxDiff, ROL with known middle options, and the two ROL with unknown middle options all together.

To make it easier to follow logic for many middle options, we’ll use J = 5¹¹. The initial block below just sets up our DGP again with that J, so feel free to skip over it.

Setting up sim study with J = 5

# Number of individuals
I <- 400
# Number of tasks per individual
Tasks <- 10
# Number of choices per task
J <- 5
# Dimension of covariate matrix
P <- 5
# Dimension of demographic matrix
P2 <- 6

# demographic matrix
W <- matrix(rnorm(I*P2), I, P2)
# Loading matrix
Gamma <- matrix(rnorm(P*P2), P, P2)

# Show W * t(Gamma) to make sure it looks right
W %*% t(Gamma)

               [,1]        [,2]         [,3]         [,4]         [,5]
  [1,] -1.216579571  2.70873392   4.76010662 -4.382515769 -0.937908838
  [2,] -1.387844888  0.17885653  -1.33873062  2.035207807 -2.473143410
  [3,] -1.974087987  2.54974157  -0.05054604  0.473104518 -0.319766538
  [4,] -0.038572162 -0.12696491  -4.25409082  2.960893246  1.404494080
  [5,] -1.973286116  1.83777720   3.11764558 -1.396922278 -3.347534740
  [6,]  0.021811383  0.13485559  -0.55249459 -2.480132938  1.291705586
  [7,] -0.695194602 -2.79267967  -4.67942726 -5.346484029 -1.963929107
  [8,]  0.172634748  0.36643230  -2.25234689  1.168163720  1.368807721
  [9,] -0.077259388 -0.59355347  -3.31794394  5.211427489 -0.505453554
 [10,]  0.360707691 -2.45773636  -8.34096274  3.798035117  0.199426789
 [11,]  2.750266966 -1.82184235   0.48818279  4.393474270  1.757772220
 [12,] -0.679066320  1.31976351  -2.30763550  1.889675707  2.789800339
 [13,]  1.356207516 -0.77427314   6.61269431 -3.815191591 -2.175364307
 [14,]  0.766761499 -0.34906173   4.66706502  1.696090489 -2.510595159
 [15,] -0.996255136  0.67730930   0.16777477  1.597980539 -1.645200273
 [16,]  3.340993722 -2.72091034  -2.83588751  1.817107395  4.048968295
 [17,]  0.351522713 -1.02381470  -7.91429219 -0.228757153  4.374795113
 [18,]  1.240716479 -2.03866079  -5.71359396  4.494115105  0.583872144
 [19,]  0.058400326  0.04797676  -0.33386413  1.875720886 -0.491040485
 [20,]  0.977746395 -0.83279447  -3.17861331 -3.326265046  3.883902047
 [21,] -1.582000171  2.31009955   4.24136353 -1.551488120 -1.241483096
 [22,] -0.514563884  0.57884167   0.47157613  2.094924274 -0.644648908
 [23,] -2.406167342  3.98769875   1.32632682 -3.453997963  2.294923969
 [24,]  1.183554922 -2.37631996   2.49497064 -5.268382542 -2.512550749
 [25,]  0.189164826 -0.17422938   0.07025894  2.769579833 -1.050418189
 [26,] -1.060474995  1.25949950  -4.62147738  0.029562404  3.640131139
 [27,] -1.429353764  1.12960654  -0.20956355 -0.968738060 -1.075509890
 [28,]  2.468893040 -2.27091117  -2.04595658  5.960297725  1.714721598
 [29,] -0.918459168  0.59887699  -3.84139586  0.994688172  2.235764032
 [30,] -1.686497217  1.20746453  -5.77564749  0.768937050  1.798329665
 [31,] -1.188033561  0.26178445  -2.40960031 -7.464547340  1.120329318
 [32,]  1.862992252  0.44641916   1.03165331  7.273582356  1.929000017
 [33,] -1.127011252  0.91077583  -0.07401327  4.276614591 -0.781270603
 [34,] -0.181426645  2.52195817   7.55051619 -2.290981931 -0.329052100
 [35,]  2.085501801 -2.58577267   2.52995115 -4.481252700 -0.342047830
 [36,]  0.701709411  0.16019537   1.43646575  0.278467356  0.510498483
 [37,] -0.671479819  3.47314033   5.36198511 -1.884920932  2.215536088
 [38,] -2.454845069  1.12652917  -6.99112053 -0.392976205  0.478820446
 [39,]  1.730178777 -0.89022809  -1.05486418  1.620248890  3.278902366
 [40,] -1.663282979  1.39908216  -5.21125371 -2.110969958  3.436185997
 [41,]  1.797840016 -1.86535129   2.52832408 -4.538488314 -0.079499900
 [42,]  1.038719451  0.26498961   3.26465934  0.599748301 -0.033729624
 [43,]  2.297279750 -1.58429404  -0.23695673  1.194527815  2.776576843
 [44,] -0.071852717  0.98513223  -4.43777834  2.495348624  4.374783717
 [45,] -1.543892918  1.97199815  -0.16631939  0.487586063 -0.283565322
 [46,]  0.426131807  1.97736196  -0.12517925  4.346262148  3.219351178
 [47,]  2.265043041 -2.09171520  -2.00845806 -0.068751460  1.797022682
 [48,]  0.567383299 -0.29024877   0.29224939  3.676598240 -0.621641327
 [49,] -0.176345124  0.44107718   4.89232266  1.126347219 -3.471926737
 [50,] -0.844702953  0.82349247  -0.49221873  0.307763571 -0.623136310
 [51,] -2.053192939  1.97084844   1.01393733 -1.086254711 -1.565288540
 [52,]  0.229281338  1.41228892  -0.90349963  4.346329529  3.349979294
 [53,]  0.984497495 -1.03927697  -4.25912636  2.922199666  2.442073417
 [54,] -3.091970408 -0.37952713  -5.18425564 -5.550823338 -3.153001017
 [55,] -0.131057602  1.20946429   0.82598928 -0.378467734  1.800731826
 [56,] -0.213103075 -1.65817438  -7.01258138  2.978287225  1.032440847
 [57,]  0.061283927 -1.69156643  -4.14396688  0.070556213 -0.722425082
 [58,] -0.049374049  2.57088819   8.30599698  0.601156956 -1.283911885
 [59,]  0.230893788  1.62922104   0.47314110  1.782321613  2.990915322
 [60,] -0.512382496  2.97945064  10.02622119 -0.832924450 -1.240592177
 [61,] -0.366173333  0.11557829  -3.87356417 -3.233602444  2.639237665
 [62,]  0.869895687  0.76576412  -0.43681650  4.001248944  3.176463018
 [63,]  0.473061889 -1.33483784   0.09255220 -1.490319783 -1.070715599
 [64,] -2.356754497  3.15982518  -0.99853340  0.735565698  0.936753056
 [65,]  0.201285345  1.41296717   8.27357369 -1.729195660 -2.741729085
 [66,] -0.401052793  0.86582842   4.54618171 -0.609548824 -1.579299240
 [67,] -2.863339391  1.71691761  -2.78441312 -0.955181086 -2.655232845
 [68,]  0.495327367  2.11081603   2.97899529  1.879336443  2.974633301
 [69,] -1.850276454  3.08091326   7.85568998 -2.790804301 -3.194963747
 [70,]  2.423489585 -3.07931789  -0.42103632 -2.065518997  0.143903961
 [71,]  0.587496372 -1.70882465  -0.38448361 -2.391568619 -1.013952113
 [72,] -2.223867231  3.84043486  -0.42335340  1.339277799  2.881342856
 [73,] -0.393369848  1.08856200   0.04034497 -1.685582733  2.328821949
 [74,] -0.201767104  0.44338176   1.24175513  0.246185021 -0.464638520
 [75,] -0.125690238  1.59907552   4.11853652 -3.969937744  0.674956172
 [76,]  2.573570274 -2.34340084  -5.26258972  3.990105924  3.116787863
 [77,] -2.367190587  2.90678698   2.14895611 -2.924219309 -0.294068412
 [78,]  0.223591758 -1.98841403  -1.81488639  0.441586505 -2.120663842
 [79,]  1.702386507 -1.65439290  -2.24642333  1.899063074  1.869864722
 [80,]  0.623283342 -2.37705031  -8.44270273 -0.320001760  2.197975787
 [81,] -1.135596092  1.93515262   2.26114534  0.187201581  0.053034881
 [82,]  0.992150277 -0.24651493   2.72928530  4.394954888 -0.835769166
 [83,]  0.888353672 -1.08083208   2.31564466 -3.506438335 -0.689231950
 [84,] -0.459142320 -1.96891848  -4.27152444 -1.673405034 -2.078976051
 [85,]  0.387188000 -1.59618300  -3.06711152 -1.776670111 -1.013240084
 [86,]  0.885856928  0.72663880   1.62617603  2.411960084  2.159386219
 [87,] -0.684864533 -0.52311775  -1.20279391 -1.366036983 -0.844771471
 [88,] -1.401212210  1.84508195  -3.17911088  1.012386180  2.650271632
 [89,]  4.381356302 -4.38588832  -2.11858977  4.122712243  1.924108586
 [90,] -2.468107679  3.08092735   0.73055995 -0.696269187  0.248239755
 [91,]  1.544479345 -1.85119287   2.16967643 -1.169893642 -0.130018928
 [92,] -2.474136352  1.86369601   2.15308368 -5.887644964 -3.080379017
 [93,]  1.397847763 -4.57065704  -9.73244829  2.421842468 -0.196002338
 [94,]  0.277850377  1.89699566   4.13484494  1.285801323  1.516673935
 [95,]  1.250203436 -1.88929158   1.13290415 -0.002303534 -1.038764649
 [96,] -0.520735336 -0.99436594  -4.46646725 -3.325818707  1.530780110
 [97,] -0.324184690  1.46353273   8.31126950 -4.562419476 -1.593689494
 [98,]  2.122994966 -3.51104869  -2.83203901 -0.025147839 -0.602620051
 [99,] -0.909697998  0.91702620   5.66086627  4.051043931 -4.876377716
[100,]  6.274110792 -6.90420006   1.42984899  2.353994272 -1.258770227
[101,] -3.438792759  1.56280262  -2.66053873 -3.811375308 -3.053074780
[102,] -1.694283234  0.85737290   0.91128272 -5.916134525 -0.433233922
[103,]  0.767207857 -0.89209510  -0.67522547 -0.845311186  0.339411959
[104,] -2.019318045  0.57742565  -1.93601609 -6.291971534  0.780497784
[105,] -1.929664677  4.55633777   7.38181491 -1.545012882  0.598312053
[106,]  0.682596844 -1.36802393  -1.71570849  2.368412757 -0.824857464
[107,]  0.876771519  0.99417037   1.88084167  5.150760527  1.850837612
[108,] -0.478065674  0.06722172   0.72906151 -1.793858351 -0.998983931
[109,]  3.148440674 -4.31802122   0.61109511  2.957494584 -2.521667099
[110,]  2.332325771 -3.41243063  -2.76281048 -0.308891070  0.517807936
[111,]  0.336806962  0.74171587   4.00147266 -0.117554886  0.218758744
[112,]  0.814612027 -0.48055430   0.38107544  1.051405077  0.079344009
[113,]  0.871134226 -1.14226771  -4.50155972 -1.529495949  2.723019057
[114,]  0.137997442  1.27417976  -0.52469212  0.935740984  2.552960788
[115,] -0.648607764  0.11710088   0.02121035 -0.038607345 -1.216682510
[116,] -0.411028381 -0.90852190  -3.52641943  1.598940129 -0.516388614
[117,] -0.533756717  1.09532411  -0.88668492 -1.491040358  0.735452080
[118,] -1.980191898  1.88182486  -2.70821992 -0.018012907  0.951254003
[119,]  1.484721941 -2.72138813  -2.99899864  3.714298943 -1.517649519
[120,] -0.279801031  0.15004612   0.78295515 -2.799122401  0.529306273
[121,]  2.521325745 -1.35081784   2.96996259  6.224844877 -0.655644812
[122,] -1.480463796  0.50282676  -1.32456803 -1.239108858 -0.886213725
[123,] -1.624486450  1.94913686   2.37968703  2.420097545 -1.428893752
[124,] -1.527492091  0.97270068  -2.92029912  0.152333844 -0.004997683
[125,]  0.110434607  0.81278610   3.24195170 -1.064361658  1.410746492
[126,] -0.080000347  0.42826225   2.59413550 -0.459459342 -0.586792140
[127,]  0.547175594 -2.01856994  -2.89988873 -1.832675271 -0.825523035
[128,] -0.220517419  0.82742692  -3.81056790  3.203059592  3.124891282
[129,] -0.440891654  2.52287285   5.39249625  1.078463945 -0.135110820
[130,] -2.090212047  1.47588635  -1.14286569  0.364415762 -1.159887965
[131,]  1.850219452 -3.52137692  -3.48196465 -3.311741555 -0.621943492
[132,]  1.812685259 -3.66820240   0.13990832  0.178667514 -3.709772974
[133,]  2.149510053 -1.61372123  -1.36428389 -0.080034562  3.544327896
[134,] -3.512844529  1.92659074  -5.61335705 -2.810951778  0.590229368
[135,] -2.896442106 -0.16840446  -5.03397850 -4.010771531 -3.729821885
[136,] -1.462333934 -1.63124024 -10.53648443 -3.664277060  1.420516356
[137,]  0.123193012  0.62010426   2.22876695  5.535365497 -1.178103868
[138,] -0.447352794 -0.49550143  -2.34934094  0.043305837 -1.050683271
[139,]  0.777457961 -0.18369536   5.21940396  0.677642015 -2.921129496
[140,] -2.312645218  3.31992302   2.23428037 -0.430973044  0.627576312
[141,] -1.500505811  1.34191676  -1.39205461  4.276350043 -1.388592031
[142,] -1.644217751  0.89374531   1.85432683 -9.325877489  0.029215124
[143,]  2.492114502 -2.71347132   1.15043127  0.913041660 -0.202423836
[144,]  0.408491669 -0.45222392   2.49349757 -2.531882415 -0.583093674
[145,]  0.182001420 -0.09272287  -1.43622392  1.943633243  0.425441117
[146,] -0.200408115 -2.52057517  -2.66381884 -5.411873939 -2.069164138
[147,]  0.035500446 -0.74051771  -3.92860362 -1.140613560  1.611409944
[148,] -0.646451266  0.88598255   3.01446603 -0.261463123 -2.195667977
[149,] -0.604713208  0.28230107  -0.81855013 -1.644805266  0.613821271
[150,] -2.107599687  1.90104058   1.88421609 -2.502571513 -2.060945188
[151,] -0.445518203 -0.55285880   0.38645170 -4.375456386 -1.244496527
[152,]  0.611360026 -0.83432704   3.36008137 -0.192024893 -1.369100762
[153,]  2.187110654 -0.20694617   1.38949758  4.113188578  2.756647523
[154,] -0.966812387  2.28392503   4.66521209 -4.916234095  0.270545511
[155,]  3.352008963 -3.09523778   1.19336681  6.771299365 -0.196799459
[156,]  0.518269302 -0.94897508  -1.79771941 -3.346219475  1.478154285
[157,]  1.987866782 -1.24563138  -0.06903570  5.355418317  1.277002150
[158,]  0.646448796 -0.09152440  -1.87568735 -3.543165816  3.585933350
[159,] -1.753567835  0.70240524  -2.84150293 -4.589372444  0.232967254
[160,] -2.315575731  2.61993301  -1.49042411  2.540972803  0.118385287
[161,] -0.821325531  0.01187968   0.75625448 -0.987762785 -2.599761936
[162,] -1.014568008 -0.88657609  -3.68394711 -3.833012602  0.040662540
[163,]  1.551996580  0.08601909  -0.53501162 -1.717490878  3.400728365
[164,]  0.639023853  1.58916639   3.29352178 -0.001243981  1.568302699
[165,]  2.443656175 -0.66648288   2.82710839  2.960306749  1.512583562
[166,]  0.359712946  1.98484048   5.35210133 -1.663816196  1.982797708
[167,]  0.949191424 -1.47522732  -0.97214673  1.467734123 -0.659636087
[168,]  0.128632745  1.11502559   2.66918688  0.632348742  1.766940746
[169,] -1.800453294  3.35248183   3.17538074  3.272052453 -0.566755837
[170,] -1.308135066  0.33724790  -3.16685871 -1.774634812 -1.854095021
[171,] -0.723454249  1.02899463  -4.46146335  0.014200403  1.943952229
[172,] -2.063358646  1.42117869   0.18854613  3.471249906 -2.257953478
[173,] -1.301690233 -0.85373707  -5.09426395 -2.655800117 -1.550760916
[174,]  1.020703164 -1.36402370   3.18561817 -1.036293935 -2.420110210
[175,]  2.813693451 -4.28902809  -2.98174492  0.387300463 -0.728057270
[176,] -0.420182732  1.04435321   4.88886138  0.716124067 -1.448363331
[177,]  0.905492114 -3.25620231  -2.29877884  1.293200856 -3.298427077
[178,]  1.100650877 -0.64417128  -4.98579805  0.378166681  3.709196577
[179,] -1.137135392  2.25819965   3.61087018 -1.888399094  0.245940323
[180,] -2.280841595  1.41123787   1.74442659 -3.949310420 -2.305304559
[181,] -3.453812367  3.07655760   3.03993861 -6.202806321 -2.289204640
[182,]  2.196224048 -1.93258726  -1.75273338  0.832002819  2.840487887
[183,] -0.466221996  0.93065594  -0.06278044 -0.286582429  0.387354535
[184,] -2.736953838  2.44698209  -3.34028763 -1.568613996  0.178060678
[185,]  0.338263316  0.27677515  -1.88495779  3.976503405  1.011432756
[186,]  0.352564211 -1.04960926   0.06498237 -4.971245841 -0.274408408
[187,]  0.353004909  1.30518613   3.36813818  2.097370887  1.687837945
[188,] -1.244044162  1.39750194   2.70351108 -4.474848748  0.081778980
[189,] -1.700211152  0.44036192   2.15526818 -1.917432605 -4.003434265
[190,] -2.035177750  2.80155670   1.73064581 -1.214388917 -0.342158093
[191,] -5.547532550  3.93837340  -0.89920337 -7.335106519 -2.717992031
[192,] -1.460221692 -0.88501486  -5.87681497 -2.942708166  0.161341351
[193,]  0.961722529  1.38991107   8.80286438  2.616002913 -2.557853005
[194,]  2.214347745 -3.30445066  -1.22762469  1.522488858 -1.473894369
[195,]  1.620065716 -0.24054816   2.64295364  4.218583696  0.185096720
[196,] -1.227408842  0.88577361   0.26581336 -2.368281471 -0.881839493
[197,]  0.155338667 -0.14250861  -0.04645994  1.851947399 -0.037051356
[198,] -0.459061858  0.10462473  -1.87390984 -6.411140822  1.400857892
[199,] -0.694520689 -0.75712013  -7.36435276 -2.112017736  2.686898461
[200,] -3.344892271  3.47399781   2.79600716 -7.833607529 -0.478585229
[201,] -0.624971473 -0.32298286  -0.85635111 -3.924491443 -1.753538802
[202,] -0.884240085  4.08358830   8.52976885  2.728276595 -0.711998265
[203,]  0.700727008  0.15711997   0.61410121  2.551085992  0.827746481
[204,]  0.154278809 -1.79824642  -1.61683584 -5.037611373  0.374662705
[205,] -1.553180327  0.06141386   0.02890125 -4.136029963 -2.306099645
[206,] -3.181219859 -0.10972461  -8.20897666 -7.181658517 -0.258450399
[207,] -0.935133061  1.64282905   1.38589823 -0.729840238 -0.939128474
[208,]  1.831114475 -3.18325489  -4.13411906 -3.052935685  1.611732248
[209,]  1.144287310 -1.18933032   0.25302956 -2.663265708  0.772068708
[210,] -1.192231085  2.89415632   2.80890942 -1.706775755  0.919187866
[211,] -2.127404572  2.92488162   3.99033240 -6.393980384 -0.489896206
[212,]  1.145706132 -3.05915541  -4.10092831  0.698566290 -1.634464284
[213,]  0.750871800 -0.19491301   2.21466881  1.186712509  0.056557752
[214,] -3.629641646  0.34568417  -7.70226572 -4.141671149 -1.450793015
[215,]  1.277073903 -0.47998170   2.44622180 -2.375967798  1.399450367
[216,] -1.323042253  2.70358618   6.13059942 -4.514457224  0.078674836
[217,] -1.135750942  2.59192130  -2.10956157  0.602795480  4.253393005
[218,] -3.368382702  2.75990269  -2.67294345  1.068057641 -1.064839655
[219,] -1.988724767  0.75419550  -6.53140106  0.568113529  1.095115740
[220,]  1.042723619  0.19362201   2.42654374 -0.379203170  1.471769839
[221,]  0.820984114  0.83863596   6.04199664 -2.950220617  0.560736353
[222,] -0.210110178  1.71668452   0.44095658  1.227485429  1.827329030
[223,] -1.898618493  0.58316976   0.89084458 -5.425560570 -2.469987647
[224,] -1.620885336  2.92562800  -0.17190273 -2.266670565  3.181762083
[225,]  0.610960272  0.24091255   4.74621105 -1.569507610 -1.107369937
[226,] -3.324280691  3.95035733  -2.84756708 -2.478605530  1.822246870
[227,]  1.756425567 -2.31421852   1.97047380  1.608573183 -1.981915765
[228,]  1.085629522 -0.86739460   1.14967813 -0.851353677  1.213899980
[229,]  0.323756244  0.32831042   0.51457663  0.242940047  0.756571048
[230,] -0.311368031 -2.01731011  -3.51802420  0.200516501 -3.112481140
[231,] -1.518342585  1.44484379   2.13177564 -2.986985935 -0.783918987
[232,] -0.391705239  0.95189901  -3.92140869  0.431564688  3.921098644
[233,]  0.721321652 -0.16985848   0.26390365  4.392385471  1.083101715
[234,] -1.402181909  2.16169081   0.34730757  0.490046231  1.681774948
[235,]  2.797610502 -3.42982768  -0.99541882  2.240535190  1.952869078
[236,]  1.318374734 -0.94690250  -1.46075166 -1.716701396  2.807917815
[237,] -0.061978116 -0.48865435  -2.71922956  1.658286057 -0.208860255
[238,] -0.890303261 -1.70362355  -4.48914722 -0.030959319 -1.903352757
[239,]  1.386722376 -1.58175522   0.25129130  3.572145386 -1.291378785
[240,] -1.248373856  0.12209949  -7.40771738  3.517288952  1.783639576
[241,]  1.781306979 -1.89436762   1.42316317  2.100995962 -1.019021475
[242,]  0.867286598 -0.62990603   3.39713229 -3.880467978 -0.960963954
[243,]  0.701634129 -1.95931361  -5.74527247  4.010412586  0.733770740
[244,] -0.134160594  1.93498642   2.25789212 -0.663664345  3.136978786
[245,]  1.117007661  0.29218005   2.81089264  1.214930094  0.237181962
[246,] -0.163571270  1.93487838   1.48168037 -0.729958627  2.467306500
[247,]  1.603204798 -0.73166534   0.45776690  1.914726717  1.787368387
[248,]  0.627389803  0.08311801   1.92324503  2.658739210 -0.128082588
[249,] -0.292158526  0.14471305  -3.09897158 -1.574686221  2.229693399
[250,]  0.328585858 -1.10220297  -3.24398541 -0.829486390  1.412927069
[251,]  1.926458310 -0.84355496   3.46890080 -0.439768296  0.549839717
[252,] -0.725891219  0.40191246   2.98260968 -2.044653413 -3.133427408
[253,] -0.607710450 -2.53413954  -5.23505164 -4.637132882 -2.274741628
[254,] -2.103698405  2.45940378   3.02075799  1.675357997 -3.414551325
[255,] -0.225586982 -0.74327401  -1.10489255 -2.423700186 -1.133981234
[256,] -1.604960926  1.36813407  -1.60813907 -0.949498150  0.709507901
[257,]  1.436250841 -2.83379540  -2.91200010  0.867792631 -0.735860864
[258,]  1.164890750 -2.08355821  -3.23762097  5.661142963 -0.661255696
[259,]  1.988364179 -0.77004509   5.15399151 -0.903078216  0.167115604
[260,] -1.225337766  1.88887393   1.74264731  0.200001622  1.176704936
[261,]  0.446192467 -2.94805584  -1.43018066  2.470910491 -3.684637601
[262,]  2.689912114 -3.55287803  -0.85687441  1.478310197 -1.916540828
[263,]  0.497087041  0.54601295   5.93671045  4.739477965 -3.166752723
[264,] -1.922235722  3.61427560   4.12208927 -3.997272387  0.870646609
[265,]  0.479593685  1.58293686  -1.62042309  1.642360320  3.652031839
[266,] -1.556813174 -0.12538057  -8.30717223 -0.696003898  2.028122085
[267,] -0.108715879 -0.52979050  -2.50209847 -0.267681968 -0.572181294
[268,] -0.637718044 -0.54873407  -2.18622631  0.853850723 -1.294362187
[269,]  0.841023435  1.02531000   8.07091569 -1.049451519 -1.317437270
[270,] -1.285653245  0.45332510  -2.70923676 -0.217927181  0.135499640
[271,]  1.873753743 -1.04406539   1.96500080  7.457488606 -1.005369102
[272,]  0.609731558  0.28139867  -0.61275710 -0.166049408  2.733017433
[273,] -0.085750174  1.57430475   0.42435688 -1.776067221  4.190441367
[274,]  1.995328491 -2.74723881  -0.90621214  2.655604905 -0.636786698
[275,] -1.311980016 -0.11403795  -2.79788257  0.288676841 -2.523783377
[276,]  1.157722347 -1.58713396   1.32378126 -2.308566633 -0.488762765
[277,] -0.680674484 -0.53362250   0.61458752 -1.950861031 -3.011816539
[278,] -1.508111302  1.84634875  -0.24850290 -1.388872437 -0.708279462
[279,]  0.917503289  0.96166425   6.45625313 -0.182669522 -0.220772529
[280,] -0.131286923 -0.65849683  -1.66142904  1.015546067 -1.419135638
[281,] -2.443279508  2.66289820   1.28109032  0.364107106 -1.592415823
[282,]  2.077899985 -0.76467584   1.47480264  0.643675012  1.893380333
[283,] -0.677435882 -0.29153531  -3.72835302 -6.037722023  2.001394917
[284,] -1.279819023  1.88379607   3.05194091 -0.490761649 -0.965105154
[285,]  5.098147669 -3.04861785   2.83630296  4.004104008  1.592949264
[286,] -1.903442479  2.67486798   0.53053863 -1.888524436  0.533943337
[287,] -0.633666489  1.45793600   0.45360214 -0.400256970  2.784971434
[288,]  1.617040194 -0.61473643   0.72557184  2.590827271  0.195715452
[289,] -0.010099896  0.70009949  -2.04898449 -2.372586626  2.940754443
[290,] -0.758466253 -0.48994911  -6.68373824 -0.728579617  1.152924537
[291,]  1.235268927 -2.30048119  -2.03212533  2.378610230 -0.342972200
[292,] -1.234680132  1.87490211   7.83511682 -3.137706948 -3.077840850
[293,] -1.065001551 -0.54736858  -4.00568004 -3.120151894 -0.375687936
[294,] -0.088914755 -1.49473951  -0.62687055 -0.829804933 -2.209448012
[295,] -2.452828385  2.87682621   2.79190494 -4.384403451  0.447660634
[296,] -1.815351822  1.10199581   3.05352825 -4.649083891 -2.726953307
[297,]  1.807925994 -0.82115147   2.74251540 -0.188802930  1.020885383
[298,] -0.304032193 -0.16499358  -2.83130623  1.465999505  0.225503227
[299,]  0.595156556  0.27531047   1.74968004  5.562974983 -1.218959658
[300,] -0.763819064 -0.61622846 -11.71098900  1.589104858  3.932059637
[301,]  0.639933505 -1.81503418  -2.43195844 -0.217595465 -1.126743674
[302,]  3.155992792 -4.57037450  -3.64536760  4.577310069  0.175771257
[303,] -1.008786047  1.85252057   2.41666541  2.494205033 -1.644322533
[304,]  3.518387518 -3.99715383  -2.17088824  4.668459463  0.864170284
[305,]  1.717473997 -0.83654832   2.82867799 -2.553167247  2.378354565
[306,] -0.619120218  0.33686585   1.87208476  1.272511313 -1.727945095
[307,] -1.327377273  0.86678874  -0.47689028 -3.685773952 -0.528349320
[308,] -0.065950168  1.27758478   1.09711040 -0.961355781  1.694910998
[309,]  2.225957982 -2.16242976   0.81520148  1.791493691 -0.655169707
[310,]  2.372180630 -2.08577892   1.07245185  6.944274274 -0.690476707
[311,] -2.130678842  0.23813717  -7.51433918  2.291649218  0.400852748
[312,]  2.055920153 -1.92132979   4.74415103  1.870492910 -2.619274650
[313,] -1.943762615  3.77706214   9.63962012 -2.527886084 -2.678707584
[314,]  1.238412229 -0.68566472   2.55744451 -4.252015758  1.573804912
[315,]  2.646871703 -1.40081743   1.11923652  1.509227382  3.663564547
[316,] -3.486481935  0.88476557  -9.58250491 -3.637342623 -0.216403930
[317,] -1.602333107  2.78259843   2.78994948  0.990263816  0.313001420
[318,] -0.307297235 -0.76359911   3.91262819 -2.436084599 -3.771644110
[319,]  0.288983224  0.46790521   1.57282948 -1.754741300  1.440563446
[320,] -2.421072687 -0.55206787  -4.56546659 -2.178445148 -2.468911057
[321,] -1.034854908  1.33252193  -2.08407344 -3.341912773  2.444401857
[322,]  1.633508075 -2.35602794  -1.44990418  2.626797324  0.111481807
[323,]  0.751325653  0.31340754   2.75955008  3.420124310  0.278456093
[324,]  0.206698627  1.24130440   6.19058401 -0.936110397 -0.398753881
[325,] -3.999313478  4.23136601   4.23354049 -1.109452465 -3.766558150
[326,] -1.221379278  1.26541092   3.35560689  0.792217212 -1.856982287
[327,]  0.821707720 -1.85461571   0.04879872 -1.448323303 -0.693378083
[328,] -0.457003201 -0.09517645   0.05509922  1.716085229 -0.994111704
[329,]  1.518681942 -2.16611196   0.03196824  0.402142846 -2.149190593
[330,] -0.078764424  0.27447215   1.56343127 -3.223689459  0.228276401
[331,]  0.980184655  0.78893354   3.79765664  3.312416250  0.445732298
[332,]  2.175787136 -1.59825511   2.27583867  2.675454583 -0.941815234
[333,] -0.114824474 -1.46088070  -4.01256181 -4.266632911  0.596618046
[334,] -2.921921643  3.18307789  -0.10599158 -6.069361405  0.610689254
[335,] -1.712972554  3.80925754   4.86194068  3.218113121 -0.781627330
[336,]  0.619001087  0.41338944   4.95658515 -1.864141925  0.515941006
[337,] -0.606500971  2.25167419  11.16270929 -1.375347895 -2.775000999
[338,]  2.217958423 -0.41597303   3.85393115  9.174158953 -1.050081050
[339,]  2.085405179 -2.77380391  -2.10605650  0.517492202  0.652220023
[340,]  4.726543562 -2.23943913   5.78876139  3.278022792  2.123194126
[341,] -2.193801620  2.65936214   0.22722917  2.884829879 -0.057843555
[342,]  1.255253026 -2.62868053  -1.92930055 -5.822755593  0.629723738
[343,]  0.592399016 -1.40224298  -3.47257230  3.027561170 -0.716042834
[344,]  0.119966099 -1.96907792  -5.54392885  0.759160936 -1.262958907
[345,] -0.807795329  1.82916917   1.37180367 -0.700631781  0.612713211
[346,]  0.341164978 -0.22263240   1.28607058  3.159470306 -1.695398467
[347,] -0.467760081  1.37992781  -0.94029216 -1.174450403  2.502551226
[348,] -3.621849483  3.93868838   2.67555854  2.924786944 -3.018185572
[349,] -0.913433944  0.62292682  -5.38362835 -0.850070228  2.900332115
[350,]  2.369180780 -0.85399411  -1.56552181  1.395627624  4.604564518
[351,] -1.141166310  0.95922259  -1.20961032  2.810004557 -1.511761224
[352,] -1.957407586  2.16001640   0.98674034  0.057197504 -1.254935414
[353,] -0.059016915 -0.16286272   0.90454263  6.127529440 -3.587639624
[354,]  0.808748491 -4.21643192  -8.72608317 -1.279902761 -0.427942708
[355,] -0.309413314  1.58825665  -0.49292967 -1.210027705  2.167031756
[356,] -0.768731216 -0.36588248  -4.31732159 -0.609272753  1.072888543
[357,]  0.349265394 -0.63556690   3.03726666 -3.261147769 -1.440456428
[358,]  2.806627413 -2.84371830   1.26787649  2.411274187 -0.074980811
[359,]  0.205507415 -1.52870237  -2.46121902  0.215091466 -0.797082997
[360,]  1.149334108 -0.35483037   3.55963929  1.520816293 -1.541252796
[361,] -0.155231927  1.60073649   5.24594146 -5.101487678  0.566700250
[362,]  0.231085596  0.58465030  -0.31466202  0.133151575  2.456692048
[363,]  3.642889860 -1.24894240   6.77847665  2.612198913  0.060178000
[364,]  0.774405977  2.68802212   8.15859749  4.250771280  0.261355750
[365,] -1.842675006  0.96892505  -7.17223221  0.229415332  2.089404412
[366,]  1.798663460 -2.20271984  -0.40030975 -3.556785981  0.346707086
[367,] -3.198417808  2.28570286  -0.38243932 -4.793175982 -1.654825708
[368,]  3.356017638 -2.35434104   2.22901307 -0.224784910  3.168800151
[369,] -3.553630052  4.31139430  -1.78588198 -1.002743066  1.227130059
[370,] -0.053067836 -0.19655825   0.36167390 -0.389082192  0.186775247
[371,] -2.368857135  1.78347438  -2.77362299 -3.875264111  1.109180794
[372,]  1.904740742 -0.85812158   5.94417551  5.310516904 -2.106634047
[373,] -3.040036703  3.99715013   3.41579992 -5.886436531  0.258341193
[374,] -0.106209707  0.44482084  -2.47023872 -2.133275165  4.163105678
[375,] -1.033846218  0.61666593  -2.99052465  5.670501210  0.650497560
[376,]  0.622611719 -0.32746819  -1.11740659 -3.374015641  2.497470897
[377,]  1.317136510  1.04065530   5.07221361  1.845771103  3.061267426
[378,] -1.304341112  0.34187656   0.91246585 -3.000075358 -3.844585775
[379,] -0.291800494  1.67054416   5.56087556 -0.243651127 -0.751876164
[380,] -0.007734198  0.70998738  -1.00341784 -0.767996485  1.194602777
[381,] -0.456955104 -0.02919313   3.31263972 -0.290289319 -3.290770359
[382,]  0.881149914 -0.08303452   6.65092611  0.307703104 -1.776583847
[383,]  0.557504127 -0.42963564  -2.87069730  5.210829241  1.158852013
[384,]  2.264710077 -2.51315957   1.13151312 -0.421917337 -0.481671154
[385,] -1.653540298  3.16512800   4.51837513 -6.796621649  0.567640704
[386,] -1.506421789  2.40131247   0.50289595 -5.357770945  3.222042011
[387,]  0.375983176  1.09372189   3.08551851 -2.546430501  1.389371108
[388,] -1.841582929  1.30092823   5.11999852 -5.904575086 -3.193939475
[389,] -0.648636743  1.22265140  -1.34515387  1.266765545  2.162066617
[390,] -2.551332163  2.51976416  -0.49036230 -4.012865729 -0.635121968
[391,] -1.000019881  1.42347471   3.56327428  2.905751797 -1.897738426
[392,] -0.655573167 -0.57785985  -1.78788252 -2.973397562 -0.203445633
[393,] -1.984303232  0.44053236  -6.05663671  0.851814514 -1.305492880
[394,] -1.914193717  3.22628735   1.51544383 -2.402711887  1.437765858
[395,] -1.075368511 -0.57129014  -0.51760917 -4.396319609 -1.566183356
[396,] -3.259144372  2.89437278  -5.10500864  0.773307666  0.813983812
[397,]  3.251940596 -3.04554432  -2.40066236  5.767109852  2.916763413
[398,]  0.307578711 -0.64129445  -5.87574163 -0.753198949  2.493047187
[399,]  1.033377609  0.63334602   1.72819503  0.960919556  3.077539059
[400,]  1.270341353 -2.39520487  -0.27432219 -2.138915024 -2.486122999

# Correlation of decisionmaker random slopes
Omega <- cor(matrix(rnorm(P*(P+2)), P+2, P))

# Scale of decisionmaker random slopes
tau <- abs(rnorm(P, 0, .5))

# Covariance matrix of decisionmaker random slopes
Sigma <- diag(tau) %*% Omega %*% diag(tau)

# Centers of random slopes
beta <- rnorm(P)

# Individual slopes
beta_i <- MASS::mvrnorm(I, beta, Sigma) + W %*% t(Gamma)

# Again, quick plot to sanity check
plot(as.data.frame(beta_i))

# Create X -- let's make this a dummy matrix
X <- matrix(sample(0:1, Tasks*I*J*P, replace = T), Tasks*I*J, P)
# Each of the rows in this matrix correspond to a choice presented to a given individual
# in a given task

indexes <- crossing(individual = 1:I, task = 1:Tasks, option = 1:J) %>% 
  mutate(row = 1:n())

# Write a Gumbel random number generator using inverse CDF trick
rgumbel <- function(n, mu = 0, beta = 1) mu - beta * log(-log(runif(n)))
mean(rgumbel(1e6))

[1] 0.5779057

# Ok, now we need to simulate choices. Each person in each task compares each 
# choice according to X*beta_i + epsilon, where epsilon is gumbel distributed. 
# They return their rankings. 

# Modify the ranked_options data frame to include all necessary information
ranked_options <- indexes %>% 
  group_by(individual, task) %>% 
  mutate(
    fixed_utility = as.numeric(X[row,] %*% as.numeric(beta_i[first(individual),])),
    plus_gumbel_error = fixed_utility + rgumbel(n()),
    true_rank = rank(-plus_gumbel_error),
    true_order = order(true_rank),
    observed_order = case_when(
      true_order == J ~ J,  # Worst choice
      true_order == 1 ~ 1,  # Best choice
      TRUE ~ 3  # Tie all middle orders
    ),
    best_choice = as.numeric(true_order == 1),
    worst_choice = as.numeric(true_order == J)
  )

tt <- ranked_options %>% 
  group_by(individual, task) %>%
  summarise(start = min(row), 
            end = max(row)) %>% 
  ungroup %>%
  mutate(task_number = 1:n())

`summarise()` has grouped output by 'individual'. You can override using the
`.groups` argument.

Ok, onto new material; we’ll generate a “permutation matrix” of all the ways the J - 2 middle options can occur, which we’ll iterate over to calculate all the possible likelihoods.

For example, here with J = 5, the permutation matrix will have the following 6 permutations:

n_tied <- J - 2  
permutations <- permn(n_tied)
permutation_matrix <- do.call(rbind, permutations)
n_permutations <- nrow(permutation_matrix)

print(permutation_matrix)

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    3    2
[3,]    3    1    2
[4,]    3    2    1
[5,]    2    3    1
[6,]    2    1    3

Calculating these outside R is convenient in two senses. First, it’s more efficient- Stan isn’t really designed or equipped to work with permutations like this as they get more involved. Second, it can be helpful for intuition to mentally step through this matrix of permutations and consider the calculation needed for each entry, imagining it growing rapidly with the number of unknown middle options.

Here’s the Stan input data and model:

ROL with unknown middle preferences

stan_data_rol_ties <- list(
  N = nrow(X),
  T = nrow(tt),
  I = I, 
  P = P, 
  P2 = P2, 
  K = J, 
  rank_order = ranked_options$observed_order,
  X = X, 
  X2 = W, 
  task = tt$task_number, 
  task_individual = tt$individual,
  start = tt$start, 
  end = tt$end,
  n_tied = n_tied,
  permutations = permutation_matrix,
  n_permutations = n_permutations
)

rol_with_unknown_middle_options <- "
functions {
  real rank_logit_ties_lpmf(int[] y, vector delta, int[,] permutations, int n_permutations, int n_tied) {
    int K = rows(delta);
    vector[K] sorted_delta = delta[y];
    real out = 0;
    
    // Handle known best
    out += sorted_delta[1] - log_sum_exp(sorted_delta);
    
    // Compute normalizing factor once
    real normalizing_factor = log_sum_exp(sorted_delta[2:K]);
    
    // Handle tied middle ranks
    real perm_sum = 0;
    for (p in 1:n_permutations) {
      real perm_ll = 0;
      for (i in 1:n_tied) {
        int idx = permutations[p, i];
        perm_ll += sorted_delta[1+idx];
      }
      perm_sum += exp(perm_ll - n_tied * normalizing_factor);
    }
    out += log(perm_sum / n_permutations);
    
    return out;
  }
}

data {
  int N; // number of rows
  int T; // number of individual-choice sets/task combinations
  int I; // number of Individuals
  int P; // number of covariates that vary by choice
  int P2; // number of covariates that vary by individual
  int K; // number of choices
  
  int rank_order[N]; // The vector describing the index (within each task) of the first, second, third, ... choices. 
  matrix[N, P] X; // choice attributes
  matrix[I, P2] X2; // individual attributes
  
  int task[T]; // index for tasks
  int task_individual[T]; // index for individual
  int start[T]; // the starting observation for each task
  int end[T]; // the ending observation for each task
  
  int<lower=2> n_tied;
  int<lower=1> n_permutations;
  int permutations[n_permutations, n_tied];
}

parameters {
  vector[P] beta; // hypermeans of the part-worths
  matrix[P, P2] Gamma; // coefficient matrix on individual attributes
  vector<lower = 0>[P] tau; // diagonal of the part-worth covariance matrix
  matrix[I, P] z; // individual random effects (unscaled)
  cholesky_factor_corr[P] L_Omega; // the cholesky factor of the correlation matrix of tastes/part-worths
}

transformed parameters {
  matrix[I, P] beta_individual = rep_matrix(beta', I) + X2 * Gamma' + z * diag_pre_multiply(tau, L_Omega);
}

model {
  tau ~ normal(0, .5);
  beta ~ normal(0, 1);
  to_vector(z) ~ normal(0, 1);
  L_Omega ~ lkj_corr_cholesky(4);
  to_vector(Gamma) ~ normal(0, 1);

  for(t in 1:T) {
    vector[K] utilities;
    utilities = X[start[t]:end[t]] * beta_individual[task_individual[t]]';
    rank_order[start[t]:end[t]] ~ rank_logit_ties(utilities, permutations, n_permutations, n_tied);
  }
}
"

Summing up all the possible middle options should be fairly transparent, but if you’re interested in the log_sum_exp tricks I’m using to speed this up, see this footnote¹².

To reinforce the intuition here, this is choosing to handles unknown middle options with maximum precision given we know very little about them, at (usually) great computational cost. Before we compare this to MaxDiff, and see if it is meaningfully more accurate with the same data, I want to also introduce an approximation to this full marginal likelihood that won’t choke even on big choice sets.

J! is pretty bad; is there an approximation?

As we’ve seen, calculating all possible orderings for unknown middle preferences can quickly become computationally intensive as the number of choices increases. Fortunately, there’s a clever approximation that can help us out. This approximation comes to us via a deep connection between discrete choice models and survival analysis. To gesture at the analogy here- tied times in survival models function like our unknown middle options, and the survival analysis literature has decades of experience handling this type of tie given how common they are in survival data¹³. We’ll look at my favorite of these, the Efron approximation, originally developed for tied times in Cox proportional-hazards models.

Mathematically, the Efron approximation for our rank-ordered logit with unknown middle preferences looks like this:

L(\mu) = \prod_{i=1}^N \frac{\exp(\mu_{ij})}{\sum_{j \in C_i} \exp(\mu_{ij})} \cdot \prod_{k=2}^{J-2} \frac{\exp(\mu_{ik})}{\sum_{j \in C_i} \exp(\mu_{ij}) - \frac{k-1}{J-2} \sum_{l \in M_i} \exp(\mu_{il})}

Where C_i is the full choice set for individual i, M_i is the subset of C_i containing the tied (unranked) middle alternatives, and J is the total number of alternatives in C_i. I recognize this is super dense, and that grokking it will require a bit of reading that’s outside the scope of this post; here’s a footnote with some helpful references¹⁴.

Intuitively, this approximation attempts to more continuously approximate the full marginal likelihood we just discussed without calculating every possible ranking. It does this by adjusting the denominator of the likelihood to account for the tied middle alternatives in a stepwise fashion.

Here’s the code implementation, which to me gives a lot more intuition than the equation above:

ROL with Efron approximation of unknown middle options

efron_model <- "
functions {
  real rank_logit_efron_lpmf(int[] y, vector delta, int n_tied) {
    int K = rows(delta);
    vector[K] sorted_delta = delta[y];
    real out = 0;
    
    // Handle known best
    out += sorted_delta[1] - log_sum_exp(sorted_delta);
    
    // Handle tied middle ranks using Efron approximation
    real tied_sum = log_sum_exp(sorted_delta[2:(K-1)]);  // Changed to log_sum_exp
    real normalizing_factor = log_sum_exp(sorted_delta[2:K]);
    for (i in 1:n_tied) {
      real d = (i - 1.0) / n_tied;
      real adjusted_factor = log_diff_exp(normalizing_factor, log(d) + tied_sum);
      out += sorted_delta[1+i] - adjusted_factor;
    }
    
    // Handle known worst
    // out += sorted_delta[K] - log_sum_exp(sorted_delta[2:K]);
    
    return out;
  }
}

data {
  int N; // number of rows
  int T; // number of individual-choice sets/task combinations
  int I; // number of Individuals
  int P; // number of covariates that vary by choice
  int P2; // number of covariates that vary by individual
  int K; // number of choices
  
  int rank_order[N]; // The vector describing the index (within each task) of the first, second, third, ... choices. 
  matrix[N, P] X; // choice attributes
  matrix[I, P2] X2; // individual attributes
  
  int task[T]; // index for tasks
  int task_individual[T]; // index for individual
  int start[T]; // the starting observation for each task
  int end[T]; // the ending observation for each task
  
  int<lower=2> n_tied; // number of tied ranks
}

parameters {
  vector[P] beta; // hypermeans of the part-worths
  matrix[P, P2] Gamma; // coefficient matrix on individual attributes
  vector<lower = 0>[P] tau; // diagonal of the part-worth covariance matrix
  matrix[I, P] z; // individual random effects (unscaled)
  cholesky_factor_corr[P] L_Omega; // the cholesky factor of the correlation matrix of tastes/part-worths
}

transformed parameters {
  matrix[I, P] beta_individual = rep_matrix(beta', I) + X2 * Gamma' + z * diag_pre_multiply(tau, L_Omega);
}

model {
  tau ~ normal(0, .5);
  beta ~ normal(0, 1);
  to_vector(z) ~ normal(0, 1);
  L_Omega ~ lkj_corr_cholesky(4);
  to_vector(Gamma) ~ normal(0, 1);

  for(t in 1:T) {
    vector[K] utilities;
    utilities = X[start[t]:end[t]] * beta_individual[task_individual[t]]';
    rank_order[start[t]:end[t]] ~ rank_logit_efron(utilities, n_tied);
  }
}
"

The beauty of the Efron approximation is that it works remarkably well in practice, especially compared to alternatives. As discussed for example in Hertz-Picciotto and Rockhill (1997), the Efron method tends to outperform other approximations, particularly when dealing with moderate to heavy ties (read: unknown middle rankings).

Evaluating ROL with Unknown Middle Preferences against MaxDiff and ROL

So now we’ve got 3 candidate models for best/worst choice data- MaxDiff, ROL calculating the full marginal likelihood for the middle options, and ROL using the Efron approximation. I’ll also add in the ROL without ties to give a sense of the information lost when just collecting best and worst choices, though of course my main goal is compare the quality of the 3 models on best/worst choice data.

The next chunk just sets up the comparison between the four models on data with J = 5, so feel free to jump to the summary plot after it.

4 way Comparison Setup

compiled_maxdiff_model <- stan_model(model_code = best_worst)
compiled_rol_ties_model <- stan_model(model_code = rol_with_unknown_middle_options)

Warning in readLines(file.path(makevar_files)): incomplete final line found on
'C:/Users/timma/Documents/.R/Makevars.win'

compiled_efron_model <- stan_model(model_code = efron_model)

Warning in readLines(file.path(makevar_files)): incomplete final line found on
'C:/Users/timma/Documents/.R/Makevars.win'

compiled_rol_model <- stan_model(model_code = ranked)

# Prepare data for Stan models
stan_data_maxdiff <- list(
  N = nrow(X),
  T = nrow(tt),
  I = I,
  P = P,
  P2 = P2,
  K = J,
  choice = ranked_options$best_choice,
  worst_choice = ranked_options$worst_choice,
  X = X, X2 = W, 
  task = tt$task_number, 
  task_individual = tt$individual,
  start = tt$start, end = tt$end
)

stan_data_efron <- list(
  N = nrow(X),
  T = nrow(tt),
  I = I,
  P = P,
  P2 = P2,
  K = J,
  rank_order = ranked_options$observed_order,
  X = X, X2 = W, 
  task = tt$task_number, 
  task_individual = tt$individual,
  start = tt$start, end = tt$end,
  n_tied = J - 2  # number of tied ranks
)

stan_data_rol <- list(
  N = nrow(X),
  T = nrow(tt),
  I = I,
  P = P,
  P2 = P2,
  K = J,
  rank_order = ranked_options$true_order,  # Use true_order for full ranking
  X = X, X2 = W, 
  task = tt$task_number, 
  task_individual = tt$individual,
  start = tt$start, end = tt$end
)

fit_maxdiff <- sampling(compiled_maxdiff_model, data = stan_data_maxdiff, iter = 800)

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

fit_rol_ties <- sampling(compiled_rol_ties_model, data = stan_data_rol_ties, iter = 800)
fit_efron <- sampling(compiled_efron_model, data = stan_data_efron, iter = 800)
fit_rol <- sampling(compiled_rol_model, data = stan_data_rol, iter = 800)

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

normalize_utilities <- function(utilities) {
  (utilities - mean(utilities)) / sd(utilities)
}

process_results <- function(fit, dataset_name) {
  as.data.frame(fit, pars = "beta_individual") %>%
    gather(Parameter, Value) %>%
    mutate(
      individual = str_extract(Parameter, "[0-9]+(?=,)") %>% parse_number(),
      column = str_extract(Parameter, ",[0-9]{1,2}") %>% parse_number()
    ) %>%
    group_by(individual) %>%
    mutate(Value_normalized = normalize_utilities(Value)) %>%
    group_by(individual, column) %>%
    summarise(
      median = median(Value_normalized),
      lower = quantile(Value_normalized, 0.05),
      upper = quantile(Value_normalized, 0.95),
      .groups = 'drop'
    ) %>%
    ungroup() %>%
    mutate(
      True_value = as.numeric(t(beta_i)),
      Dataset = dataset_name
    ) %>%
    group_by(individual) %>%
    mutate(True_value_normalized = normalize_utilities(True_value)) %>%
    ungroup()
}

results_maxdiff <- process_results(fit_maxdiff, "MaxDiff")
results_rol_ties <- process_results(fit_rol_ties, "ROL with ties")
results_efron <- process_results(fit_efron, "Efron Approximation")
results_rol <- process_results(fit_rol, "ROL (no ties)")

all_results <- bind_rows(results_maxdiff, results_rol_ties, results_efron, results_rol)

# Calculate and print RMSE; couldn't be bothered to find one in a package
calculate_rmse <- function(data) {
  sqrt(mean((data$True_value_normalized - data$median)^2))
}

rmse_comparison <- data.frame(
  Model = c("MaxDiff", "ROL with ties", "Efron Approximation", "ROL (no ties)"),
  RMSE = c(
    calculate_rmse(results_maxdiff),
    calculate_rmse(results_rol_ties),
    calculate_rmse(results_efron),
    calculate_rmse(results_rol)
  )
) %>%
  arrange(desc(RMSE))


fourway_comparison_plot <- ggplot(all_results, aes(x = True_value_normalized, y = median, color = Dataset)) +
  geom_linerange(aes(ymin = lower, ymax = upper), alpha = 0.3) +
  geom_point(alpha = 0.5) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  labs(x = "Normalized True Utility",
       y = "Normalized Estimated Utility",
       title = "Comparison of Model Performances",
       subtitle = "With 90% credibility intervals") +
  scale_color_manual(values = c("MaxDiff" = "blue", 
                                "ROL with ties" = "green",
                                "Efron Approximation" = "red",
                                "ROL (no ties)" = "purple")) +
  theme_minimal() +
  facet_wrap(~Dataset, ncol = 2) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))

Before we show the comparison plot, here’s a minor technical caveat. The natural scale of the utilities the two unknown handling models produce is different than that of the MaxDiff and ROL without ties models. To plot them usefully against one another, I thus have to normalize the utilities. This makes it a bit harder to compare these plots to the last few, but unfortunately its necessary for a valid comparison.

print(fourway_comparison_plot)

So roughly from the plot: ROL without unknowns > ROL with full marginal likelihood for unknowns ~ ROL with Efron Approximation for Unknowns > MaxDiff. Let’s sanity check that by looking at RMSEs.

print(rmse_comparison)

                Model      RMSE
1             MaxDiff 0.4370070
2       ROL with ties 0.2507489
3 Efron Approximation 0.2499740
4       ROL (no ties) 0.1625324

I flagged this briefly above, but to say it again: J = 5 has few unknown middle options, meaning these are fairly favorable conditions for the Efron approximation. We’d expect higher J to widen the gap between the full marginal likelihood and the Efron approximation.

So practically, the takeaway here is that you can do better than modeling with MaxDiff even if you designed your discrete choice experiment to just collect best/worst choice data. If you have the time and/or your number of alternatives isn’t super large, you can do better with the full marginal likelihood model. If you want a workable and fast approximation, the Efron model still looks meaningfully better than MaxDiff.

Beyond giving a better option for existing MaxDiff data, I’ll admit I had a secondary motive in wanting to explore these “unknown middle choices” style models. I’ve sometimes seen MaxDiff discussed as an imperfect modeling strategy that becomes desirable because it gives access to the best/worst choice data collection strategy.

To flesh out that claim, the idea is that MaxDiff has known flaws that aren’t ideal, but a study only needing to collect best and worst options hits a sweet spot of respondent experience, simplicity to setup, and statistical efficiency (via being able to repeat many choice sets). Better a slightly worse model on significantly more informative-dense data, since design trumps analysis.

This is probably the type of argument for MaxDiff I’m most sympathetic to, in the sense that I do genuinely do believe choice experiments place a lot of burden on respondents in a way that may bias results. Perhaps selecting only the best and worst option meaningfully helps with that. Further, I might even find it plausible that if you need to compare a larger than ideal set of options, there might be value in asking a bunch of best/worst of 10 choice sets instead of trying to do tons of best of threes to enable the full ROL with ties option I mentioned above. Even if you think this way though, the ROL with unknown middle options models I’ve just shown are likely the more accurate model for the design you find significant value in. You can’t justify using MaxDiff via the designs it enables.

Conclusions

In summary:

1. There are some significant psychological and mathematical limitations entailed in using a MaxDiff model of discrete choice.
2. These issues with the MaxDiff model are more than theoretical, and can cause significant harm to the quality of model predictions.
3. If you have control over the design of the choice experiment, you can do much better than MaxDiff by designing the choice sets in a connected graph, and then using a rank-ordered logit. This gets most of the benefit of ROLs more generally, while preserving the relatively low burden respondent experience of just choosing a best and worst choice.
4. Even if you don’t have control over the design of the experiment, you can still outperform MaxDiff using ROL models that correctly handle unknown middle choices. Crucially, the existence of these models negates the argument for using MaxDiff based solely on its compatibility with best/worst choice designs.

In writing this post, I’ve aimed to both give intuition and concrete examples of models that I would’ve benefited from when I was initially trying to find alternatives to the MaxDiff formulation of discrete choice. My hope is gathering both the objections to MaxDiff and alternatives in one place can make it a bit easier for others to pivot to the better models available.

There are plenty of ways you can extend what I’ve done in this post, and I don’t want to claim this is in any way some definitive model for discrete choice. Just to name a few- one could improve the priors in the models, design much more efficient graphs of options, compare the efficiency of connected graphs to ROL with unknowns for large choice sets via simulation, or even implement these models in something like Blackjax with a Gibbs sampler to make them blazing fast. My primary goal here has been to show that the MaxDiff likelihood isn’t a great foundation for discrete choice.

Thanks for reading!

Footnotes

1. The smiley face is threatening: it’s feels like it’s there to make sure MaxDiff stays down :-D

2. If you want more of an in-depth introduction, here are some resources I’d recommend. First, Jim Savage’s blog post series is great, and does a fantastic job of explaining how our assumptions about choicemaking map onto the math. I also benefited from reading Glasgow’s Interpreting Discrete Choice Models, which builds up the basics of choice models a bit more slowly. Finally, if you want something that goes much more in detail, Kenneth Train’s Discrete Choice Models with Simulation goes into great mathematical detail about the models.

3. I won’t rehash more subtle implications of this basic model here; see the footnote above for references. Instead, I’ll mostly pull up assumptions as they become relevant for discussing MaxDiff.

4. Other options are possible, but less common, and would take us too far afield, so I’ll skip explaining the choice of Gumbel versus other distributions here.

5. President Biden is great, it’s just hard to compete with the comically malignant and incompetent by being solid and stabilizing.

6. Unlike the blog post about left-handed kangaroos drinking fosters I’m working on, this post was not primarily motivated by the urge to ask DALL-E to create this.

7. I won’t re-explain the full details of the model here, since the original source already does a fantastic job. If you haven’t already read the original post series, it’s well worth taking the time to do so now, as there are a lot of moving pieces here- setting sensible priors for this model, incorporating covariates logically, and setting this up to run efficiently in Stan.

8. I will note that if you’re following along at home, sampling variability means that for some realizations of this DGP + these models at much lower N, MaxDiff performs just slightly better. But “barely better in a fraction of runs” isn’t really a good outcome from gathering twice as much data.

9. A great place to start reading about efficient design of discrete choice experiments is the literature review section of Traets et al, 2020. For more on the desirable properties of designs like the BIBD I actually chose, Raghavarao’s 1988 book Constructions and Combinatorial Problems in Design of Experiments is wonderful. Part of what I find appealing about this broader literature is that results not just from the ’80s, but stretching back a century or more remain relevant for building good choice experiments. While the last decade or so has seen numerous new techniques pop up that squeeze out additional bits of accuracy in the right conditions, the core challenge remains interwoven with the combinatorics of experimental design that were recognized by the earliest statisticians.

10. Just for fun, notice what happens if there’s only 1 middle term. It’s your friendly neighborhood complete rank-ordered logit, since there’s only 1 way 1 option can occur!

11. To be transparent, note that J = 5 isn’t a totally neutral comparison ground between the full marginal likelihood model and the Efron approximation. The Efron will get somewhat worse as the number of middle options grows.

12. Woo optimization gang! Ok, so the idea here is similar to this block in the ROL without unknowns model, where at each choice we normalize by subtracting out out += tmp[i] - log_sum_exp(tmp[i:]), but it has to function differently given the unknowns. More specifically, the right normalizing factor for all the middle J -2 options is log_sum_exp(sorted_delta[2:K]) for each possible permutation, not each possible middle option. So as the number of permutations grows, we can save a bunch of compute mostly through pre-computing the normalizing factor, but also by doing n_tied * normalizing_factor instead of subtracting n_tied times per permutation.

13. See the Allison/Christakis paper above for a more fleshed out version of this comparison. It’s very elegant, but I don’t think there’s a way to do that analogy justice without a massive detour in this post. Similarly, I’d love to walk through the reasoning that leads to the Efron approximation and how it compares to competitors, but this post is already pretty long.

14. I really appreciate this set of slides for pumping intuitions about different approaches to approximating the middle options. These course notes are also super helpful, especially if you want to see the Cox proportional hazards model discussed more in depth before you try to understand the tie handling in those terms. Finally of course, the original paper is great, but very very dense.

Reuse

Citation

BibTeX citation:

@online{timm2024,
  author = {Timm, Andy},
  title = {Better Discrete Choice Modeling Through the Rank Ordered
    Logit},
  date = {2024-08-06},
  url = {https://andytimm.github.io/posts/doing_maxdiff_better/better_maxdiff.html},
  langid = {en}
}

For attribution, please cite this work as:

Timm, Andy. 2024. “Better Discrete Choice Modeling Through the Rank Ordered Logit.” August 6, 2024. https://andytimm.github.io/posts/doing_maxdiff_better/better_maxdiff.html.