Model for population density estimation

Question

A database of (population, area, shape) can be used to map population density by assigning a constant value of population/area to each shape (which is a polygon such as a Census block, tract, county, state, whatever). Populations are usually not uniformly distributed within their polygons, however. Dasymetric mapping is the process of refining these density estimates by means of auxiliary data. It is an important problem in the social sciences as this recent review indicates.

Suppose, then, that we have available an auxiliary map of land cover (or any other discrete factor). In the simplest case we can use obviously uninhabitable areas like waterbodies to delineate where the population isn't and, accordingly, assign all the population to the remaining areas. More generally, each Census unit $j$ is carved into $k$ portions having surface areas $x_{ji}$, $i = 1, 2, \ldots, k$. Our dataset is thereby augmented to a list of tuples

$$(y_{j}, x_{j1}, x_{j2}, \ldots, x_{jk})$$

where $y_{j}$ is the population (assumed measured without error) in unit $j$ and--although this is not strictly the case--we may assume every $x_{ji}$ is also exactly measured. In these terms, the objective is to partition each $y_{j}$ into a sum

$$ y_j = z_{j1} + z_{j2} + \cdots + z_{jk} $$

where each $z_{ji} \ge 0$ and $z_{ji}$ estimates the population within unit $j$ residing in land cover class $i$. The estimates need to be unbiased. This partition refines the population density map by assigning the density $z_{ji}/x_{ji}$ to the intersection of the $j^{\text{th}}$ Census polygon and the $i^{\text{th}}$ land cover class.

This problem differs from standard regression settings in salient ways:

1. The partitioning of each $y_{j}$ must be exact.
2. The components of every partition must be non-negative.
3. There is (by assumption) no error in any of the data: all population counts $y_{j}$ and all areas $x_{ji}$ are correct.

There are many approaches to a solution, such as the "intelligent dasymetric mapping" method, but all those I have read about have ad hoc elements and an obvious potential for bias. I am seeking answers that suggest creative, computationally tractable statistical methods. The immediate application concerns a collection of c. $10^{5}$ - $10^{6}$ Census units averaging 40 people apiece (although a sizable fraction have 0 people) and about a dozen land cover classes.

Answer 1

You might want to check work of Mitchel Langford on dasymetric mapping.

He build rasters representing population distribution of Wales and some of his methodological approaches might be useful here.

Update: You might also have a look at work of Jeremy Mennis (especially these two articles).

Answer 2

Interesting question. Here is a tentative stab at approaching this from a statistical angle. Suppose that we come up with a way to assign a population count to each area $x_{ji}$. Denote this relationship as below:

$z_{ji} = f(x_{ji},\beta)$

Clearly, whatever functional form we impose on $f(.)$ will be at best an approximation to the real relationship and thus the need to incorporate error into the above equation. Thus, the above becomes:

$z_{ji} = f(x_{ji},\beta) + \epsilon_{ji}$

where,

$\epsilon_{ji} \sim N(0,\sigma^2)$

The distributional error assumption on the error term is for illustrative purposes. If necessary we can change it as appropriate.

However, we need an exact decomposition of $y_{ji}$. Thus, we need to impose a constraint on the error terms and the function $f(.)$ as below:

$\sum_i{\epsilon_{ji}} = 0$

$\sum_i{f(x_{ji},\beta)} = y_j$

Denote the stacked vector of ${z_{ji}}$ by $z_j$ and the stacked deterministic terms of ${f(x_{ji},\beta)}$ by $f_j$. Thus, we have:

$z_j \sim N(f_j,\sigma^2 I) I({f_j}' e = y_j) I((z_j-f_j)' e = 0)$

where,

$e$ is a vector of ones of appropriate dimension.

The first indicator constraint captures the idea that the sum of the deterministic terms should sum to $y_j$ and the second one captures the idea that the error residuals should sum to 0.

Model selection is trickier as we are decomposing the observed $y_j$ exactly. Perhaps, a way to approach model selection is to choose the model that yields the lowest error variance i.e., the one that yields the lowest estimate of $\sigma^2$.

Edit 1

Thinking some more the above formulation can be simplified as it has more constraints than needed.

$z_{ji} = f(x_{ji},\beta) + \epsilon_{ji}$

where,

$\epsilon_{ji} \sim N(0,\sigma^2)$

Denote the stacked vector of ${z_{ji}}$ by $z_j$ and the stacked deterministic terms of ${f(x_{ji},\beta)}$ by $f_j$. Thus, we have:

$z_j \sim N(f_j,\sigma^2 I) I({z_j}' e = y_j)$

where,

$e$ is a vector of ones of appropriate dimension.

The constraint on $z_j$ ensures an exact decomposition.