Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
72,641 characters · 23 sections · 42 citation commands
csranks: An R Package for Estimation and Inference Involving Ranks
In economics and other social sciences, it is often desirable to rank populations according to some performance measure or to rank observations before running regressions. For instance, it may be desired to rank neighborhoods according to some measure of intergenerational mobility, countries according to some measure of academic achievement, or hospitals according to patients' average waiting times (Mogstad:2023aa). A prominent example for regressions involving ranked observations is the study of intergenerational mobility in which the slope coefficient of a rank-rank regression is a popular measure of the persistence in socioeconomic status across generations (Deutscher:2023oo and Mogstad:2023uu).
In this paper, we introduce the \proglang{R} package \pkg{csranks} and show how it can be used to perform inference on ranks as well as in regressions involving ranks.
First, we review the statistical methods proposed by Mogstad:2023aa and Mogstad:2021bb for the construction of confidence sets for ranks. In this context, there are several populations (e.g., neighborhoods, countries, hospitals) that we want to rank according to some estimated performance measure. Since the performance measure is estimated (e.g., because it is computed on a random sample of data from the population), the statistical uncertainty in the performance measure transfers into statistical uncertainty in the ranking of the populations according to these (estimated) performance measures. Mogstad:2023aa propose (i) marginal confidence intervals for the rank of a single population, (ii) simultaneous confidence intervals for the ranks of all populations, and (iii) confidence intervals for the $\tau$-best populations. We show how each of these can be computed using functions from the \pkg{csranks} package. We then provide an empirical illustration using data from the PISA study for ranking countries according to their students' scholastic performance.
Second, we review the statistical methods proposed by Chetverikov:2023aa for inference on regressions involving ranks. The main specification, which is popular in empirical work in economics, is a rank-rank regression in which both the independent and the dependent variable are first transformed into ranks and then one of the ranked variables is regressed on the the other, possibly including further (non-ranked) covariates. Inference in such regressions is nonstandard because both the independent and dependent variables have been estimated. The OLS estimator of the regression coefficients is asymptotically equivalent to a U-statistic. Chetverikov:2023aa show that the estimator is asymptotically normal and derive the asymptotic variance. In addition, this paper provides asymptotic normality results for several related regression specifications that involve ranked variables. We show how the \pkg{csranks} package can be used to compute all of these estimators and how to perform inference on the regression coefficients, e.g., by computing standard errors and confidence intervals.
Many statistical tests and software packages involve ranks because they consider rank-based statistics, e.g. Wilcoxon's statistic, which are not the subject of this paper. The only statistical software related to inference on ranks that we are aware of is the the \proglang{R} package \pkg{ICRanks} by Al Mohamad et al., based on mohamad2017simultaneous, mohamad2017improvement and al2022simultaneous. It implements a number of alternative methods for construction of confidence sets for ranks. However, the package's scope is restricted to simultaneous confidence sets in the case when the performance measures are independent and follow a Gaussian distribution. In contrast, the methods in the \pkg{csranks} package do not require these two assumptions and, in addition, are weakly more powerful as they implement stepwise improvements. For inference in rank-rank regressions, applied researchers typically employ standard covariance estimators for OLS regressions such as the homoskedastic or heteroskedasticity-robust estimators implemented in the \proglang{R} package \pkg{sandwich} (zeileis2020various,zeileis2004econometric). As shown by Chetverikov:2023aa, these estimators do not lead to valid inference in rank-rank regressions and the \pkg{csranks} package implements the valid inference procedures proposed by Chetverikov:2023aa. The \proglang{R} packages \pkg{copula} (copula1,copula2,copula3,copula4) and \pkg{rvinecopulib} (rvinecopulib) implement methods for estimation and inference on parameters of copulas. In a special case, the slope of the rank-rank regression is equal to Spearman's rank correlation, a feature of the copula of the independent and dependent variables in the regression, but in general the parameters considered in \pkg{csranks} cannot be written as a feature of the copula and thus \pkg{copula} and \pkg{rvinecopulib} cannot be used for inference in rank-rank regressions.
Suppose we want to create a ranking of $p$ populations, e.g. countries, political parties or hospitals, according to some performance measure $\theta_1,\ldots,\theta_p$. The rank of a population can be defined in different ways. First, ranks can be defined such that the population with the largest performance measure is assigned rank 1, the second largest is assigned rank 2 and so on. We will refer to this as a “decreasing” ranking as the rank of population $j$ is decreasing in its own performance measure $\theta_j$. Alternatively, the population with the largest performance measure could be assigned the rank $p$, the second largest is assigned the rank $p-1$ and so on (“increasing” ranking). Second, if there are ties among the performance measures $\theta_1,\ldots,\theta_p$, then one needs to decide which rank to assign to tied populations. For instance, consider three populations with $\theta_1=10$, $\theta_2=\theta_3=20$. In a decreasing ranking, population 1 is assigned rank $3$. The populations 2 and 3 are tied with the largest performance measure. So, both could be assigned rank 1, both could be assigned rank 2, or any value between 1 and 2.
Let $\theta := (\theta_1,\ldots,\theta_p)'$. A general definition of an increasing rank for population $j$ is
where $\omega\in[0,1]$ is a parameter that describes how ties are handled. If none of the other populations are tied with population $j$ ($\theta_j\neq\theta_k$ for all $k$), then the rank of $j$ is equal to $1+ \sum_{k=1}^p 1\{\theta_k< \theta_j\}$ and does not depend on $\omega$. A decreasing rank is obtained by multiplying all performance measures by $-1$:
Sometimes, it is useful to scale the integer ranks defined above back to the $[0,1]$ interval. A common way to do this in practice is to divide the integer ranks by $p$. For instance, the increasing (“fractional”) rank for population $j$ is then
which can be expressed as a weighted average of the empirical cdf $\hat{F}_{\theta}(\theta_j)$ and $\hat{F}_{\theta}^-(\theta_j)$. When $\omega=1$, then the rank corresponds to a popular definition of the rank in applied work, namely the rank of $j$ is equal to the empirical cdf evaluated at population $j$'s performance measure.
For concreteness, in this section we consider the ranking of $j=1,\ldots,p$ countries according to how well they educate their children (as in our empirical illustration in Section (ref)). The true performance measures for the countries are $\theta_1, \ldots, \theta_p$. These are not observed directly. Instead, for each country we observe data of sample size $n$ from which we compute estimators $\hat{\theta}_1,\ldots,\hat{\theta}_p$ of the performance measures $\theta_1, \ldots, \theta_p$. These estimators may be sample averages of children's test scores, for example. However, the methods below are theoretically justified for the general case in which $\hat{\theta}:=(\hat{\theta}_1,\ldots,\hat{\theta}_p)'$ is a consistent estimator of some parameter $\theta=(\theta_1, \ldots, \theta_p)'$ as long as we can construct a consistent estimator $\hat{\Sigma}$ of the $p\times p$ asymptotic covariance matrix $\Sigma$, with $(j,k)$-element denoted by $\hat{\sigma}_{jk}$, such that
as $n\to\infty$. In this section, we focus on decreasing ranks, i.e., (ref), with $\omega=0$. The goal is to use data from each country to form confidence sets $R_{n,j}$ that cover the rank of country $j$, i.e., $R_j^{\theta}$, with probability (approximately) no less than a prespecified level (e.g., $0.95$).
The goal in this subsection is to construct a two-sided confidence set $R_{n,j}$ for the rank of a particular country $j$ that satisfies\footnote{In fact, the construction described in this section satisfies a stronger property, $\liminf_{n \rightarrow \infty} \inf_{P\in\mathbf{P}} P\left\{R_j^{\theta}(P) \in R_{n,j}\right\} \geq 1 - \alpha$, i.e., asymptotic uniform coverage. For details, see Mogstad:2023aa.}
for some pre-specified confidence level $1-\alpha$. This requirement means that the set $R_{n,j}$ covers the true rank of country $j$ with asymptotic probability no less than $1-\alpha$. The construction is based on simultaneous confidence sets for the differences of performance measures as in Mogstad:2023aa and Mogstad:2021bb. For concreteness, we explain one particular approach based on the parametric bootstrap which exploits the asymptotic normality in (ref), but other constructions are possible; see Mogstad:2023aa. To this end consider the confidence set
where $\hat{se}_{jk}^2 := \hat{\sigma}_{jj} + \hat{\sigma}_{kk}-2\hat{\sigma}_{jk} $ is an estimate of the variance of $\hat{\theta}_j-\hat{\theta}_k$ and $c_{{\rm symm}, n,j}^{1-\alpha}$ is the $(1-\alpha)$-quantile of $$\max_{k\colon k\neq j} \frac{|\hat{\theta}_j - \hat{\theta}_k - (\theta_j-\theta_k)|}{\hat{se}_{jk}}.$$ This quantile can be simulated using a parametric bootstrap based on (ref) as follows. Generate $m$ draws of normal random vectors $Z:=(Z_1,\ldots,Z_p)' \sim N(0,\hat{\Sigma}))$. The desired quantile $c_{{\rm symm}, n,j}^{1-\alpha}$ can then be approximated by the empirical ($1-\alpha$)-quantile of the $m$ draws of $\max_{k\colon k\neq j} |Z_j-Z_k|/\hat{se}_{jk}$.
Under weak conditions, the confidence sets for the differences simultaneously cover all true differences involving country $j$:
Collect the countries $k$ whose differences with $j$ have a confidence set $C_{{\rm symm}, n,j,k}$ that lies entirely below zero, $$N_j^- := \{ k \colon k\neq j\text{ and } C_{{\rm symm}, n,j,k} \subseteq \mathbf R_-\}, $$ and similarly $$N_j^+ := \{ k \colon k\neq j\text{ and } C_{{\rm symm}, n,j,k} \subseteq \mathbf R_+\}. $$ Thus $N_j^-$ contains all countries $k$ that have a significantly larger performance measure than $j$, while $N_j^+$ contains all the countries $k$ that have a significantly smaller performance measure than $j$. If the true performance measures of countries $k$ in $N_j^-$ ($N_j^+$) were indeed all larger (smaller) than that of country $j$, then the rank of country $j$ could not be better than $|N_j^-|+1$ and not be worse than $p-|N_j^+|$. Thus, the set
would contain the true rank of country $j$. Of course, the confidence sets for the differences cover the true differences only with probability approximately no less than $1-\alpha$, so $R_{n,j}$ covers the true rank of country $j$ only with probability approximately no less than $1-\alpha$. In conclusion, for the construction described in this subsection, (ref) implies that $R_{n,j}$ is a confidence set for the rank $R_j^{\theta}$ satisfying (ref) as desired.
It is possible to improve the simple construction of $R_{n,j}$ above by inverting hypotheses tests of
versus its negation, for all $k$ that are not equal to $j$. After testing this family of hypotheses, one then counts the number of hypotheses that were rejected in favor of $\theta_j<\theta_k$ and in favor of $\theta_j>\theta_k$. The first number plus one is then used as lower endpoint and the second number subtracted from $p$ is then used as upper endpoint for $R_{n,j}$. This confidence set satisfies (ref) provided that the procedure used to test the family of hypotheses controls the mixed directional familywise error rate (mdFWER) at $\alpha$, i.e., $$\lim_{n\to\infty} \text{mdFWER} \leq \alpha,$$ where mdFWER is the probability of making any mistake, either a false rejection or an incorrect determination of a sign; see Mogstad:2023aa for details.
A small modification of the above construction of a marginal confidence set for the rank of a single country delivers two-sided confidence sets $R_{n,j}$ for the ranks of all countries $j=1,\ldots,p$ such that all true ranks are covered simultaneously, i.e.,
We start with confidence sets for the differences $C_{{\rm symm}, n,j,k} $ as in (ref) except that the critical value $c_{{\rm symm}, n,j}^{1-\alpha}$ is now defined as the $(1-\alpha)$-quantile of $$\max_{(j,k)\colon k\neq j} \frac{|\hat{\theta}_j - \hat{\theta}_k - (\theta_j-\theta_k)|}{\hat{se}_{jk}},$$ where the max is taken over all pairs $(j,k)$ such that $j\neq k$, so the critical value is independent of $j$. As above this critical value can be approximated by the ($1-\alpha$)-quantile of the $m$ draws of $\max_{(j,k)\colon k\neq j} |Z_j-Z_k|/\hat{se}_{jk}$. Then, the confidence set for country $j$, $R_{n,j}$, is computed as in (ref) using the definitions of $N_j^-$, $N_j^+$ as above except that the confidence sets for the differences, $C_{{\rm symm}, n,j,k} $, are replaced by the new ones described here.
Stepwise methods can be used to improve this simple construction of simultaneous confidence sets similarly to the stepwise improvements described for the marginal confidence sets.
In this section, we are interested in constructing confidence sets for the $\tau$-best countries, defined as $$R_0^{\tau-\rm{best}} := \{j\in\{1,\ldots,p\} : R_j^{\theta} \leq \tau \}.$$ This is the set of countries that have the $\tau$ largest performance measures, the “top-$\tau$” countries. If there are no ties, then this set contains exactly $\tau$ countries. If there are ties, then the set may contain more countries.
We want to learn which countries could be in this set, i.e. among the top-$\tau$. To this end we construct a (random) set $R^{\tau-\rm{best}}_n$ satisfying
Such a confidence set contains all the countries that cannot be rejected to be among the top-$\tau$. By construction, the confidence set contains at least $\tau$ countries, but typically more.
Let $R_{n,j}$, $j=1,\ldots,p$, be simultaneous lower confidence bounds on the ranks of all journals, i.e., each $R_{n,j}$ has upper endpoint equal to $p$ and (ref) is satisfied. Such one-sided confidence sets for the ranks can be constructed similarly as the two-sided confidence sets described in Section (ref), except that the two-sided confidence sets for the differences are replaced by one-sided confidence sets; see Appendix (ref) for details. Then,
is a confidence set satisfying (ref). Mogstad:2023aa propose a different, more direct approach to constructing confidence sets for the $\tau$-best, which in simulations has been shown to produce shorter confidence sets, but is computationally more demanding. The \pkg{csranks} package currently only implements the simpler construction in (ref), which is referred to as the projection confidence set.
Confidence sets for the $\tau$-worst can be constructed as confidence sets for the $\tau$-best among $-\theta(P_1),\ldots,-\theta(P_p)$.
Consider the special case in which the data come from a poll in which a random sample of respondents are asked to choose one of $p$ political parties. Suppose we want to rank the parties according to $\theta_1,\ldots,\theta_p$, the shares of the total population that support them. Let $X_j$ denote the number of times party $j$ has been chosen by respondents in the poll. Then, $X:=(X_1,\ldots,X_p)'$ is distributed according to the multinomial distribution with parameters $n$, the number of respondents, and $\theta:=(\theta_1,\ldots,\theta_p)'$, the vector of multinomial probabilities.
As above the goal is to form marginal and simultaneous confidence sets for the rank of each party. Mogstad:2021bb propose a construction that exploits the multinomial structure of the setup and show that the resulting confidence set $R_{n,j}$ covers the true rank of party $j$ {\it in finite samples}:
for any sample size $n$. This is in stark contrast to the asymptotic coverage property of the more general construction in (ref).
For a given country $j$, the construction of the confidence set is based on hypotheses tests of the family $$H_{k,l}\colon \theta_k\leq \theta_l$$ for all pairs $(k,l)$ such that $k\neq l$ and one of the two is equal to $j$. Then, let $$R_{n,j} := \left\{|\text{Rej}_{j}^-| + 1, \ldots, p-|\text{Rej}_{j}^+|\right\},$$ where
indicate the set of hypotheses that are rejected in favor of $\theta_{j}<\theta_k$ and
the set of hypotheses that are rejected in favor of $\theta_{j}>\theta_k$. If one has a test of this family of hypotheses that controls the familywise error rate (FWER) at $\alpha$, then the resulting confidence set $R_{n,j}$ satisfies $P\{R_j^{\theta} \in R_{n,j}\} \geq 1-FWER \geq 1-\alpha$ as desired. So, the remaining task is to provide tests of the individual hypotheses $H_{k,l}$ such that, overall, the FWER is controlled.
This is achieved in two steps. First, consider testing a single hypothesis for a given pair $k$ and $l$. Let $S_{k,l}:=X_k+X_l$. One can show that the conditional distribution of $X_k$ given $S_{k,l}=s$ is binomial based on $s$ trials and success probability $\theta_k/(\theta_k+\theta_l)$. This is intuitive because, given that parties $j$ and $k$ together have been chosen by $S_{k,l}=s$ respondents, the distribution of $X_k$ is a binomial experiment about how often out of these $s$ trials $k$ was chosen. Therefore, testing a single hypothesis is equivalent to testing an inequality for a binomial probability, for which a test with finite sample validity can easily be constructed (Lehmann:2005p3350). Mogstad:2021bb propose the p-value
for testing the individual hypothesis $H_{k,l}$. In the second step, the individual p-values are then combined using the Holm procedure so as to control the FWER for the family of hypotheses $H_{k,l}$ with $(k,l)$ such that $k\neq l$ and one of the two is equal to $j$.
Confidence sets that simultaneously cover the true ranks for all parties are constructed in a similar fashion except that now one needs to test the family of hypotheses $H_{k,l}$ for all pairs $(k,l)$ such that $k\neq l$.
Slope coefficients in rank-rank regressions are popular measures of intergenerational mobility, for instance in regressions of a child's income rank on their parent's income rank. In this section, we review recent results by Chetverikov:2023aa providing the asymptotic theory for coefficients in regressions involving ranks and the inference methods they propose.
In this section, we are interested in regression models of the form
where, for some $\omega\in[0,1]$, $$R_X(x) := \omega F_X(x) + (1-\omega)F_X^-(x), $$ $F_X$ is the cdf of a random variable $X$, and $F_X^-(x):=P(X<x)$. $R_Y(y)$ is defined analogously based on the same value of $\omega$ as in $R_X(x)$. $W$ is a vector of regressors, $\rho$ and $\beta$ are coefficients of interest.
Suppose we have an i.i.d. sample $\{(Y_i,X_i,W_i)\}_{i=1}^n$ from the distribution of $(Y,X,W)$. In this case, $R_X(x)$ is the probability limit of
and the increasing fractional rank defined in (ref) satisfies $R_i^X = \hat{R}_X(X_i) $. Analogously define $\hat{R}_Y(y)$ so that the increasing fractional rank for $Y$ satisfies $R_i^Y = \hat{R}_Y(Y_i) $. We can then estimate the coefficients $\rho$ and $\beta$ from an OLS regression of $R_i^Y$ on $R_i^X$ and $W_i$:
Chetverikov:2023aa show that, under weak conditions, this estimator is consistent and asymptotically normal,
and derive the expression of $\Sigma$. Importantly, the expression of the asymptotic variance $\Sigma$ differs from the probability limits of commonly used variance estimators such as the homoskedastic and Eicker-White variance estimators that software implementations\footnote{For instance, the \code{lm()} command in \proglang{R} or the \code{regress} command in Stata.} of standard OLS regressions report. This is because both the dependent and the independent variable in the rank-rank regression are estimated. The commonly used variance estimators ignoring this additional estimation error thus lead to invalid standard errors and confidence sets. Chetverikov:2023aa show that, in fact, these standard errors may be too large or too small depending on the shape of the copula of $Y$ and $X$. Therefore, the invalid standard errors may lead to conservative or misleading inference.
In the special case, in which $X_i$ and $Y_i$ are both drawn from a continuous distribution and $W_i$ includes only a constant, the OLS estimator is equal to Spearman's rank correlation. Otherwise, it is not.
The \code{lmranks()} function implements the OLS estimator (ref) and standard errors, p-values and t-values based on a consistent estimator of the correct asymptotic variance $\Sigma$. The estimate of the asymptotic variance $\Sigma$ can be calculated with the \code{vcov()} method applied to an \code{lmranks()} object. It is used internally for other methods for \code{lmranks()} objects, such as \code{summary()} or \code{confint()}.
Considerable care was taken for the implementation of this method to be computationally efficient and scalable. The achieved complexity of the implemented algorithm is linearithmic in terms of the number of observations. Appendix (ref) contains technical details of the implementations.
In this section, we briefly describe some variants of the rank-rank regression which are used in empirical work in economics and also implemented in the \code{lmranks()} function.
\paragraph{Rank-rank regressions with clusters.} We consider a population (e.g., the U.S.) that is divided into $n_G$ subpopulations or “clusters” (e.g., commuting zones). We are interested in running rank-rank regressions separately within each cluster. The ranks, however, are computed from the distribution of the entire population (e.g., the U.S.). Such regressions are common in studies of intergenerational mobility (e.g., Chetty:2018iu), for instance.
Specifically, we consider the model
where $G$ is an observed random variable taking values in $\{1,\ldots,n_G\}$ to indicate the cluster to which an individual belongs. $(G,X,W,Y)$ have distribution $F$ and we continue to denote marginal distributions of $X$ and $Y$ by $F_X$ and $F_Y$. $F_X^-$, $F_Y^-$, $R_X(x)$, and $R_Y(y)$ are also as previously defined, so that $R_X(X)$, for instance, is the rank of $X$ in the entire population, not the rank within a cluster. So, in the model (ref), the coefficients $\rho_g$ and $\beta_g$ are cluster-specific, but the ranks $R_Y(Y)$ and $R_X(X)$ are not. In consequence, $\rho_g$ cannot be interpreted as the rank correlation within the cluster $g$.
Let $\{(Y_i,X_i,W_i,G_i)\}_{i=1}^n$ be a random sample from the distribution of $(Y,X,W,G)$. The coefficients $\rho_g$ and $\beta_g$ for cluster $g$ can be consistently estimated by first constructing the ranks $R_i^X$ and $R_i^Y$ and then running an OLS regression of $R_i^Y$ on $R_i^X$ and $W_i$ using only observations from cluster $g$ (i.e., for which $G_i=g$). Denote by $\hat\rho:=(\hat\rho_1,\ldots,\hat\rho_{n_G})'$ and $\hat\beta:=(\hat\beta_1',\ldots,\hat\beta_{n_G}')'$ the vectors of all cluster-specific OLS estimators of $\rho:=(\rho_1,\ldots,\rho_{n_G})'$ and $\beta:=(\beta_1',\ldots,\beta_{n_G}')'$. Then, Chetverikov:2023aa show that
and derive the expression of the asymptotic variance $\Sigma$.
\paragraph{Regression of a general outcome on a rank.} In this case, we consider a regression model with a general, non-ranked dependent variable $Y$ and a ranked independent variable:
Let $\{(Y_i,X_i,W_i)\}_{i=1}^n$ be an i.i.d. sample from the distribution of $(Y,X,W)$. Chetverikov:2023aa show that the OLS estimator of a regression of $Y_i$ on $R_i^X$ and $W_i$ is asymptotically normal as in (ref) and derive the expression of the asymptotic variance $\Sigma$.
\paragraph{Regression of a rank on a general regressor.} In this case, we consider a regression model with a ranked dependent variable and a general, non-ranked independent variable:
Let $\{(Y_i,X_i,W_i)\}_{i=1}^n$ be an i.i.d. sample from the distribution of $(Y,X,W)$. Chetverikov:2023aa show that the OLS estimator of a regression of $R^Y_i$ on $W_i$ is asymptotically normal, $$\sqrt{n}(\hat\beta-\beta)\to_d N(0,\Sigma) $$ and derive the expression of the asymptotic variance $\Sigma$.
The package \pkg{csranks} comprises 14 functions. These functions implement methods for construction of confidence sets for ranks and inference in rank-rank regressions, described in Section (ref).
The central functions for constructing confidence sets for ranks (as described in Sections (ref)-(ref)) and confidence sets for the $\tau$-best/worst (as described in Section (ref)) are \code{csranks()}, \code{cstaubest()} and \code{cstauworst()}. They all require a vector of estimates \code{x} and an estimate of their covariance matrix \code{Sigma} (if the estimates are independent, the user should pass a diagonal matrix). The nominal coverage of confidence set can be specified with the argument \code{coverage} and the number of bootstrap samples is set with the argument \code{R}. \code{csranks()} also accepts a boolean \code{simul} argument, which specifies whether the returned set should be marginal (\code{FALSE}) or simultaneous (\code{TRUE}) confidence set.
The \code{csranks_multinom()} function is similar to \code{csranks()}, but designed for the special case in which the data is multinomial (as described in Section (ref)). In the multinomial case, the computation of the covariance matrix of the estimates \code{x} requires only knowledge of \code{x} and thus the function does not require the argument \code{Sigma}. The arguments \code{coverage} and \code{simul} play an identical role as in \code{csranks()}. The only new argument, \code{multcorr}, specifies the method used for correction of the p-values for multiple testing (\code{Bonferroni} or \code{Holm}).
For both \code{csranks()} and \code{csranks_multinom()} an S3 \code{plot} method is implemented. Additionally, there are utility functions \code{irank()} and \code{frank()} used to compute integer and fractional ranks, and \code{irank_against()} and \code{frank_against()} to compute ranks of one vector based on values in another reference vector.
The main function for inference in regressions involving ranks (as described in Section (ref)) is \code{lmranks()}. It is designed to be as similar to the well-known \proglang{R} function \code{lm()} (for linear regressions) as possible. As for \code{lm()}, the most important argument of \code{lmranks()} is the \code{formula} argument, which specifies the model using the \proglang{R} formula syntax (RManual). A typical model has the form \code{response terms} where \code{response} is the (numeric) response vector and \code{terms} is a series of terms which specifies a linear predictor for the response. Terms can be added with \code{+}, removed with \code{-}, and a colon \code{:} is used to specify an interaction.
A new functionality in \code{lmranks()} is that the user can specify variables in the formula to be ranked (i.e. their values be replaced by their ranks) before running the regression. This is achieved by wrapping the response or one of the terms with \code{r()}. A typical rank-rank regression model with ranked response \code{Y}, ranked regressor \code{X}, an intercept and non-ranked regressors \code{W1} and \code{W2} is specified with a formula \code{r(Y) r(X) + W1 + W2}.
The \code{weights} argument is not supported due to lack of theory on weighted rank-rank regression and the \code{subset} and \code{na.action} arguments are not supported due to the order in which the \proglang{R} function \code{model.frame} processes the arguments. It first evaluates the \code{formula}, and then applies the \code{subset} and \code{na.action} arguments (RManual). Since those arguments remove observations from the dataset, it affects the distribution of ranks in the resulting dataset and thus is not permitted. The user is therefore required to subset the data and handle the \code{NA} values on their own, before passing data to \code{lmranks()}.
Many functions defined for \code{lm()} also work correctly with \code{lmranks()}. These include \code{coef()}, \code{model.frame()}, \code{model.matrix()}, \code{resid()}, \code{predict()}, \code{update()} and others. On the other hand, some would return incorrect results if they treated \code{lmranks()} output in the same way as \code{lm()}'s and have been disabled. These functions in most cases require the number of degrees of freedom of the model, and for rank-rank regressions it is not yet clear how to calculate them. The central contribution of this package are \code{vcov()}, \code{summary()} and \code{confint()} implementations using the correct asymptotic theory for regressions involving ranks.
Sometimes, the dataset is divided into clusters and one is interested in running rank-rank regressions separately within each cluster, where the ranks are not computed within each cluster, but using all observations pooled across all cluster. This is the model in (ref). This regression model can be written as a rank-rank regression in which the regressors are multiplied by cluster indicators. For $q$ regressors and $n_G$ clusters we get $q n_G$ columns -- one for each regressor-cluster pair. This expansion is conveniently achieved by using functions already implemented in base \proglang{R}. In the \proglang{R} formula, an interaction operator has to be used, and the rest is done in \code{model.matrix}. In \code{lmranks()}, a typical rank-rank regression with clusters specified in the variable \code{G}, ranked response \code{Y}, ranked regressor \code{X}, an intercept and non-ranked regressors \code{W1} and \code{W2} is specified with a formula \code{r(Y) (r(X) + W1 + W2):G}.
The following example illustrates how the \pkg{csranks} package can be used to quantify the statistical uncertainty in the PISA ranking of countries. Over the past two decades, the Organization for Economic Co-operation and Development (OECD) have conducted the PISA study. The goal of this study is to evaluate and compare educational systems across countries by measuring 15-year-old school students’ scholastic performance on math, science, and reading. Each country that participates in a given year draws a sample of students to be tested. The OECD then processes the test results so as to produce a score for each country and publishes league tables ranking countries by their scores.
In this example, we use publicly available data from the 2018 PISA study to examine in which countries school students do best and worst at math.
First, we load the required libraries and the dataset \code{pisa}, which is part of the \pkg{csranks} package. It contain the three test scores and accompanying standard errors for each country (“jurisdiction”):
The PISA study’s math scores are stored in \code{math_score} and their standard errors in \code{math_se}. The following graph shows the raw math scores with 95% marginal confidence intervals:
Figure (ref) shows the resulting graph. The function \code{irank()} can be used to produce integer ranks based on these math scores:
Japan is ranked first (i.e., best), Korea is ranked second and so on. Since the math scores are estimates of countries’ true achievements, the ranks assigned to these countries are also estimates, rather than the true ranks. Just like the test scores, the ranks therefore also contain statistical uncertainty. Various functions in the \pkg{csranks} package implement methods for the quantification of this uncertainty, which were described in section (ref).
Suppose, that the researcher is interested in finding out which countries could be among the top-5 in terms of their true math score. One can answer this question by constructing the $\tau$-best confidence sets, described in section (ref). In \pkg{csranks}, it is implemented in function \code{cstaubest()}. It requires as an argument an estimate of the covariance matrix of the test scores. In this example, it is assumed the estimates from the different countries are mutually independent, so the covariance matrix is diagonal. The function \code{cstaubest()} can then be used to compute a 95% confidence set for the top-5:
The confidence set contains 12 countries: with probability approximately 0.95, these 12 countries could all be among the top-5 according to their true math score. According to the estimated test scores, the countries Japan, Korea, Estonia, Netherlands, and Poland are the top-5 countries. However, due to the statistical uncertainty in the ranking, there is uncertainty about which countries are truly among the top-5.
Suppose that the researcher is interested in a single country, for example the United Kingdom. She would like to learn where its true ranking may lie. A marginal confidence set, described in section (ref) and implemented in the function \code{csranks()} with parameter \code{simul=FALSE}, is a way to answer this question:
\code{CS_marg\$L} and \code{CSmarg\$U} contain the lower and upper bounds of the confidence set for the rank of the United Kingdom.
Based on the estimated math scores, the United Kingdom is ranked at 13-th place. However, due to statistical uncertainty in the ranking, its true rank could be anywhere between 7 and 23, with probability approximately 95%.
Finally, suppose the researcher is interested in the entire ranking of countries. Simultaneous confidence sets (described in Section (ref)) for the ranks of all countries quantify the statistical uncertainty in the entire ranking. They are implemented in the function \code{csranks()} with parameter \code{simul=TRUE}:
The resulting graph is shown in Figure (ref). The simultaneous confidence sets indicate substantial statistical uncertainty about ranks in the middle of the ranking. For instance, the confidence set for the true rank of Germany has a lower bound of 7 and and upper bound of 24. At the top and the bottom of the ranking, the statistical uncertainty is smaller. For instance, with approximately 95% probability, the true rank of Colombia is 37. The confidence sets for Mexico and Chile are also very tight and only contain two values.
The following example illustrates how the \pkg{csranks} package can be used for estimation and inference in rank-rank regressions. These are commonly used for studying intergenerational mobility.
The dataset used in this example is the \code{parent_child_income} dataset that is part of the \pkg{csranks} package. It is a simulated dataset using a data-generating process calibrated to the National Longitudinal Survey of Youth 1979 from the U.S. Bureau of Labor Statistics. It includes data about parents' (column \code{c_faminc}) and children's (\code{p_faminc}) family income, as well as individual characteristics (\code{gender} and \code{race}: \code{"hisp"} (Hispanic), \code{"black"} or \code{"neither"}).
First, we take a quick look at the dataset:
A popular approach to measuring income mobility is to estimate a rank-rank regression of child's income (\code{c_faminc}) on a constant and parent's income (\code{p_faminc}). The \code{lmranks()} function implements this regression. The model can be specified through a formula in which the variables to be ranked are marked by \code{r()}.
This regression specification takes each child's income, computes its rank among all children's incomes, then takes each parent's income and computes its rank among all parents' incomes. Then the child's rank is regressed on the parent's rank using OLS. The \code{summary()} method computes standard errors, t-values and p-values according to the asymptotic theory developed in (Chetverikov:2023aa).
One can also run the rank-rank regression with additional covariates, e.g.:
In some economic applications, it is desired to run rank-rank regressions separately in subgroups of the population, but compute the ranks in the whole population. For instance, we might want to estimate rank-rank regression slopes as measures of intergenerational mobility separately for males and females, but the ranking of children's incomes is formed among all children (rather than form separate rankings for males and females).
Such regressions can be run using the \code{lmranks()} function with interaction notation:
In this example, we have run a separate OLS regression of children’s ranks on parents’ ranks among the female and male children. However, incomes of children are ranked among all children and incomes of parents are ranked among all parents. The standard errors, t-values and p-values are implemented according to the asymptotic theory developed in (Chetverikov:2023aa), where it is shown that the asymptotic distribution of the estimators now need to not only account for the fact that ranks are estimated, but also for the fact that estimators are correlated across gender subgroups because they use the same estimated ranking.
One can also create more granular subgroups by interacting several characteristics such as gender and race:
Finally, we compare the confidence intervals for rank-rank regression coefficients produced by \code{lmranks()} with those of the naive approach which computes the ranks, then runs a regression of the child's income rank on the parent's income rank using \code{lm()}, and then reports the confidence intervals based on the standard errors from \code{lm()}.
Figure (ref) shows the resulting graph comparing the confidence sets obtained from \code{lmranks()} (denoted by “csranks”) with those of the naive use of \code{lm()} (denoted by “naive”) that ignores the estimation error in the ranks. “Whole sample” refers to the confidence sets for intercept and slope when a rank-rank regression is run on the whole sample. The other rows show the confidence sets for the intercept and slope for each group in the grouped rank-rank regression.
The point estimates are the same for both methods, but the confidence intervals of the naive method are not valid. This is because the usual OLS formulas for standard errors do not take into account the estimation uncertainty in the ranks. This leads to different confidence intervals. The differences in the confidence intervals are not particularly large in this dataset, but Chetverikov:2023aa provide further examples in which the differences are considerable.
The results in this paper were obtained using \proglang{R} 4.2.1 with the \pkg{csranks} 1.2.2 and \pkg{ggplot2} 3.4.3 packages. \proglang{R} itself and all packages used are available from the Comprehensive \proglang{R} Archive Network (CRAN) at \url{https://CRAN.R-project.org/}.
The authors gratefully acknowledge financial support from the European Research Council (Starting Grant No. 852332). Chetverikov, Mogstad, Romano, Shaikh, Wilhelm developed the statistical methods presented in this paper. Morgen and Wilhelm developed the software and wrote the article.