Rank-ordered data - dealing with increased randomness among lower ranks What is the best way to analyze rank-ordered data when there are signs that respondents were less diligent/able to assign lower ranks? Is it sufficient to introduce a dummy for lower/earlier ranks into the model or should one model this directly?
I have a set of survey data containing rank-orderings. More specifically, a large sample of respondents was asked to indicate which groups of people should first get access to a COVID-19 vaccine once it becomes available. They ranked 7 different groups from highest to lowest priority; ties were not permitted, rankings had to be complete and the initial ordering of the groups was randomized across respondents. 
In my analysis so far, I followed Allison & Christakis (1994, "Logit Models for Sets of Ranked Items", Sociological Methodology). Descriptive analyses show a quite consistent rank-ordering across a number of demographics and non-parametric tests (Friedman & Wilcoxon) indicate significant differences between pairs of choices.
The problem is: It turns out that respondents seem to have been less diligent in assigning the lower ranks, or at least less able to do so (not unrealistic, IMO). Following Allison/Christakis (pp. 216-218), I created a dummy for higher (1-4) and lower (5-7) ranks and included it as an individual-specific covariate into the model. The interactions between the dummy and the alternatives were significant and the coefficient estimates for later choices also show signs of greater randomness (the standard deviation between the coefficients for early choices is 3 times the SD for late choices).
My question is: What is the best way to deal with this? If I understand Allison/Christakis correctly, I should treat later choices as ties by, it seems to me, keeping the early/late dummy in the model - or, in Stata at least, code all later choices as 0 and use the incomplete() option for rologit. 
One alternative is to use models that allow for unobserved heterogeneity in ranking capabilities (Fok. et al., "A Rank-Ordered Logit Model with Unobserved Heterogeneity in Ranking Capabilities", Journal of Applied Econometrics), but I am not aware where this would be implemented (not in Stata and the gmnl package for R by Sarrias & Daziano cannot deal with ranked data, as far as I can tell). A final alternative, but hardly a great one, would be to use only the most-preferred choice as the outcome and use the regular conditional logit (or multinomial logit, as it is also known).
Any better suggestions or pointers to relevant packages would be greatly appreciated!
 A: I think you might be well served by applying the multistage models of Plackett and Luce, later extended by Benter. Briefly, these models are constructed as sequences of conditional multinomial distributions.
In more detail, I recently proposed a method that further extends these models and so I have been doing quite a bit of reading about them. Here are the standard Plackett-Luce modeling assumptions as I write out in my paper.

For $i=1,\ldots,n$, the $i$th ranker's ordered list of $\ell_i$ items is denoted by ${\bf x_i} = \{x_{i1},x_{i2},\ldots,x_{i\ell_i}\}$, with $x_{is}\in\{1,\ldots,v\}$ and $s=1,\ldots,\ell_i$ indexing each stage. If the lists are complete, then $\ell_i\equiv v$ for all lists; if they are partial, then $\ell_i \equiv \ell < v$ for all $i$, where $\ell$ is artificially chosen and external to the modeling process; if they are ragged, then $\ell_i \leq v$ for each $i$, with potentially different values of $\ell_i$ for each $i$.

and later

The $i$th ranker generates an ordered list of length $v$ from among a pre-specified, fixed-length set of items, starting with his/her/its most-preferred item. Define $\mathcal{O}_{is}$ to be the set of items yet to be ranked just before the $s$th stage:
\begin{align}
\mathcal{O}_{is} = \begin{cases} 
\{1, \ldots, v\}, & s=1\\
\{k: k \not\in \{x_{is'}\}_{s'<s}\}, & s>1
\end{cases}\Bigg\},\label{ois}
\end{align}
and let $1_{[X]}$ be 1 when the statement $X$ is true and 0 otherwise. The Plackett-Luce (PL) probability that item $k\in\{1,\ldots,v\}$, is ordered $s$th is $\Pr(x_{is} = k|\mathcal{O}_{is}) = 1_{[k\in\mathcal{O}_{is}]}\exp(\theta_k)/\sum_{j\in \mathcal{O}_{is}}\exp(\theta_j)$, i.e. proportional to $\exp(\theta_{k})$ until it gets ordered, and zero afterwards. There are $v$ parameters, $\Theta = \{\theta_1,\theta_2,\ldots,\theta_v\}$. Of these, $v-1$ are identified, and without loss of generality, we may assume that $\min_j\{\theta_j\}\equiv0$.

An important extension that I think is appropriate for your situation where rankers are more ambivalent at later stages comes from Benter, who proposed to dampen the weights towards zero so that at later stages (later ranks), differences in the log-likelihood are smaller. Let a dampening function $\delta(s)$ map the set of integers $s\in\{1,\ldots,v-1\}$ to the interval $(0,1]$, with $\delta(1)\equiv 1$ for identifiability. From my paper again:

...the  Benter-Plackett-Luce (BPL) model for the probability of selecting item $k$ at the $s$th stage conditional on the choices from the previous $s-1$ stages is $\Pr(x_{is} = k|\mathcal{O}_{is}) = 1_{[k\in\mathcal{O}_{is}]}\exp(\theta_k\delta(s))/\sum_{j\in \mathcal{O}_{is}}\exp(\theta_j\delta(s))$, for $k=1,\ldots,v$ and $s = 1,\ldots,\ell_i$. To be estimated are the $v-1$ identified parameters in $\Theta$ plus the number of parameters in the chosen functional form of $\delta(\cdot)$

I have not yet put my R code into a package, but the scripts for fitting BPL models are on my github repo, along with a few examples and vignettes, which should hopefully be useful for you.
If you're interested, what I did in my paper to extend the BPL models was equip this BPL log-likelihood with an $L_0$ variable selection penalty, so that when maximizing the penalized log-likelihood, some of the item weights ($\theta_i$) get forced to zero, and you can thus obtain a sparse consensus list that doesn't necessarily include all of the items that were ranked.
EDIT Since your data are comprised of complete rankings, you don't need the $\theta_0$ parameter that I introduce. All parameters are logged, so the way to drop $\theta_0$ from the likelihood would be to set it equal to $-\infty$, which you do by setting fixed = matrix(-Inf,dimnames = list(c("0"))) in the call to penalized_rank_path.
