# Estimate population quantiles from subpopulations' quantiles

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. • Couple questions: do you know the size of each subpopulation? Also, some quantiles may not be available for some sample sizes due to discreteness. Are you assuming we just have SOME set of quantiles for each partition, but those may be different for each? – user31668 Feb 8 '14 at 1:32
• Yes, we know the size of each subpop. Assume that the subpops are large enough to avoid the discreteness problem. The set of quantiles are the same for each subpop (e.g., all quintiles, all quartiles, all deciles, etc.) – J. Miller Feb 9 '14 at 1:53

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_{} \le x_{} \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_{}=-\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)
lines(c(-3, quantiles), c(0.6, 0.6), col="Gray", lty=2)
lines(c(quantiles,quantiles), 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")

• After thinking through this for a while, you've convinced me that this approach is the best we can do if we want to stay nonparametric. I'm still also playing with the idea of fitting distributions to each subpopulation and then creating a mixture distribution, but in the meantime this approach clearly shows the bounds of the problem. Thanks for the graphs and R snippets too. – J. Miller Feb 13 '14 at 22:25
• That's a good comment: by adopting a parametric model you can do better--sometimes much better--than the results of interval arithmetic, which are (almost by definition) as pessimistic as one can be (and thereby produce the widest possible intervals of uncertainty). – whuber Feb 13 '14 at 22:27

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.

• I greatly appreciate this answer and the effort involved. However, unfortunately I did not state that the subpops were drawn randomly ... I said "arbitrarily". Each subpop may in fact have its own bias. I will clarify this in the question. – J. Miller Feb 10 '14 at 18:28
• @J.Miller in general arbitrary means random or without explicit purpose. Are you saying they were opportunistically collected? If they are not representative of the overall population then they can't be used for inference about your overall population, sicne the "data" are effectively "cherry picked"...imagine trying to make inferences about human height when you only have samples from baseketball teams, bushmen tribes, and a few european cities....you won't really know much unless you know how prevalent each subtype is. – user31668 Feb 10 '14 at 18:32
• "arbitrarily" virtually rules out any specific procedure; it is no information at all. – Nick Cox Feb 10 '14 at 18:33
• @J.Miller see Nick's response as well. You can't do anytihing with truly aribtrary samples. However, any sample will be biased, relative to the true population by virtue of being a sample. What is really of concern is whether or not all the samples are systematically biased, as a group, so that taking more samples will not tell you more about the population. – user31668 Feb 10 '14 at 18:42
• Can't we use the fact that the subpops are partitions? I.e., the sum of their sizes equals the population size -- they completely "cover" the population. – J. Miller Feb 10 '14 at 18:49