Estimate population quantiles from subpopulations' quantiles

Question

Suppose there is a population partitioned arbitrarily into a set of subpopulations that completely cover the original population. Assume that for some variable, we know each subpopulation's quintiles only (or some other fixed set of quantiles), but not its distribution. However, the entire population's percentiles (and of course distribution) are unknown.

Is there a way to estimate the quantiles of the whole population?

This may sound like an odd situation, but it happens for many published economic statistics where access to the full data is not possible, only summaries of it.

EDIT (CLARIFICATION): You may not assume that each subpopulation is drawn randomly. In fact, each is likely to have its own bias.

EDIT2: Additional details... (1) You may assume the variable is nonnegative, and further that we have a decent estimate for its minimum value across the population. (2) Rough or "unrefined" estimates are OK.

FINAL UPDATE: Expanding on @whuber's suggestions and code snippets, I thought I'd add a nice example with some real-world data (hourly wages, if you're curious) that shows the final population upper/lower bounds, along with the known population quantiles, which happened to be available for this case. I've also made some conservative assumptions about the min and max values of each subpopulation, which helps tighten up the intervals a bit.

Answer 1

Because the subpopulations are arbitrary and not random, the best you can do is to use interval arithmetic.

Let the quantiles of a subpopulation be $x_{[1]} \le x_{[2]} \le \cdots \le x_{[m]}$ corresponding to percentiles $100 q_1, 100 q_2, \ldots, 100 q_m,$ respectively. This information means that

1. $100 q_i \%$ of the values are less than or equal to $x_{[i]}$ and

2. $100 q_{i-1} \%$ of the values are less than or equal to $x_{[i-1]},$ whence

3. $100(q_i - q_{i-1}) \%$ of the values lie in the interval $(x_{[i-1]}, x_{[i]}].$

In the case $i=1$, take $q_0 = 0$ and $x_{[0]}=-\infty.$ Similarly take $q_{m+1}=1$ and $x_{[m+1]} = \infty.$

Consider the set of all possible distributions consistent with this information. Let $F$ be the CDF of one of them and suppose $x\in (x_{[i-1]}, x_{[i]}]$ for some $i \in \{0, 1, \ldots, m\}.$ From the preceding information we know

$$q_{i-1} \le F(x) \le q_{i}.$$

The set of all possible CDFs therefore forms a "p-box" filling up these intervals. For example, let the quartiles be $\{-1, 0, 1\}$. The corresponding p-box lies between the upper (red) and lower (blue) curves. A possible distribution $F$ consistent with this p-box is shown in black.

The horizontal gray line shows how quantiles can be read off this plot: the 60th percentile, shown, must lie between $0$ and $1$ given that the 50th percentile is at $0$ and the 75th percentile is at $1$. The solid part of the gray line depicts the interval of possible values of the 60th percentile.

When presented with information of this sort for separate populations of sizes $n_1, n_2, \ldots, n_k$, having associated distributions $F_i, i=1, 2, \ldots, k,$ the distribution for the total population will be the weighted average of the $F_i$:

$$F(x) = \frac{n_1 F_1(x) + n_2 F_2(x) + \cdots + n_k F_k(x)}{n_1 + n_2 + \cdots + n_k}.$$

Because we do not know $F$, we replace it by the p-boxes obtained from the available information and use interval arithmetic to perform the computation. Interval arithmetic in this case is simple: when a value $u$ is known to be in an interval $[u_{-}, u_{+}]$ and $v$ is known to lie in $[v_{-}, v_{+}],$ then certainly $u+v$ is in $[u_{-}+v_{-}, u_{+} + v_{+}]$ and a constant positive multiple $\alpha u$ is in $[\alpha u_{-}, \alpha u_{+}].$ And that's all we can say.

---

For example, suppose we have the following quantile information for subpopulations of sizes $n_i = 5, 4, 7$:

1. Subpopulation 1 has quartiles at $-2, 0, 1$,

2. Subpopulation 2 has quintiles at $-1, 0, 1/2, 3/2,$ and

3. Subpopulation 3 has tertiles at $-4/3, -1/3.$

The resulting p-box computed using interval arithmetic is shown here:

Its 60th percentile (shown by the dashed gray line) must lie between $-4/3$ and $1$, but that is all we know for certain. The distribution of the collective population of $5+4+7=15$ individuals will have a CDF lying somewhere between the upper and lower bounds.

---

R code to compute and manipulate p-boxes is relatively straightforward to write because R supports step functions (the piecewise constant functions that form the envelopes of empirical p-boxes). The hard work is performed by the functions f (which converts quantile specifications into p-boxes) and mix (which forms positive linear combinations of p-boxes).

#
# Create a pair of functions giving the p-box of a set of quantiles.
#
f <- function(quantiles, quants=seq(0, 1, length.out=length(quantiles)+2)) {
  n <- length(quants)
  return (list(lower=stepfun(quantiles, quants[-n]), 
               upper=stepfun(quantiles, quants[-1])))
}
#
# Figure 1: show the p-box for a single population.
#
g <- f(quantiles <- qnorm(c(1,2,3)/4) / qnorm(3/4))
curve(g$upper(x), from=-3, to=2, ylim=c(0,1), n=1001, col="Red", lwd=2,
      ylab="Probability", main="Quartiles {-1, 0, 1}")
curve(g$lower(x), add=TRUE, n=1001, col="Blue", lwd=2)
curve(pnorm(x * qnorm(3/4)), add=TRUE)
lines(c(-3, quantiles[2]), c(0.6, 0.6), col="Gray", lty=2)
lines(c(quantiles[2],quantiles[3]), c(0.6, 0.6), col="Gray")
#
# Figure 2: show how to combine p-boxes using interval arithmetic.
#
quantiles <- list(c(-2, 0, 1), c(-1, 0, 1/2, 3/2), c(-4/3, -1/3))
weights <- c(5, 4, 7); weights <- weights / sum(weights)
mix <- function(x, components, weights) {
    matrix(unlist(lapply(components, function(u) u(x))), ncol=length(weights)) %*% weights
}
g.upper <- lapply(quantiles, function(q) f(q)$upper)
g.lower <- lapply(quantiles, function(q) f(q)$lower)
curve(mix(x, g.upper, weights), from=-5/2, to=2, ylim=c(0,1), 
      ylab="Probability", main="P-box for Three Subpopulations", n=1001,  col="Red")
curve(mix(x, g.lower, weights), add=TRUE, n=1001, col="Blue")
abline(h=0.6, lty=2, col="Gray")

Answer 2

Since your samples collectivey constitiute the entire population, you can re-aggregate to get some ideas about the overall distribution.

Let $N_s$ be the number of subpopulations, $n_i$ be the size of subpopulation $i$,$N_q$ be the number of quantiles, and $q_{ij}$ be the set of quantiles for subpopulation $i$, and let $p_{ij}$ be the associated percentile for that quantile (e.g. the first quartile would have p=0.25) where $p_{kj}=p_{lj}$ .

To get an estimate of the overall population's percentiles, you can aggregate your subpopulation data as follows:

1. Form the set of paired observations $(q_{ij},n_{ij})$ for each subpopulation, where $n_{ij} = p_{ij}n_i$ for each subpopulation. Let $N=\sum n_{i}$
2. Now, you need to form the Kaplan Meier estimate of the cumulative distribution function. Since the Kaplan meier estimator is for right censored data, you will need to multiply your quantiles by -1 to get "reversed quantiles", $r_{ij}$.
3. Group all quantiles together and form an ordered set of reversed quantiles, $r_{(k)}$ where k=1 corresponds to the smallest reversed quantile and k=N corresponds to the largest.
4. You will be estimaing the survival function, S(t), at each $r_{(k)}$. The number "at risk" (from the Kaplan meier estimate) at $r_{(k)}$ will be N minus the sum of the $n_{ij}$ whose reversed quantiles are less than (but NOT equal to) $r_{(k)}$. The number of "deaths" will be the the $n_{ij}$ whose reversed quantile equals $r_{(k)}$.
5. Use the formulas for S(t) in the wikipedia link to estimate the survival function using the above terms.
6. Now, the estimated CDF can be calculated at each actual ordered quantile, $q_{(k)}$: $\hat F(q_{(k)}) = 1 - S(r_{(N-k+1)})$.

This is a nonparametric estimate of the left-censored CDF. It will not be able to give an estimate below the lowest quantile.

A more refined estimate can be made if you are willing to assume a particular distribution for the population (or perhaps a mixture of several distributions). In that case, you can fit a distribution via maximum likelihood:

1. Use the empirial CDF from (6) to guess at the general features of the parametric distribution you want to use (e.g., normal, lognormal, gamma, exponential etc). Let $F(x|\theta)$ be the CDF for your assumed distribution given a set of parameters, $\theta$. For example, the normal distribution has $\theta = (\mu, \sigma^2)$.
2. Form the "censored" likelihood function for your data as follows: $L(\{\# q_{(k)},q_{(k)})\}|\theta) = \prod\limits_{\{(\# q_{(k)},q_{(k)})\}} (F(q_{(k)}|\theta))^{\#q_{(k)}}$. Where the $q_{(k)}$ are the ordered quantiles as before, and $\# q_{(k)}$ is the sum of $n_{ij}$ with $q_{ij}$ less than or equal to $q_{(k)}$
3. Now, you need to maximize L as a function of its parameters. To make the math easier, you can take the logarithm of L to get the LogLikelihood, which turns the product into a sum
4. The set of parameters, $\theta$ that maximizes the likelihood will give you a fitted distribution to your "censored" data. Now, you can estimate any desired percentiles from the fitted distribution.

Those are two possible approaches to your data, depending on how "refined" you need your esimates to be..of couse, with greater precision comes more potential bias if your data do not actually follow your assumed distriubtion.