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.
79,190 characters · 17 sections · 71 citation commands
Spatial Econometrics for Misaligned Data
Spatial data analysis has become increasingly popular in the social sciences. In many applications, data sets providing the specific location of households, firms, villages, or other economic units are matched by location to data sets with geographic features such as rainfall, temperature, soil quality, ruggedness, or air pollution in order to analyze the impact of such environmental variables on economic outcomes.\footnote{maccini2009under, miguel2004economic, and shah2013drought study the impact of rainfall. dell2014we survey applications using weather data. fabregas2017institutions use measurements of soil nutrients in some locations to make fertilizer recommendations in others. nunn2012ruggedness use terrain ruggedness for identification, and chay2003impact study the impact of air pollution. } Such data underpins important economic research, including policy responses to droughts, smog outbreaks, poor harvests, and other events. A typical issue is that the matched data sets will be misaligned. That is, the respective geographical locations of the observations in the matched data sets do not generally coincide. For instance, a researcher might observe crop outputs from a sample of farms in a large area as well as measurements of rainfall collected over the same area from several weather stations. The locations of the weather stations and the farms will generally not coincide, resulting in misaligned data sets.
The approaches commonly used in social sciences to address the misalignment problem yield inefficient estimates and incorrect confidence intervals, and are consistent only under implausible assumptions. Popular approaches for analyzing such data sets involve imputing a level of the regressor for each misaligned observation of the outcome variable and then proceeding with standard regression analysis. It is common to impute using either the value of the nearest location, the value of the nearest location instrumented with more distant locations maccini2009under, or a distance-weighted average of nearby locations shah2013drought. These methods impute the regressors in an initial step before considering the outcome data. Hence we refer to them as “two-step” or “plug-in” methods. As detailed below and discussed in the simulations and application, two-step methods used in the literature are often inconsistent and are generally inefficient because they do not use all the relevant information.
That being said, efficient estimation is challenging in social science applications because the covariance structure of the regression errors will typically be unknown, which raises a fundamental identification issue since, for instance, the Gaussian maximum likelihood estimator requires stipulating the covariance of the regression errors for point estimation pouliot2016missing.
Specification of the covariance model for the regression errors is arguably the most onerous assumption the researcher needs to make. We may worry that the regression error covariance does not have the nice structure common to geophysical variables, e.g., since we ineluctably omit variables in the regression function. Consider, for example, pseudo-replicates (also referred to as random effects or clustered errors): it may be that some regions and their crops are affected by toxic emissions, but we do not know which regions. Our inability to identify which observations share this regional shock prevents us from identifying the correct covariance model. This motivates an estimation scheme relieved of that specification burden.
Indeed, the researcher may or may not have a suitable model for the (spatial) covariance of regression errors. We thus require methodology that will dispense with the requirement to model the regression error covariance structure while preserving as much statistical power as possible.
We propose an easy-to-implement two-step method, based on a Krig-and-regress approach, that does not require stipulation of the regression error covariance and produces standard errors that account for uncertainty due to first stage estimation.
We further argue that a suitable choice of moments, excluding the regression error covariance, identifies the regression coefficients, and we construct a minimum-distance estimator based on these moments. Jointly estimating all coefficients, this estimator brings about efficiency gains.
These two methods should be thought of as complementary. As further detailed below, while the two-step Krig-and-regress estimator is substantially simpler and more transparent, efficiency gains may obtain from the minimum-distance estimator, especially on gridded data.
Reproducing and extending the cross-validation exercise of madsen2008regression, we find that both estimators are competitive, but the more efficient minimum-distance estimator outperforms the Krig-and-regress estimator when the locations are regular, i.e., a subset of a lattice.
We reanalyze the influential data set of maccini2009under and find that the analysis benefits from the use of our methods, as those yield some statistically and economically significant changes in the value of key parameter estimates.
We are interested in the regression coefficient $\beta$ in the spatial regression problem
where $\mathbf{Y}=Y(x)$ is an $N$-tuple $(Y(x_1), \dots, Y(x_N))^T$ and
is drawn from a stationary random field\footnote{We further speak of $\left\{ R(x):x\in\mathbb{R}^{d}\right\} $ as a Gaussian random field if, for any choice of vector of locations $(x_{1},...,x_{n})$, the random vector $(R(x_{1}),...,R(x_{n}))$ is distributed multivariate normal. The practical usefulness of Gaussian random fields to model rainfall data has long been established, see in particular phillips1992comparison and tabios1985comparative. Nowhere will we, however, require the correct specification of Gaussianity.} $R(\cdot)$ evaluated at geographic locations $x_{i}\in\mathscr{\mathscr{D}}\subset\mathbb{R}^{2}$, $i=1,...,N$. The error vector $\epsilon \in \mathbb{R}^{N}$ has mean zero and unknown covariance matrix $\Sigma$ and is independent of $R$. The outcome vector $\mathbf{Y}$ and the matrix of controls $F$ are observed, and the vector of regression coefficients of the controls $\gamma$ is unknown. The matrix of control variables $F$ may include covariates such as age, location dummies, etc. The difficulty is that $\mathbf{R}_{\mathrm{true}}$ is not observed. However, the $M$-tuple
with $x_{i}^{*}\in\mathscr{\mathscr{D}}\subset\mathbb{R}^{2}$, $i=1,...,M$, is observed. That is, although the outcome variable data $Y(x)$ (e.g., crop yields at farm locations) is not sampled at the same locations as the independent variable data $R(x^{*})$ (e.g., rain measured at fixed weather stations), it is $R$ evaluated at the same locations as that of the outcome variable, that is $R(x)$, which enters the regression function.
The random field $R(\cdot)$ has mean function $m:\mathscr{\mathscr{D}}\rightarrow\mathbb{R}$ and covariance function $K_{\theta}:\mathscr{\mathscr{D}}\times\mathscr{\mathscr{D}}\rightarrow\mathbb{R}_{+}$, where $\theta$ indexes a parametric model for the covariance function. The mean and variance of the observables are then
and
respectively, where $\mathbf{K}=K_{\theta}(x,x)=V_{\theta}(R(x))\in \mathbb{R}^{N \times N}$, $\mathbf{\bar{K}}=K_{\theta}(x,x^{*})=Cov_{\theta}(R(x),R(x^{*}))\in \mathbb{R}^{N \times M}$ and $\mathbf{K^{*}}=K_{\theta}(x^{*},x^{*})=V_{\theta}(R(x^{*}))\in \mathbb{R}^{M \times M}$ for some $\theta\in\Theta$.
In the absence of rainfall measurements at the locations of outcomes, the identifying assumption is that we have a parametrized covariance function, and know $Cov_{\theta}(\mathbf{R}_{\mathrm{true}},\mathbf{R}^{*})$ up to the value of a small-dimensional parameter vector $\theta$, which we estimate consistently. This allows, for instance, the construction of a best linear unbiased predictor for unobserved rainfall.
For our purposes, it will generally be the case that $m$ is constant with respect to the location $x$, and thus the mean parameter of the random field $R$ can be absorbed in the constant vector (for the intercept) in the matrix of control variables $F$. Hence, we are concerned throughout with the mean and variance
and
where the coefficient of interest, $\beta$, only appears in the covariance.
This article pertains to different segments of the literature. First and foremost, we provide methods for applied researchers. As detailed above, even careful applied work maccini2009under, shah2013drought relies on ad hoc approaches because econometric methodology has not caught up to the needs of applied economists. This paper intends to address those methodological needs. It speaks to a well-established literature on generated regressors pagan1984econometric, murphy1985estimation, and looks at problems in which the regressors to impute are amenable to imputation using the best linear predictor, e.g., their law well approximated by a Gaussian random field. Gaussian random fields and best linear prediction, also known as Kriging, are well-studied in geostatistics gelfand2010handbook, stein2012interpolation where, however, interest is concentrated on interpolation. This paper relates the two literatures and leverages geostatistical methodology and results to provide robust and accurate methods for economists carrying out regression analysis with misaligned data.
The central limit theorem we obtain for inference with the minimum-distance estimator builds on results in spatial statistics. lahiri2002asymptotic give a central limit theorem for empirical variogram estimators when the data is on a lattice. We leverage results from lahiri2003central, who gives a family of central limit theorems for spatial statistics, in order to extend the asymptotic theory for the empirical variogram estimators to the case of irregularly spaced data -- which is increasingly common in development economics as well as other fields of applied economics.
An alternative asymptotic theory is laid out in the work of jenish2009central, jenish2012spatial. Their theory is very general, even allowing for nonstationarity. Definitions of “spatial mixing” and rate assumptions are not critically different. As suggested by the application to minimum-distance estimation on lattice data lahiri2002asymptotic, the general results in lahiri2003central are immediately relevant to the analysis of our minimum-distance estimator. That being said, we have found, for instance, that the approach of jenish2012spatial is more convenient for developing maximum likelihood asymptotic distribution theory, as exemplified in the quasi-likelihood case by qu2017qml and the well-specified case by xu2015maximum.
The remainder of the article is divided as follows. Section (ref) presents and discusses key concepts for the analysis. Section (ref) presents Krig-and-regress and minimum-distance methods dealing with misaligned data. Section (ref) studies the comparative performance of the considered estimators in the cross-validation exercise of madsen2008regression. Section (ref) reanalyzes the misaligned data set in maccini2009under. Section (ref) discusses and concludes. Technical material is deferred to the Appendix.
Two-step methods for regression analysis of misaligned data consist in first predicting the misaligned covariates at the outcome locations where they are not observed, thereby generating an aligned data set, and then proceeding to spatial regression analysis with this generated data set. The first step, which consists of predicting the missing independent variables, requires the choice of an interpolation method. Nonparametric methods, such as approximation by the average of a given number of nearest neighbors, may be used. However, when the misaligned variable can be modeled as following, or approximately following, the law of a Gaussian random field, Kriging generally affords the researcher more accurate interpolation gelfand2010handbook.
Kriging (named after the South African mining engineer D. G. Krige) consists in using the estimated best linear unbiased predictor for interpolation. It can be developed as follows stein2012interpolation. The random field of interest, $R$, is assumed to follow the model \[ R(x)=s(x)^{T}\rho+\varepsilon(x), \] $x\in\mathscr{D}\subset\mathbb{R}^{2}$, where $\varepsilon$ is a mean zero random field, $s$ is a known function with values in $\mathbb{R}^{p}$ and $\rho$ is a vector of $p$ coefficients. We observe $R$ at locations $x_{1}^{*},x_{2}^{*},...,x_{M}^{*}$. That is, we observe $\mathbf{R}^{*}=\left(R(x_{1}^{*}),R(x_{2}^{*}),...,R(x_{M}^{*})\right)$ and need to predict $R(x_{0})$. With $\rho$ known, the best linear predictor (BLP) is \[ s(x_{0})^{T}\rho+\mathbf{k}^{T}\mathbf{K}^{-1}(\mathbf{R}^{*}-\mathbf{S}\rho), \] where $\mathbf{k}=Cov(\mathbf{R}^{*},R(x_{0}))$, $\mathbf{K}=Cov(\mathbf{R}^{*},\mathbf{R}^{*T})$ and $\mathbf{S}=(s(x_{1}^{*}),s(x_{2}^{*}),...,s(x_{M}^{*}))^{T}$. Of course, the mean parameter $\rho$ is, in general, unknown. If $\rho$ is replaced by its generalized least-squares estimator, $\hat{\rho}=\left(\mathbf{S}^{T}\mathbf{K}^{-1}\mathbf{S}\right)^{-1}\mathbf{S}^{T}\mathbf{K^{-1}}\mathbf{R}^{*}$ (under the assumption that $\mathbf{K}$ and $\mathbf{S}$ are of full rank), we obtain the best linear unbiased predictor (BLUP) for $R(x_{0})$. Again, in general, the covariance structure will be unknown, and $\mathbf{k}$ and $\mathbf{K}$ will be replaced by estimates $\hat{\mathbf{k}}$ and $\hat{\mathbf{K}}$. The resulting plug-in estimator will be called the estimated BLUP (EBLUP). Prediction with the BLUP and EBLUP are both referred to as Kriging. As far as this article is concerned, the covariance structures will always be a priori unknown, and Kriging will refer to prediction with the EBLUP.
There are many choices for the covariance functions gelfand2010handbook, and we present three of them in the isotropic case in which only the distance $d$ between the covariates determines their covariance: the exponential covariance function
the Gaussian covariance function
and the Mat\'ern covariance function
where $\mathcal{K}_{\nu}$ is the modified Bessel function of the second kind of order $\nu$ abramowitz1964handbook. All functions have positive parameters $\theta_{1}$ and $\theta_{2}$, which are the sill and range, respectively. The sill parameter should be thought of as controlling the scale of the covariance, and the range should be thought of as controlling how fast the covariance decays over distance. The Mat\'ern function has an additional parameter $\nu$, which controls smoothness.
jiang1997derivation and stein2012interpolation present an alternative derivation of the BLUP as the best predictor, under normality, based on all error contrasts. An excellent theoretical treatment of the topic can be found in stein2012interpolation. cressie2015statistics and diggle2007springer offer a more applied treatment of the topic. matheron1962traite is a classic reference.
Kriging, or the best linear prediction of missing variables, naturally extends regression, or the estimation of best linear prediction of the outcome variable, to the misaligned case. Indeed, we want to estimate the best linear predictor\footnote{The best linear predictor is defined analogously to the conditional expectation. We may define the conditional expectation as the best predictor $f(R^{*})$ minimizing $E\left[\left(R-f(\mathbf{R}^{*})\right)^{2}\right]$. Likewise, the best linear predictor is the linear function $\alpha^{T}\mathbf{R}^{*}$ minimizing $E\left[\left(R-\alpha^{T}\mathbf{R}^{*}\right)^{2}\right]$ over all $\alpha\in\mathbb{R}^{M}$. } $E^{*}\left[Y\left|R\right.\right]=\beta_{0}+\beta R$ (omitting all other control variables for simplicity of exposition) but only observe $Y$ and $R^{*}$. However, note that
hence a consistent estimate of $E^{*}\left[R\left|R^{*}\right.\right]$ will deliver a consistent estimate of the best linear predictor $E^{*}\left[Y\left|R^{*}\right.\right]$ and, in particular, of $\beta$.
To be sure, although two-step methods commonly used in the literature (e.g. nearest neighbors, two-stage least squares) are not in general consistent, Krig-and-regress is a consistent two-step method. However, Krig-and-regress does not make efficient use of the data; in the first stage, only the covariance $Cov\left(\mathbf{R}^{*},\mathbf{R}^{*}\right)$ is used to estimate the rainfall covariance $\theta$ and build the EBLUP, but $\theta$ is also informed by the second stage statistic $Cov\left(\mathbf{Y},\mathbf{R}^{*}\right)=\beta\mathbf{\bar{K}}$. Minimum-distance estimation offers a natural alternative and uses both moments to jointly estimate $\theta$ and $\beta$.
One of our main motivations is to build principled methods for misaligned regression that are pivotal with respect to $\Sigma$, the covariance of the regression errors. If one wants to circumvent altogether the need to model the covariance of the regression errors, the two-step bootstrap of Subsection (ref) delivers consistent estimation and valid inference. If one has a model for the covariance of the regression errors but does not want to rely on it for point estimation, the minimum-distance estimator may be employed. With the latter estimator, covariance parameters appearing in both the first and second stage problems are estimated only once, jointly on the first and second stage data, and specification of the covariance of regression errors, $\Sigma$, is only required for inference. If one wants to avoid specifying $\Sigma$ even for inference with the minimum-distance estimator, a heuristic approximate Bayesian computation (ABC) approach to inference is provided which does not require stipulating $\Sigma$. Table (ref) summarizes when each method should be preconized. The R package SpReg implementing both methods is available on the author's webpage.\footnote{https://sites.google.com/site/guillaumeallairepouliot/research}
The main motivation for resorting to a two-step method such as Krig-and-regress is the desire to avoid specifying a model for the covariance of the regression errors, $\Sigma$. Indeed, the researcher may conclude that modeling the covariance of the regression errors is too restrictive. For instance, specification of different but equally credible covariance structures may yield tangibly different maximum likelihood estimators pouliot2016missing. Furthermore, researchers may preconize the two-step method for its simplicity and ease of implementation.
Explicitly, the Krig-and-regress point estimation method proceeds in two steps:
Correct inference requires standard errors that take into account the uncertainty brought about by the imputation of the missing covariates. That is, one needs a two-step estimator for which neither estimation nor inference requires knowledge of $\Sigma$ and whose standard errors account for the uncertainty in the first step.
Versions of this problem have come up in the literature under many guises pagan1984econometric. A very general case is addressed by murphy1985estimation who provide an asymptotic covariance formula with a positive-definite correction term accounting for the variation due to the estimation of the imputed regressors. However, such two-step standard errors again require the stipulation of the covariance of the regression errors, $\Sigma$.\footnote{There is a delicate conceptual point here. In the frequentist framework with a DGP corresponding to (ref), the variance of the Krig-and-regress estimate $\hat{\beta}_{KR}$, obtained by generating regressors $\hat{R}^{*}$ via Kriging and then regressing $Y$ on $\hat{R}^{*}$, depends on $\Sigma$. Consequently, we must entertain a different but still credible data generating process (DGP) in order to produce pivotal inference methodology.}
maccini2009under circumvent this issue by relying on a two-stage least-squares approach. However, identification under this approach is hard to argue, and estimates of the regression coefficients need not be consistent (see Subsection (ref)).
madsen2008regression work out standard errors for the regression coefficient estimated using the Krig-and-regress method. They provide a protocol for estimating the unconditional (on the realization of the random field for the misaligned regressor) variance of the regression coefficient. They find, in their application, that the produced unconditional standard errors differ only mildly from the OLS standard errors. Crucially, the standard errors they provide do not account for the estimation of the covariance parameter for rainfall, i.e., the uncertainty due to first-stage estimation.
Our concern, in contrast to that of madsen2008regression, is to provide confidence intervals that take into account the uncertainty due to the estimation of the imputed regressor. Note that in the Krig-and-regress method, accounting for the variation due to the estimation of the imputed regressor is tantamount to accounting for the variance due to the estimation of the covariance and mean parameters of the random field of the misaligned regressor (say, rainfall).
Since the motivation for using a two-step method is to avoid the modeling burden of specifying the covariance structure for the regression errors, we ought to produce standard errors that do not require evaluating $\Sigma$.
It is plausible that the residual errors of the best linear predictor \[ Y(X)-E^{*}[Y(X)|R(X)], \] where $E^{*}$ is the best linear prediction operator, are spatially correlated (for instance, through an omitted variable such as pollution, which concentrates differently in different areas). Uncertainty assessments of the Krig-and-regress coefficient estimates, if we do not condition on the realization of $R$ and $Y$, are thus bound to rely on the estimation of $\Sigma$.\footnote{Unless the estimation is done with standard errors so conservative that they are useless in practice.} A different yet plausible stochastic model must thus be entertained.
Consider the identification strategy in the context of our main application. Let $\mathscr{D}\subset\mathbb{R}^{2}$ be the geographic domain under study. Let $R$ and $Y$ be the random fields for rainfall and the outcome variable from a geolocated household survey (say, height), respectively. Let $X^{*}\subset\mathscr{D}$ be the locations of the rainfall measurement stations. Let $\hat{R}=E^{*}[R|R^{*}]$, where $E^{*}$ is the best linear prediction operator and $R^{*}=R(X^{*})$ is the observed rainfall. Let $X\subset\mathscr{D}$ be the locations of surveyed households.
The key observation is that if all the uncertainty in the second step arises from the resampling (with replacement) of the surveyed households, then the observations $(Y(X_{i}),\hat{R}(X_{i}))$ are independently distributed, conditional on the realization of $R$ and $Y$. That is, conditional on the realization of the random field of rainfall and the outcome variable at all households (but unconditional on which household is randomly drawn with replacement from the population), the observations are independent and identically distributed.
Remarkably, by modeling the sources of randomness as detailed in the previous paragraph, one can provide correct inference without relying on knowledge of $\Sigma$. This is, of course, in contrast to standard frequentist inference, which is unconditional on the observed data and where uncertainty in the regression coefficients is thought to capture the variation arising from repeated samples of $(Y,R,X)$.\footnote{For simplicity of exposition, we omit the additional covariates in this section.}
Indeed, we can treat the problem as one of survey sampling and estimate the linear regression coefficient of $Y(\chi_{\mathrm{hh}})$ on $R(\chi_{\mathrm{hh}})$ over the whole population (say, of Indonesian households) for the given realizations of the random fields $Y$ and $R$, where $\mathcal{X}_{\mathrm{hh}}$ is the set of locations of all households.
If the surveyed households are drawn with replacement from the full population, then the corresponding observations $(Y_{i},R_{i})$ sampled with replacement from $(Y(\chi_{\mathrm{hh}}),R(\chi_{\mathrm{hh}}))$ will be independent and identically distributed.
Crucially, this modeling approach can furthermore accommodate inference which takes into account first-stage uncertainty. Indeed, variance due to estimation in the first stage is captured by allowing the locations of the rainfall stations to be random, while still conditioning on the realization of the random field of rainfall.
We need to be clear as to what kind of process we have in mind to accommodate the misaligned case. The variation in the outcome variable comes from the random selection of the survey households, whose locations are collected in $\mathbf{X}_{\mathrm{hh}}\subset\mathcal{X}_{\mathrm{hh}}$. This is a natural assumption as it mimics the original data collection process of the main application maccini2009under, which was itself a survey. We assume a similar survey sampling scheme of the locations $\mathbf{X}_{\mathrm{rain}}$ of the rainfall stations from the set of all possible rainfall locations, $\chi_{\mathrm{rain}}$. We consider that $\chi_{\mathrm{rain}}$ contains all the locations where the weather stations could have been and that their locations collected in $\mathbf{X}_{\mathrm{rain}}$ were selected independently and uniformly at random from $\chi_{\mathrm{rain}}$; i.e., we condition on $R(\chi_{\mathrm{rain}})$ and it is the location of the weather stations which is different if we resample a new data set from the DGP.\footnote{Sampling of the weather station locations need not be uniform or over a finite set; pivotal inference obtains as long as we condition on $R(\cdot)$, and locations may be sampled according to any distribution $f$ over the domain.}
The key point is that, since the sampling of the locations of the realizations of $Y$ and $R$ are independent, the outcome data variation does not inform the variation of the interpolated values.
The exact nature of the target parameters depends on our conceptualization of the DGP. We may still believe that the data arose from a law under which the best linear predictor in population has the form $E^{*}[Y|R]=\beta_{0}+R\beta$, and carry out inference conditional on the realization of the random fields $Y(\cdot)$ and $R(\cdot)$. Of course, the target, conditional on the realizations of the random fields, is the regression coefficient of $Y(\chi_{\mathrm{hh}})$ on $R(\chi_{\mathrm{hh}})$, which we call $\beta_{\chi_{\mathrm{hh}}}$.
However, under large domain asymptotics and standard ergodicity assumptions, we consistently estimate the population best linear predictor coefficient, $\beta$. The intuition for why conditional estimation delivers a consistent estimate of an unconditional quantity is straightforward. We can think of chopping off from an arbitrarily large random field an increasing number of themselves increasingly large and increasingly pairwise distant random fields; these increasingly many random fields will be asymptotically independent from each other, thereby delivering independent replications of random fields drawn from the underlying DGP. Therefore, estimated coefficients will benefit from the consistency properties that obtain under usual frequentist asymptotics.
We believe the assumption of resampling with replacement is innocuous. Certainly, sampling without replacement describes more accurately the sampling protocol of the survey. Nevertheless, the survey size is so small compared to the population size that both sampling methods (with and without replacement) yield the same observables with high probability.
We detail the estimator which consistently estimates the population regression coefficient and does not require specification of $\Sigma$ for neither point estimation or inference.
Given $\theta$, and conditional on the realization of the random field of rainfall as well as the outcome variable for each household of the population, $\hat{\beta}$ only depends on which households are drawn (randomly, with replacement) to be part of the survey. This variation is captured by resampling using the bootstrap. This naturally suggests a two-step bootstrap procedure in which $\theta$ is first drawn, accounting for variation in $R(\mathbf{X}_{\mathrm{rain}})$, to determine $\hat{R}(\theta)$, thus capturing the uncertainty due to the estimation of $\theta$.
Instead of bootstrapping the rainfall data and estimating multiple times a maximum likelihood estimate for the coefficients of the covariance function of rainfall, we rely on the heuristic large sample argument described above (applied to the first step), and use the much more convenient asymptotic distribution of the maximum likelihood estimator $\hat{\theta}_{\mathrm{mle}}$, obtained from training only on rainfall data. The full procedure is described in pseudocode as follows:
For each $j=1,...,J,$
Quantiles from the set of bootstrap draws $\{\beta^{(j)}\}_{j=1,...,J}$ can be used to form confidence intervals, and the average can be used to give a point estimate. See Section (ref) of the Appendix for the case of random fields $R(\cdot)$ with mean component $s(\cdot)^{T}\rho$.
A simple identification argument naturally suggests a procedure delivering pivotal point estimation. Upon inspection of (ref), we find that $\beta$ is identified from $V_{R^{*}}=\mathbf{K}^{*}$ and $V_{YR^{*}}=\beta\bar{\mathbf{K}}^{T}$ alone, and hence identification does not require modeling the covariance structure of the regression errors. Specifically, using rainfall data only, we obtain an estimate of the covariance parameter $\theta$ and thus obtain an estimate of $\bar{\mathbf{K}}$. Since $\beta\bar{\mathbf{K}}$ is directly identified from the covariance between the outcome data and the observed rain, $\beta$ is identified. This naturally invites a procedure that will rely on this identification observation to produce a robust estimate of $\beta$.
We develop a minimum-distance estimator. As opposed to two-step methods, all parameters showing up both in the first and second stage problems are estimated once, jointly on both first and second stage data. In order to conduct inference, we develop limit distribution theory, which is detailed in Appendix (ref). The resulting asymptotic covariance matrix depends on the covariance of the regression errors, which may be unpalatable to some users. In order to carry out pivotal inference, we make a novel use of likelihood-free inference methods.
We state the minimum-distance estimator and its limit distribution in terms of variograms to dovetail the results of lahiri2002asymptotic. Note that the minimum-distance estimator and theory articulated in terms of covariances instead of variances obtains analogously (see, for instance, Subsection (ref)).
Let $\gamma_{R^{*}}(h;\phi)=V_{\phi}(R(x)-R(x+h))$ and $\gamma_{YR^{*}}(h;\phi)=V_{\phi}(R(x)-Y(x+h))$, where $\phi=(\beta,\theta)$, be the variogram of $R^{*}$ and the covariogram of $Y$ with $R^{*}$, respectively. Note that $\gamma_{YR^{*}}(h;\phi)=(1+\beta^{2})V_{\phi}\left(R(x)\right)-2\beta Cov_{\phi}\left(R(x),R(x+h)\right)$. Let $\hat{\gamma}_{R^{*}}(h)$ and $\hat{\gamma}_{YR^{*}}(h)$ be nonparametric estimators of $\gamma_{R^{*}}(h;\phi)$ and $\gamma_{YR^{*}}(h;\phi)$, respectively. Note that the method defined in terms of covariances can be implemented analogously.
Let $\{h_{1},...,h_{K_{R^{*}}}\}$ and $\{c_{1},...,c_{K_{YR^{*}}}\}$ be finite sets of lag vectors in $\mathbb{R}^{2}$ such that $\hat{\gamma}_{R^{*}}(h_{i})$ is defined for all $i=1,...,K_{R^{*}}$ and $\hat{\gamma}_{YR^{*}}(c_{j})$ is defined for all $j=1,...,K_{YR^{*}}$.\footnote{The lags can be the default lags of the directional empirical variogram estimator from a geostatistical package, such as $\mathrm{\mathtt{gstat}}$ in R.} Let \[ g_{n}(\phi)=\left(\hat{\gamma}_{YR^{*}}(c_{1})-\gamma_{YR^{*}}(c_{1};\phi),\dots,\hat{\gamma}_{YR^{*}}(c_{K_{YR^{*}}})-\gamma_{YR^{*}}(c_{K_{YR^{*}}};\phi),\right. \] \[ \left.\hat{\gamma}_{R^{*}}(h_{1})-\gamma_{R^{*}}(h_{1};\phi),\dots,\hat{\gamma}_{R^{*}}(h_{K_{R^{*}}})-\gamma_{R^{*}}(h_{K_{R^{*}}};\phi)\right)^{T}. \] For some positive-definite weighting matrix $B_{n}$, define the minimum-distance estimator
\[ \hat{\phi}_{\text{M-D}}=\arg\min_{\phi\in\Phi}g_{n}(\phi)^{T}B_{n}g_{n}(\phi), \] for a convex support $\Phi=\mathbb{R}\times\Theta$. Then $\hat{\beta}_{\text{M-D}}$, the estimate of $\beta$, does not depend on any specification of the covariance structure of the regression errors. Different choices of $B_{n}$ will correspond to different traditional estimators; $B_{n}=I$ yields the ordinary least-squares estimator, $B_{n}=\mathrm{diag}(b_{n,1}(\phi),...,b_{n,K_{YR^{*}}+K_{R^{*}}}(\phi))$, for some choice of weights $b_{n,i}$, $i=1,...,K$, gives the weighted least squares estimator, and $B_{n}(\phi)=\Sigma_{g}^{-1}(\phi)$, where $\Sigma_{g}(\phi)$ is the asymptotic covariance matrix of $g_{n}(\phi)$ is the generalized least-square version of the minimum-distance estimator. Our suggested rule of thumb is to use the efficient $B_{n}(\phi)=\Sigma_{g}^{-1}(\phi)$ evaluated at the Krig-and-regress estimates.
Another attractive feature of this estimator is its flexibility. The vector of moments can be extended to accommodate other conditions, perhaps motivated by economic theory.
In order to carry out inference using the proposed minimum-distance estimator, we need asymptotic distribution theory for the statistic, which is the empirical variogram (defined below). lahiri2002asymptotic develop such theory for data on a regular lattice, and give the asymptotic distribution of the minimum-distance estimator as a corollary. lahiri2003central proves a series of useful central limit theorems for spatial statistics, some of which can be leveraged to extend the asymptotic theory for the empirical variogram to the case of irregular data. In contrast to the inferential framework of Section (ref), the inference here does not rely on survey sampling and is unconditional on the realization of the random fields. Intuitively, the asymptotic framework considers the realization of a single “arbitrarily large” random field, and sampling variation obtains because arbitrarily distant observations are effectively independent.
To approximate the variogram $E\left[\left(\varepsilon(x)-\varepsilon(x+h)\right)^{2}\right]$ of a given random field $\varepsilon$, we define the empirical variogram gelfand2010handbook, \[ \hat{\gamma}(h)=\frac{1}{\left|N_{n}(h)\right|}\sum_{(s_{i},s_{j})\in N_{n}(h)}\left(\hat{\varepsilon}(s_{i})-\hat{\varepsilon}(s_{j})\right)^{2}, \] where $\hat{\varepsilon}(s)$ is an estimate of the random component $\varepsilon(s)$ and the bin $N_{n}(h)$ is the set of pairs of observations separated by a vector close to $h\in\mathbb{R}^{2}$. If, instead, the summation is over $\left(\hat{\varepsilon}(s_{i})-\hat{\varepsilon}'(s_{j})\right)^{2}$ for distinct random fields $\varepsilon$ and $\varepsilon'$, then we speak of the covariogram between both random fields and of its empirical estimator.
The limit distribution theory for $g_{n}$ and $\hat{\phi}_{\text{M-D}}$ can be obtained in the pure- and mixed-increasing domains with the so-called stochastic design lahiri2003central. Explicitly, the sampling region, denoted $\mathcal{R}_{n}$, is for each $n$ a multiple of a prototype region $\mathcal{R}_{0}$, defined as follows. The prototype region satisfies $\mathcal{R}_{0}^{*}\subset\mathcal{R}_{0}\subset\bar{\mathcal{R}}_{0}^{*}$, where $\mathcal{R}_{0}^{*}$ is an open connected subset of $(-1/2,1/2]^{2}$ containing the origin. Let $\left\{ \lambda_{n}\right\} _{n\in\mathbb{N}}$ be a sequence of positive real numbers such that $n^{\epsilon}/\lambda_{n}\rightarrow0$ as $n\rightarrow\infty$ for some $\epsilon>0$. Then the sampling region is defined as \[ \mathcal{R}_{n}=\lambda_{n}\mathcal{R}_{0}. \] To avoid pathological cases, we will assume that the boundary of $\mathcal{R}_{0}$ is delineated by a smooth function. This assumption can be modified and weakened to adapt to other domains (e.g., star-shaped), see lahiri2003central.
Furthermore, we speak of a stochastic design because the data is not placed on a regular lattice, and observation locations must be modeled otherwise. They are modeled as follows. Let $f(x)$ be a continuous, everywhere positive density on $\mathcal{R}_{0}$, and let $\{X_{n}\}_{n}$ be a sequence of independent and identically distributed draws from $f$. Let $x_{1},...,x_{n}$ be realizations of $X_{1},...,X_{n}$, and define the locations $s_{1},...,s_{n}$ of the observed data in $\mathcal{R}_{n}$ as \[ s_{i}=\lambda_{n}x_{i},\ i=1,...,n. \]
In the stochastic design, pure-increasing asymptotics require that $n/\lambda_{n}^{2}\rightarrow C$ for some $C\in(0,\infty)$ as $n\rightarrow\infty$. Mixed-increasing asymptotics require that $n/\lambda_{n}^{2}\rightarrow\infty$.
First, we obtain a central limit theorem for the statistic entering the minimum-distance objective function. For simplicity of exposition, take $g_{n}$ as defined above but in which enter as statistics only variograms in $R$, and let $\{h_{1},...,h_{K_{R^{*}}}\}$, $K_{R^{*}}\in\mathbb{N}$, be the full set of lag vectors. The result is trivially extended to accommodate as statistics the covariogram in $R$ and $Y$.
Define the mixing coefficient \[ \alpha(a;b)=\sup\left\{ \tilde{\alpha}(T_{1},T_{2}):d(T_{1},T_{2})\geq a,\ T_{1},T_{2}\in\mathcal{S}_{3}(b)\right\} , \] where $\mathcal{S}_{3}(b)=\left\{ \bigcup_{i=1}^{3}D_{i}:\sum_{i=1}^{3}\left|D_{i}\right|\le b\right\} $ is a collection of disjoint unions of three cubes $D_{1},D_{2},D_{3}$ in $\mathbb{R}^{2}$, $d(T_{1},T_{2})=\min\left\{ \left\Vert x_{1}-x_{2}\right\Vert :x_{1}\in T_{1},x_{2}\in T_{2}\right\} $, \[ \tilde{\alpha}(T_{1},T_{2})=\sup\left\{ \left|P\left(A\cap B\right)-P(A)P(B)\right|:A\in\sigma\left\langle R(s)|s\in T_{1}\right\rangle ,B\in\sigma\left\langle R(s)|s\in T_{2}\right\rangle \right\} , \] and $\sigma\left\langle R(s)|s\in T\right\rangle $ is the $\sigma$-field generated by the variables $\left\{ R(s)|s\in T\right\} $, $T\subset\mathbb{R}^{2}$.
Suppose there exists a non-increasing function $\alpha_{1}(\cdot)$ with $\lim_{a\rightarrow\infty}\alpha_{1}(a)=0$ and a non-decreasing function $g(\cdot)$ such that \[ \alpha(a,b)\le\alpha_{1}(a)g(b),\ a>0,b>0. \]
Consider $N_{u,n}(h_{k})=\left\{ (i,j)\in N_{n}(h_{k}):j\le j'\ \mathrm{for\ all}\ (i,j')\in N_{n}(h_{k})\right\} $, a set of ordered pairs in $N_{n}(h_{k})$ with unique starting locations $i$, and define the difference sets $N_{r,n}(h_{k})=N_{n}(h_{k})\backslash N_{u,n}(h_{k})$ and $N'(h_{k})=\{1,...,n\}\backslash\left\{ i:(i,j)\in N_{u,n}(h_{k})\ \mathrm{for\ some\ }j\right\} $. We give theory for two-dimensional random field with unknown constant mean. The general case is stated and proved in the Appendix.
The assumptions on the mixing rates are standard (see, for instance, lahiri2002asymptotic and lahiri2003central). The assumption on the bin accuracy $E\left[\left|(\varepsilon(s_{i})-\varepsilon(s_{j}))^{2}-(\varepsilon(s_{i})-\varepsilon(s_{i}+h_{k}))^{2}\right|\right]=o(n^{-1/2})$, while it is realistic for applications such as ours where sampled locations are meant to be spread out somewhat evenly, is strong. It is, however, intrinsic to the matter at hand; the bias in the moment condition must vanish in the $\sqrt{n}$-asymptotics. This theoretical condition connects with our practical experience in simulations; a careful choice of lags for which even small bins $N_{n}(h)$ will contain a large enough number of pairs makes for a noticeably more accurate asymptotic approximation of the distribution. Further note that, as discussed in Section 3, the coverage is close to nominal in simulation using real data.
With the limit distribution of the statistic in hand, the central limit distribution of the minimum-distance estimator obtains under additional identifying assumptions.
The main distributional result may now be stated. Let $g_{j}(\theta)$ be the gradient of $g_{n}(\theta)$ with respect to $j^{\mathrm{th}}$ coordinate of $\theta$.
The density of the observation locations $f$ has an intuitive impact on the asymptotic covariance. As one would expect, if the observations are well spread geographically, this makes for a lower variance because $Q=\int_{R_{n}}f^{2}(x)dx$ is smaller. Correspondingly, cluttered data arranged as a few clusters provides worse information, and the variance is greater for it. Estimation of the asymptotic variance is discussed in the Appendix.
The asymptotic variance of the minimum-distance estimator depends on the covariance matrix of regression errors, $\Sigma$. That is, point estimation is pivotal with respect to $\Sigma$, but large sample inference is not. The reason for this is fairly intuitive. The stability of the point estimate is informed by variability of $R$, which is itself informed by the variability of $Y$, which can only be recovered if it is disentangled from the variability of the error term.
The pivotal point estimation motivates the use of a Monte Carlo sampler using the minimum-distance objective as a basis for its acceptance criteria. Specifically, we want a Monte Carlo procedure that samples a proposed coefficient $\phi^{*}=(\beta^{*},\theta^{*})$ when moments estimated or simulated at $\phi^{*}$ are close to their observed, empirical counterpart. Two approaches naturally suggest themselves: an approximate Bayesian computation (ABC) approach forneron2018abc and a Metropolis-Hastings approach chernozhukov2003mcmc.
Typically, ABC is used as a likelihood-free approach to inference when the likelihood, even up to a constant of proportionality (i.e., even without computing the normalizing constant), is computationally intractable. In this case, however, ABC comes in handy not because it circumvents likelihood computations but because it allows us to provide confidence intervals for the minimum-distance estimator that do not depend on the covariance of regression errors, $\Sigma$.
We develop the ABC sampler for moments built with the sample covariances, but the same approach can be employed with variograms at the cost of estimating a few nuisance parameters.
Consider the following derivation of a large sample approximation to the ABC procedure. For purposes of exposition, first consider in isolation the cross-covariance term at lag $h_{k}$ for a given $k$. The ABC procedure is trying to sample coefficients $\theta$, $m$ and $\beta$ that make differences such as \[ \frac{1}{\left|N_{n}(h_{k})\right|}\sum_{(i,j)\in N_{n}(h_{k})}\left(R(s_{i})-\hat{m}\right)Y(s_{j})-\frac{1}{\left|N_{\tilde{n}}(h_{k})\right|}\sum_{(i,j)\in\tilde{N}_{\tilde{n}}(h_{k})}\left(\tilde{R}(s_{i})-\hat{\tilde{m}}\right)\tilde{Y}(s_{j}) \] small, where $\tilde{R}$ and $\tilde{Y}$ are random fields generated from a generative model evaluated at some candidate value $(\beta^{*},\theta^{*})$, $\tilde{n}$ is the number of synthetic draws, and $\hat{m}$ and $\hat{\tilde{m}}$ are the estimated means of rainfall in the true and synthetic data sets, respectively. The goal of the ABC procedure is thus to select coefficients such that synthetic data generated according to these coefficients makes the above difference between sample moments evaluated on true and synthetic data as small as possible.
A priori, this approach may seem counterintuitive since generating synthetic draws of $\tilde{Y}$ requires stipulating $\Sigma$. However, observing that the moment estimate may be expressed as \[ \frac{1}{\left|N_{n}(h_{k})\right|}\sum_{(i,j)\in N_{n}(h_{k})}\left(R(s_{i})-\hat{m}\right)R(s_{j})\beta+\frac{1}{\left|N_{n}(h_{k})\right|}\sum_{(i,j)\in N_{n}(h_{k})}\left(R(s_{i})-\hat{m}\right)\varepsilon_{j} \] and that $\frac{1}{\left|N_{n}(h_{k})\right|}\sum_{(i,j)\in N_{n}(h_{k})}\left(R(s_{i})-\hat{m}\right)\varepsilon_{j}\rightarrow0$, $\frac{1}{\left|N_{n}(h_{k})\right|}\sum_{(i,j)\in N_{n}(h_{k})}\left(R(s_{i})-\hat{m}\right)R(s_{j})\beta\rightarrow$\\ $\beta Cov_{\theta}(R(s),R(s+h_{k}))$ as $n\rightarrow\infty$, and $\frac{1}{\left|\tilde{N}_{\tilde{n}}(h_{k})\right|}\sum_{(i,j)\in\tilde{N}_{\tilde{n}}(h_{k})}\left(\tilde{R}(s_{i})-\hat{m}\right)\tilde{R}(s_{j})\beta^{*}\rightarrow\beta Cov_{\theta^{*}}(R(s),R(s+h_{k}))$ as $\tilde{n}\rightarrow\infty$, we may instead consider the difference \[ \frac{1}{\left|N_{n}(h_{k})\right|}\sum_{(i,j)\in N_{n}(h_{k})}\left(R(s_{i})-\hat{m}\right)Y(s_{j})-\beta^{*}Cov_{\theta^{*}}(R(0),R(h_{k})) \] as a fit criterion. The full statistic is then a vector of such differences, which we may again designate by $g_{n}$ as it is the same vector statistic used for point estimation, when moments are built with covariances as opposed to semivariances.
The full fit criteria is the quadratic form $l(\phi)=g_{n}(\phi)^{T}\tilde{B}_{n}g_{n}(\phi)$, where $\phi=(\beta,\theta)$. As with point estimation, a positive-definite weighting matrix $\tilde{B}_{n}$ must be used. A naive estimate of the asymptotic variance -- using a covariance matrix proportional to the identity for the regression errors -- performs well. More in the spirit of inference by simulation using a generative model, the inverse of an empirical estimate of the variance of $g_{n}(\phi)$, evaluated on synthetic data from the generative model evaluated at the Krig-and-regress estimates, has been found to work very well and to be easy to implement -- and this estimate can likewise be used in the evaluation of the covariance matrix of the proposal distribution $q$ defined below.
The pseudocode for the procedure is as follows. Draw $\phi^{(0)}$ from $q$.\footnote{There is no burn-in period, but the sample starts at the first accepted coefficient, $\phi^{(1)}$.} For $j=1,...,J$,
where $q$ is a proposal distribution, which we pick to be a naive approximation\footnote{We propose $\phi^{*}$ according to its two-step bootstrap distribution.} to the asymptotic distribution of $\hat{\phi}$, and $\xi$ is a tolerance parameter. In the companion R package SpReg, the default setting is $\xi=0.1$, but it is best to tune the threshold $\xi$ to ameliorate mixing and by inspection of the Monte Carlo distribution $\left\{ \phi^{(1)},...\phi^{(J)}\right\} $, the acceptance ratio, and the trace plots.\footnote{We do not develop on MCMC diagnostics here, but refer the interested reader to hoff2009first, gelman2013bayesian, and pouliot2020lecture.}
As detailed in Table (ref) of Section (ref) below, the heuristic appears to be reliable, slightly over-rejecting in semi-synthetic simulations.
The suggested “asymptotic ABC” is conceptually close to a Metropolis-Hastings algorithm using the minimum-distance objective in lieu of a target/posterior distribution. Alternatively, one could exponentiate the minimum-distance objective, as suggested in chernozhukov2003mcmc. Perhaps surprisingly, this approach seems to perform worse; we obtained substantially worse mixing and coverage for a comparable tuning effort. It has, however, the advantage of not including a tolerance parameter such as $\xi$ for acceptance.
We assess the performance of our suggested methods as well as that of competing approaches on real and hybrid data.
We apply and compare the methods under study using the cross-validation exercise of madsen2008regression. We use the same data set\footnote{We would like to thank the authors for kindly providing their data for replication.} as in their article. As explained therein and further detailed in herlihy1998relationship, the data is a subset of the sample obtained for the Environmental Monitoring and Assessment Program of the Environmental Protection Agency. All samples are from rivers and streams in the American Mid-Atlantic region, and the analysis objective was to relate stream characteristics with land use variables. There are 558 observations over an area of 400,000 squared kilometers. The outcome variables $Y(x)$ are the logarithm of chloride concentration at locations $x$, and the independent variables $R(x^{*})$ are the logit transformations of the percent of watershed in forest at locations $x^{*}$.
The reference value is obtained by doing generalized least-squares on the full, aligned data set weighting with an estimate of the regression error covariance: we obtain $\hat{\beta}_{\mathrm{full}}=-0.38$. The simulation is implemented as follows; for each run, a randomly chosen half of the independent variables, the $R_{i}$'s, are “hidden”, and the outcome variables, the $Y_{i}$'s, are “hidden” for the other half of the data set, thus creating a misaligned data set. For each round of the cross-validation exercise, $\beta$ is estimated with each method, and estimates are recorded.
We see from the empirical variogram displayed in Figure (ref) that neighboring dependent variables do covary, thus allowing for useful interpolation as a first step.
Table (ref) presents hybrid data on a regular lattice. This special case is important as this is a common design. Rainfall is generated according to the exponential covariance model with parameters estimated on the aligned rivers data, but data is generated at locations on a lattice instead of the original, irregular locations. The outcome variable is generated according to the linear regression model with parameters estimated on the aligned rivers data. The mean $E[\hat{\beta}]$ is the average point estimate over all simulation draws. The simulation standard errors $\sqrt{V(\hat{\beta})}$ are the standard deviations of the simulation draws, and the output standard errors $\sqrt{\hat{V}(\hat{\beta})}$ give the average standard errors over all simulation draws. The squared-root mean-squared errors (RMSE) are defined analogously. For the two-step bootstrap, the mean is the average mean over the bootstrap samples, and the regression output standard errors are given by the average over simulation runs of the standard deviation of bootstrap estimates. The coverage corresponds to the fraction of times the confidence interval, computed with the output standard errors of the current run of the simulation, covered -0.38.
On the lattice, as detailed in Table (ref), the efficiency gains of the minimum-distance approach deliver the expected improvement in squared-root mean-squared error (RMSE) of the regression coefficient estimate. The Krig-and-regress estimator, with naive inference ignoring first-stage uncertainty, severely undercovers while the two-step bootstrap estimator has correct coverage.
Table (ref) presents the output of the cross-validation exercise using the observed data. While both the Krig-and-regress and the minimum-distance estimator do well in terms of RMSE and coverage accuracy, the Krig-and-regress estimator outperforms the Minimum-Distance estimator in both. This is explained in part to the difficulty of estimating the weighting matrix in such a small sample, which is important for to the quality of point estimation with the minimum-distance estimator.\footnote{One avenue we have found fruitful for weighing matrix estimation in small samples is to simulate data from the posited model evaluated at the Krig-and-regress estimates, and compute the variance of synthetic sample moments using the observed location and simulated observations. We leave further exploration of this approach for follow-up research.}
As noted in Table (ref), estimation of the asymptotic variance of the minimum-distance estimator requires stipulation of $\Sigma$. The pivotal ABC alternative appears to do relatively well, even though it over-rejects by about 5% in our simulations, where the interval has 95% nominal coverage.
In Under the Weather: Health, Schooling, and Economic Consequences of Early-Life Rainfall, maccini2009under estimate the effect of a rainfall shock in infancy on adult socioeconomic outcomes such as education or health. The paper merges a rainfall measurements data set with survey data, both of which have geographically located observations. The data sets are misaligned, as can be seen from Figure (ref), which plots the respective locations.
The data set for the regression is obtained by merging two misaligned data sets. The first data set contains rainfall measurements from measuring stations across Indonesia. The whole rainfall data spans the years from 1864 to 2004.
Only the years 1953-1975 are used to identify yearly variation; the other years are used to estimate long term averages. In almost every year used in the analysis, more than 300 rainfall measurement stations are active. The rainfall data comes from the Global Historical Climatology Network (GHCN), which is a publicly available data set.\footnote{Available online at http://www.ncdc.noaa.gov/oa/climate/research/ghcn/ghcn.html.}
The second data set is the third wave of the Indonesian Family Life Survey (IFLS3), which includes each surveyed individual's year and location of birth. It consists of 4,615 women and 4,277 men born outside large cities between 1953 and 1974.
The locations of the birthplaces and rainfall stations are given in Figure (ref). We can see that the data sets are misaligned. We can also see from Figure (ref) that most birthplaces are situated fairly close to one or more rainfall stations; the median distance to the closest rainfall station is 11.43 km, and the third quartile distance is 90.28 km. The median distance to the fifth closest station is 30.24 km, and the third quartile distance to the fifth closest station is 317.10 km. We show below that, at those distances, the rainfall measurements are still correlated; thus, informative interpolation is possible.
We use, as maccini2009under did, the log of the ratio of the yearly rainfall with the long run average yearly rainfall. We find in Figure (ref) that this has the benefit of eliminating the point mass at zero, and making the distribution “closer to Gaussian.”
The typical variogram fit of Figure (ref) suggests a good fit of the variogram and corroborates the assumption of a Gaussian covariance function (ref). However, stein2012interpolation warns against relying on such plots to draw definitive conclusions, and likelihood fits were likewise investigated. Altogether, this suggests that Kriging, which is the best predictor under Gaussianity, ought to produce reliable interpolation.
The first order of business is to implement the two-step estimator. As mentioned above, we carry out Kriging on the log ratio of yearly rainfall to long run yearly average. The transformed data has the additional virtue of making immediate the comparison with the regression output of maccini2009under. Remark that if the imputed covariate were a nonlinear transformation of the interpolated quantity, the estimated coefficients may not be consistent estimators of the best linear predictor coefficients.
A potential concern is the isotropy assumption, i.e., the direction of the vector giving the difference between the locations of two points in the random field does not matter. That is, for the variogram to only be a function of the distance between its arguments and not the direction of the vector from one to the other, we must assume that the direction does not impact the covariance. One way to assess this is to plot and inspect the directional variograms, a typical (for this data) example of which is shown in Figure (ref).
The directional variograms are reassuring. The empirical directional variograms seem to align well with the fitted isotropic variogram up to at least 500 km. Beyond that distance, very little data goes into estimating each point in the plots; hence they are quite variable. The appearance of a trend away from the fitted isotropic variogram can be due to the fact that these points are highly correlated.
We are interested in fitting the model
where $\mathcal{B}$ and $\mathcal{G}$ are the set of observation indices corresponding to subjects who are boys and girls, respectively. The random variables $R_{\mathrm{true},i}$ are sampled from the random field of the log of the ratio of yearly rainfall to long term average at the outcome locations, and $F$ includes location (district) dummies, season dummies, time trend, and interactions. $R_{\mathrm{true},i}$'s are not observed but $\mathbf{R}^{*}$, a vector of observations from the same random field albeit at different locations, is observed.
Remark that we use a different two-step method than maccini2009under did. Whilst we do the interpolation step with the best linear predictor, they use an instrumental variables approach. Their strategy is to estimate (ref) by running a two-stage least-squares regression in which the nearest rainfall is instrumented with slightly more distant rainfall stations (the second to fifth closest stations). We argue that this approach is problematic for two reasons. Conceptually, there is no reason to believe that instrumenting the wrong rainfall measurements with other wrong rainfall measurements will “wash away” the error term, which can be thought of as the difference between the rainfall at the outcome location and the rainfall at the nearest station (for instance, all stations could be near each other and have correlated “errors”). However, such a strong assumption is necessary for two-stage least-squares approach to deliver a consistent estimator. Practically, the data set at hand makes the first-stage regression coefficient difficult to interpret and makes this strategy tricky to implement. Indeed, the closer rainfall stations do not have the same distance or configuration around any outcome location. Furthermore, for any given outcome location, there will be few years for which the nearest rainfall station and all of the next four nearest stations will have observed measurements, so the, say, second nearest station does not correspond to the same station for every observation.
To further assess the quality of the Kriging estimates, we plot the geographically located interpolation errors from leave-one-out cross-validation in Figure (ref). The figure suggests that the residuals have low spatial correlation.
The second stage for the Krig-and-regress is done by ordinary least-squares. We find that using feasible generalized least-squares did not help in this case.
The difference between the IV and the Krig-and-regress methods is noticeable. This should not be surprising as they correspond to completely different approaches to imputing the missing covariate in the first step. The economic magnitude of the difference is relevant for policy implications; the Krig-and-regress point estimate is 0.115 higher than that of the IV approach. This means that the estimated impact of the first standard deviation in yearly rainfall, away from the long term average, on the estimated probability of declaring oneself very healthy increases from 0.4% to 3.5%.
The valid standard errors of the two-step bootstrap tend to be larger than the naive, unprincipled ones but moderately so and do not seem to pose a threat to statistical significance in this application; note that there is nothing inherently conservative about their design.
The only important disparity in the minimum-distance regression output is for completed grades of schooling. It suggests a significant effect where Krig-and-regress does not, but for a smaller point estimate than the IV estimator does in maccini2009under. We find that, for girls, a one standard deviation increase in rainfall brings about an increase of 0.15 in years of schooling, as opposed to 0.42.
As remarked in the simulation of Section (ref), the standard errors of the minimum-distance estimator can be tangibly different from those of the Krig-and-regress estimator, sometimes smaller and sometimes larger, but they were observed to be generally accurate.
We developed methods handling regression analysis with misaligned data and paid particular attention to the case in which the researcher does not want to specify the covariance structure of the regression errors, thus making maximum likelihood estimation for the full model inaccessible. First, we suggested a modified two-step method which produces confidence intervals that account for the uncertainty due to the estimation of the imputed regressor and is very easy to compute. This contribution is important because, although Krig-and-regress was already recommended for two-step estimation in the misaligned data regression problem, to the best of our knowledge, the literature did not offer standard errors that account for estimation in the first stage without requiring specification of the regression error covariance structure. Second, we proposed a one-step minimum-distance estimator and developed its limit distribution theory. We produced an asymptotic approximation formula for the covariance of the estimator as well as a likelihood-free sampling strategy that provides pivotal inference with respect to the covariance matrix of the regression errors.
First and foremost, I would like to thank my PhD advisors Gary Chamberlain, Edward Glaeser, Neil Shephard, and Elie Tamer for their guidance and support. I am indebted for their insightful comments to Alberto Abadie, Nathaniel Hendren, James Heckman, Guido Imbens, Bruce Meyer, Jann Spiess, Michael Stein, and Bryce Millett Steinberg. Zhen Xie has provided outstanding research assistance.