Biometrics 70, 648–660 September 2014

DOI: 10.1111/biom.12175

A Pseudo-Penalized Quasi-Likelihood Approach to the Spatial Misalignment Problem with Non-Normal Data Kenneth K. Lopiano,1, * Linda J. Young,2 and Carol A. Gotway3 1

Statistical and Applied Mathematical Sciences Institute, RTP, North Carolina, U.S.A. 2 Department of Statistics, University of Florida, Gainesville, Florida, U.S.A. 3 Centers for Disease Control and Prevention, Atlanta, Georgia, U.S.A. ∗ email: [email protected]

Summary. Spatially referenced datasets arising from multiple sources are routinely combined to assess relationships among various outcomes and covariates. The geographical units associated with the data, such as the geographical coordinates or areal-level administrative units, are often spatially misaligned, that is, observed at different locations or aggregated over different geographical units. As a result, the covariate is often predicted at the locations where the response is observed. The method used to align disparate datasets must be accounted for when subsequently modeling the aligned data. Here we consider the case where kriging is used to align datasets in point-to-point and point-to-areal misalignment problems when the response variable is non-normally distributed. If the relationship is modeled using generalized linear models, the additional uncertainty induced from using the kriging mean as a covariate introduces a Berkson error structure. In this article, we develop a pseudo-penalized quasi-likelihood algorithm to account for the additional uncertainty when estimating regression parameters and associated measures of uncertainty. The method is applied to a point-to-point example assessing the relationship between low-birth weights and PM2.5 levels after the onset of the largest wildfire in Florida history, the Bugaboo scrub fire. A pointto-areal misalignment problem is presented where the relationship between asthma events in Florida’s counties and PM2.5 levels after the onset of the fire is assessed. Finally, the method is evaluated using a simulation study. Our results indicate the method performs well in terms of coverage for 95% confidence intervals and naive methods that ignore the additional uncertainty tend to underestimate the variability associated with parameter estimates. The underestimation is most profound in Poisson regression models. Key words:

Berkson Error; Generalized linear models; Kriging; Spatial Misalignment.

1. Introduction Data associated with geographical units are routinely used to answer questions in both basic and applied science. Environmental epidemiologists assess the relationships between environmental exposure and human health; soil scientists consider the relationship between weather and soil chemistry; atmospheric scientists seek to understand the relationship between ground based measurements and data collected from satellites and other remote sensing technology. Regardless of the application, the variables are often observed at different locations or aggregated over different geographical units. Such datasets are said to be spatially misaligned. Although we use the terminology point-to-areal misalignment, it should be noted that points and areas are never misaligned, rather there is a change-of-support problem at hand (Gotway and Young, 2002). Two common types of misalignment problems are pointto-point and point-to-areal misalignment. In both cases, a covariate is observed at a set of points in geographic space. For point-to-point misalignment, the response is observed at a different set of points, while for point-to-areal misalignment, the response is associated with an areal unit, for example, a county. The disparate datasets must be aligned on a common spatial support before inference about the relationship between the two variables can be conducted. Additional uncer-

648

tainty is introduced when kriging or other smoothing methods are used to align the spatially referenced data (Madsen, Ruppert, and Altman, 2008; Buonaccorsi, 2009; Gryparis et al., 2009; Paciorek et al., 2009; Lopiano, Young, and Gotway, 2011, 2013; Szpiro, Sheppard, and Lumley, 2011). Although the misalignment problem has been considered for models with normally distributed responses (Madsen et al., 2008; Gryparis et al., 2009; Lopiano et al., 2011, 2013; Szpiro et al., 2011), accounting for the additional uncertainty induced from using kriging to align disparate spatial datasets with non-normal responses has received less attention. A fully model-based Bayesian approach was presented in Mugglin, Carlin, and Gelfand (2000), Banerjee and Gelfand (2002), and Zhu, Carlin, and Gelfand (2003). Gryparis et al. (2009) briefly discussed measurement error in probit regression models and conducted a simulation study for a logistic regression model. In the simulation study, the authors compared the naive estimator found by using the kriging mean as if it were the true value with an out-of-sample regression calibration (RC-OOS) estimator (Thurston, Spiegelman, and Ruppert, 2003) and an exposure simulation approach. The exposure simulation approach is intuitively appealing because it attempts to account for both the uncertainty due to using the kriging mean and that due to estimating the parameters that define the kriging mean and kriging variance. However, the approach actually © 2014, The International Biometric Society

PPQL for Spatial-Misalignment with Non-Normal Data induces bias and results in an incorrect imputation scheme. Thus, the authors recommend against the exposure simulation approach. In this work, we consider the problem where the response and a covariate are spatially misaligned, kriging is used to predict the covariate at the locations where the response is observed, and the response is assumed to be non-normal. For normally and Poisson-distributed responses using the canonical link, naive analyses that ignore the fact that the kriging predictions are not the true covariates result in unbiased estimates of the parameters. However, the standard errors associated with the parameter estimates will be underestimated if the additional uncertainty is ignored. For a generalized linear model with binomially distributed responses using the canonical link, naive estimates are biased and underestimate the uncertainty (Buonaccorsi, 2009). Our main contribution is a method to account for the additional uncertainty from using kriging predictions to align disparate datasets in a generalized linear modeling framework. More specifically, we show the additional uncertainty can be accounted for using standard techniques developed for generalized linear mixed models. In doing so, the parameter and an associated measure of uncertainty are estimated in a single procedure. For concreteness, we illustrate the proposed methodology with two examples: (1) relating Florida county asthma counts to PM2.5 concentrations after the onset of the Bugaboo scrub fire, the largest wildfire in state history; and (2) relating the probability of low birth weight to PM2.5 concentrations after the onset of the Bugaboo scrub fire. The applications here are representative of a larger class of applications related to understanding the health effects of extreme weather events such as wildfires. The performance of the proposed methodology is assessed using a simulation study. The structure of the article is as follows. In Section 2 we define the spatial misalignment modeling framework, show equivalence between point-to-point and point-to-areal misalignment and develop a pseudo-penalized quasi-likelihood estimator to estimate the regression parameters in generalized linear models with spatially misaligned data. In Section 2.2 we describe the difference between estimators that ignore the additional uncertainty and the proposed pseudo-penalized quasilikelihood estimator. In Section 3 we provide a point-to-point application and a point-to-areal application. In Section 4 we present the design of a simulation study used to assess the proposed estimator and a subset of the results. Finally, we end the article with a discussion in Section 5. 2. Proposed Methodology For the case of point-to-point misalignment, suppose the response variable is observed at n spatial locations su =



X(su )n×1 X(so )no ×1



 ∼N

C(su )n×k ξ C(so )no ×k ξ

[su,1 , . . . , su,n ], and the explanatory variable is observed at no different spatial locations so = [so,1 , . . . , so,no ]. The gener-

649

alized linear model relating the response observed at su , denoted as Y(su ) = [Y (su,1 ), . . . , Y (su,n )] , and the unobserved values of the covariate at the locations su , denoted as X(su ) = [X(su,1 ), . . . , X(su,n )] , is Y(su )n×1 ∼ Exponential Family(mean = μn×1 , scale = φ) g(μn×1 ) = β0 1n×1 + β1 X(su )n×1 , where g(·) is the link function and the subscripts indicate the dimensions of the vectors. The point-to-point spatial misalignment problem arises because Y(su ) is observed, X(su ) is unobserved, and X(so ) = [X(so,1 ), . . . , X(so,no )] is observed. For the case of point-to-areal misalignment, suppose the response variable is observed at nB areal units and the covariate is observed at no spatial locations so = [so,1 , . . . , so,no ]. Let B denote the set of n areal units. For each areal unit b = 1, . . . , nB in B, X(b) = |b|−1 s∈b X(s) ds is the average value of the covariate in unit b, where |b| is the area of unit b. The generalized linear model relating the observed values of the response at the areal units B, denoted as Y(B) = [Y (1), . . . , Y (nB )] , and the unobserved values X(B) = [X(1), . . . , X(nB )] is Y(B)nB ×1 ∼ Exponential Family(mean = μnB ×1 , scale = φ) g(μnB ×1 ) = β0 1nB ×1 + β1 X(B)nB ×1 , where the subscripts indicate the dimensions of the vectors. The point-to-areal spatial misalignment problem arises because Y(B) is observed, X(B) is unobserved, and X(so ) is observed. Although point-to-point and point-to-areal misalignment are different, the modeling structure is identical if kriging and block kriging are used, respectively, to align the disparate datasets. To see this, suppose the response and covariate are observed in a two dimensional geographic space D ⊂R2 . Let X(s) ∼ Gaussian Process(c(s)k×1 ξ k×1 , KX (·)),

(1)

where c(s) is a k × 1 vector of covariates observed at each location s ∈ D, ξ k×1 is an unknown vector of k regression coefficients, and KX (s, s ; τ x ) is a covariance function defined for any two locations s and s parameterized by an unknown vector of l parameters τ x . For example, τ x = [scale, range] may be the l = 2 parameters that define an exponential covariance structure. For the case of point-to-point misalignment, the joint distribution of unobserved values of the covariate, X(su ), and the observed values of the covariate, X(so ), is



 , X (τ x )n+no ×n+no =

X,uu X,uo X,ou X,oo

 ,

where C(su ) and C(so ) are n × k and no × k matrices of observed covariates at the locations in su and so , respectively, and the covariance matrix X (τ x ) has KX (si , sj ; τ x ) as the

650

Biometrics, September 2014

(i, j)-th element. The covariance matrix X is partitioned into the covariance matrices X,uu , the n × n covariance matrix among the unobserved locations, X,uo , the n × no covariance matrix among the unobserved and observed locations, and X,oo , the no × no covariance matrix among the observed locations. The conditional (kriging) distribution of X(su )n×1 |X (so )no ×1 is X(su )n×1 |X(so )no ×1 , ξ, τ x ∼ N(Wn×1 = C(su )ξ + X,uo −1 X,oo × (X(so ) − C(so )ξ)xn×n = X,uu − X,uo −1 X,oo X,ou ). Notice the kriging mean has a Berkson error relationship with the unobserved covariate, X(su )n×1 |X(so )no ×1 , ξ, τ x = Wn×1 + xn×1 , where xn×1 ∼ N(0n×1 , xn×n ) and the Berkson parameters are λ = {ξ, τ x }. Because of this Berkson error relationship, the generalized linear model relating Y(su ) and X(su ) becomes Y(su )n×1 ∼ Exponential Family(mean = μn×1 , scale = φ) g(μn×1 ) = β0 1n×1 + β1 Wn×1 + β1 xn×1 xn×1 ∼ N (0n×1 , xn×n ) .

(2)

Similarly, for the case of point-to-areal misalignment, block kriging can be used to align the disparate datasets. Given the observed values X(su ), the block kriging distribution is ¯ X(B)nB ×1 |X(so )no ×1 ∼ N(WnB ×1 = C(B)ξ + B,so −1 so (X(so ) − C(so )ξ), xnB ×nB = B,B − B,so −1 so so ,B ) ¯ where C(B) is a nB × k matrix the i, jth entry is the average of the jth covariate within block b, b = 1, . . . , nB , B,B is the nB × nB covariance matrix between the areal units, B,so is the nB × no point-to-areal covariance and, so is the no × no covariance matrix among the points. As in the case of kriging, the block kriging mean has a Berkson error relationship with the unobserved covariate, X(B)nB ×1 |X(so )no ×1 , ξ, τ x = WnB ×1 + xnB ×1 ,

Y(B)nB ×1 ∼ Exponential Family(mean = μnB ×1 , scale = φ) g(μnB ×1 ) = β0 1nB ×1 + β1 WnB ×1 + β1 xnB ×1



xnB ×1 ∼ N 0nB ×1 , xnB ×nB .

Yn×1 ∼ Exponential Family(mean = μn×1 , scale = φ) g(μn×1 ) = β0 1n×1 + β1 Wn×1 + β1 xn×1 xn×1 ∼ N (0n×1 , xn×n ) ,

(4)

where W represents the kriging mean and x represents the kriging variance. The parameters that define W and x are defined as the Berkson parameters λ = [ξ, τ x ]. The equivalent modeling framework allows the simultaneous development of methods to estimate the regression parameter β1 for both types of misalignment. For simplicity, with the notation established, the subscripts indicating the dimensions of the vectors and matrices will be dropped. Under the model described in (4), our goal is to estimate the parameter β1 and an appropriate measure of uncertainty that accounts for the additional error term x . Note that the relationship in (4) is similar to the specification of a generalized linear mixed model (GLMM) where x plays the role of a random effect. One key difference, however, is the multiplication of β1 by the random effect. In light of this difference, we propose adapting the methods used to estimate parameters in generalized linear mixed models to estimate β. We adapt the penalized quasi-likelihood (PQL) method developed by Breslow and Clayton (1993). Breslow and Clayton (1993) proposed a Fisher-scoring algorithm to estimate regression parameters and to predict random effects in a GLMM. For use with the Berkson error model in (4), we adapt the algorithm in two ways: (1) We account for the presence of β1 in both the fixed and random component, and (2) We use known or externally estimated Berkson parameters; that is, we assume λ and subsequently W and x are known. The adapted algorithm is: β , using the estimated (1) Initialize the estimate of β, value of β obtained from ignoring the additional error (i) term. Initialize x , x , as 0n×1 . (2) Define the working response (i)

(i) (i) (i) (i) Z(i+1) = β0 1n×1 + β1 W + β1 x

where xnB ×1 ∼ N(0nB ×1 , xnB ×nB ). Similarly, the generalized linear model relating Y(B) and X(B) becomes



equivalence, we simplify the notation and consider the following model to represent both the point-to-point misalignment problem and the point-to-areal misalignment problem where kriging and block kriging are used, respectively, to align the disparate datasets;

(3)

The models in (2) and (3) are equivalent in structure with only the definition of W and x changing. Because of the





(i) (i) (i) (i) + g g−1 β0 1n×1 + β1 W + β1 x







(i) (i) (i) (i) × Y − g−1 β0 1n×1 + β1 W + β1 x



where g is the derivative of the link function. Unlike in typical GLMMs where a known matrix is multiplied by the random effect, β1 is multiplied by the random effect. As a result, here we plug-in the working estimate (i) of β1 , β1 .

PPQL for Spatial-Misalignment with Non-Normal Data

651

(3) Estimate D, the usual diagonal GLM weight matrix, (i+1) , where as D

jj D

(i+1)

=

wj





2 (i) (i) (i) (i) (i) (i) (i) (i) V g−1 β0 + β1 Wj + β1 j,X|W g g−1 β0 + β1 Wj + β1 j,X|W

where Wj is the jth element of the vector W, wj are potential weights and V (·) is the usual variance function (see McCullagh and Nelder, 1989).



(i+1) )−1 + β 1 (i+1) = (D (4) Calculate 

(i)

2

x . In typical

GLMMs, the parameters that define x are estimated using the approximate profile quasi-likelihood function. Here, because the parameters that define x are assumed known or estimated externally, the additional step required to estimate the parameters is omitted. (5) Calculate x W = [1n×1 W] (6) Calculate

(i+1)

(i+1) )−1 (Z − W = β1 x ( β ), where (i)

(i)

(i+1) (i+1) )−1 W)−1 W ( (i+1) )−1 Z. β = (W (

When Berkson error arises from spatial misalignment, the observed X(so ) can be used to estimate ξ and τ x with EGLS and restricted maximum likelihood. That is, τˆx is found by maximizing the restricted likelihood of X(so ) and  ˆ −1 −1  ˆ −1 ˆξ EGLS = (C(so ) X,oo C(so )) C(so ) X,oo X(so ),

ˆ X,oo = X,oo (ˆ where  τ x ). Then, the algorithm above can be ˆ = {ˆξ EGLS , τˆx }. used with λ Here, as in Lopiano et al. (2013), we account for the error associated with estimating ξ. For illustration, we use the point-to-point misalignment framework. A similar argument can be shown for point-to-areal misalignment. The approximate conditional distribution of X(su )|X(so ), τ x is

(7) Repeat Steps 2–6 until convergence. The covariance of the final estimator is estimated via

−1 W)−1 . Var( β) = (W 

. ˜ = C(su )˜ξ GLS + X,uo −1 (X(so ) X(su )|X(so ), τ x ∼ N(W X,oo

− C(so )˜ξ GLS ), ˜ x = X,uu − X,uo −1  X,oo X,ou  + (C(su )−X,uo −1 X,ooC(so ))(C(so )

2.1. Estimating the Berkson (or Kriging) Parameters λ The algorithm above assumes λ is known, which is rarely, if ever, the case. As a result, the algorithm only accounts for the additional uncertainty induced by using the kriging mean as the covariate. It does not account for the variability associated with estimating the kriging mean or kriging variance which depend on the parameters λ, the Berkson parameters. However, if a dataset that depends on λ is available, the dataset can be used to obtain an estimate of λ. Then the ˆ can be used as if they are the true values in the estimates λ algorithm detailed above. Historically, estimators found using likelihood methods with estimated parameters treated as known are referred to as pseudo-likelihood estimators (Gong and Samaniego, 1981). Our approach builds directly on the penalized quasi-likelihood method for estimating parameters in generalized linear mixed models introduced by Breslow and Clayton (1993). Thus, we refer to the two-step estimator defined here as a pseudo-penalized quasi-likelihood estimator (PPQL estimator). Although a well recognized limitation of pseudo-likelihood approaches is the fact that estimating λ introduces additional uncertainty that can cause bias in parameter estimates and standard errors to be too small, if the Berkson parameters can be estimated with sufficient precision, then the additional uncertainty and bias will be negligible (Buonaccorsi, 2009). Nonetheless, we address how to account for the additional uncertainty induced from estimating the Berkson parameters.

× oo C(so ))−1  × (C(su ) − X,uo −1 X,oo C(so )) ).

(5)

The algorithm proposed in the preceding section can be refor˜ and  ˜ x , respectively. mulated by replacing W and x with W In doing so, the estimators account for the additional error associated with estimating ξ but not that associated with estimating τ x . Although it would seem reasonable to remove the dependence on τ x , and hence the variability associated with it, as well, an approximate distribution of X(su )|X(so ) does not exist in closed form and a fully Bayesian approach would be more appropriate. For a discussion of Bayesian methods in spatial misalignment problems see Mugglin et al. (2000), Gelfand, Zhu, and Carlin (2001), Banerjee and Gelfand (2002), and Zhu et al. (2003). However, as shown later through simulation, the effect of the additional uncertainty associated with estimating τ x is negligible. 2.2. Comparing Naive and Adjusted Estimators Naive estimators are the estimators calculated when one ignores the error term x and uses W as if it were the true value. Here we focus our attention to the case when the response is assumed to follow a Poisson distribution or a binomial distribution.

652

Biometrics, September 2014

2.2.1. Poisson and Binomial. Ignoring the Berkson error, the covariance of the estimated regression parameters is  −1 ˆ Var(β naive ) = (W DW) ,

(6)

where D is a diagonal matrix of the GLM weights. In the case of Poisson regression with a log link, the ith diagonal element of D is exp(β0 + β1 Xi ), and in the case of logistic regression, the ith diagonal element of D is nT

exp(β0 + β1 Xi ) , (1 + exp(β0 + β1 Xi ))2

where nT is the number of trials. The PPQL-based variance estimate is  −1 −1 ˆ Var(β PPQL ) = (W  W) ,

(7)

where  = D−1 + β12 x . Comparing (6) and (7) we see the difference in the variance estimates depends on β1 and the parameters that define x . We explore this fact in the simulation study. 3. Application In 2007, during the months of April, May, and June, the largest fire in the history of Florida and Georgia, the Bugaboo Scrub Fire, raged in southern Georgia and north central Florida. Smoke from the fire reached Atlanta, south Florida and even parts of Alabama and Mississippi (Georgia Forestry Commision, 2007). Smoke from the early stages of the fire poured into Florida from Georgia. On April 29, the smoke was relatively concentrated. However, on April 30, smoke from the fire blanketed most of north central Florida and even parts of central Florida. On May 3, then Governor Charlie Crist signed executive order 07-86 prompting a fire assistance grant from the Federal Emergency Management Agency (The State of Florida State Emergency Response Team, 2007). On May 10, a statewide health advisory was released by the Florida Department of Environmental Protection: “Caution Statement: The elderly, children, and those with respiratory diseases including asthma or heart disease should avoid prolonged outdoor activity; everyone else should limit prolonged outdoor activity” (Florida Department of Environmental Protection, 2007). 3.1. Point-to-Point Misalignment The time period consisting of weeks 3–8 of human gestation is known as embryogenesis. During embryogenesis, vital organs and limb structures develop. At the end of embryogenesis, the offspring is officially referred to as a fetus. Pregnancy is considered full term after 37 weeks of gestation (Sadler, 2011). During embryogenesis, the embryo is vulnerable to developmental defects. Because of this, we hoped to assess whether high levels of particulate matter during embryogenesis are related to an increased probability of low birth weight where

a low birth weight is defined as a fetus being born weighing less than 2,500 g. The PM2.5 value defined as the maximum amount of PM2.5 during the month of May was observed at each of 32 Air Quality System (AQS) air quality monitors in the state of Florida (Figure 1, U. S. Environmental Protection Agency (2007), Supplementary File 1). Because the monitors have different sampling frequencies, for example, daily or once every three days, the number of observations from each station varies. The minimum number of samples during the month of May was 8 and the maximum number of samples was 35. Of the 32 monitors, 18 had less than or equal to 15 samples, 3 had between 26 and 30 samples and the remaining 11 had at least 31 samples. For a complete listing of the number samples from each site see Supplementary File 1. The lack of daily measurements from all monitors is a noted limitation of our study. In 2007, infants of gestational age 37 weeks born between December 24 and December 31 experienced the majority of embryogenesis during the month of May. The birth weight data obtained from the Florida Department of Health’s Office of Vital Statistics consisted of 179 non-Hispanic white mothers who gave birth during this time period to one infant 37 weeks in gestational age. Let so denote the coordinates (latitude/longitude) of n = 32 the monitors. Let the vector X(so ) be the vector of the PM2.5 values at the locations so . Suppose X(so ) ∼ N(132×1 ξ, X,oo (τ x )).

(8)

where τ x is a 1 × 2 vector of the parameters that define an exponential covariance structure, that is, the scale and range. In comparing three standard covariance structures, spherical, Gaussian, and exponential, it was determined the exponential covariance structure performed best in terms of the Akaikie information criterion (AIC) and Bayesian information criterion (BIC). The parameters ξ and τ x were estimated based on restricted maximum likelihood and estimated generalized least squares as ˆξ = 53.28 and τˆx = [856.04, 0.36]. Let su denote the coordinates of the nu = 179 mother residential locations. The conditional (kriging) distribution of X(su )|X(so ) is X(su )|X(so ), ξ, τ x ∼ N(W = 1179×1 ξ + X,uo −1 X,oo × (X(so ) − 132×1 ξ), x = X,uu − X,uo −1 X,oo X,ou ). (9) Let Y (su,i ) be an indicator of low-birth weight for the ith infant, i = 1, . . . , 179. Suppose the relationship between Y(su ) = [Y (su,1 ), Y (su,2 ), . . . , Y (su,179 )] and X(su ) = [X(su,1 ), X(su,2 ), . . . , X(su,179 )] follows a logistic regression model. Other covariates (mother’s age, mother’s education, mother’s smoking status, mother’s alcohol use, mother’s Women, Infant, Children (WIC) status) were considered, but were not significant at a 10% significance level. The model described here has only PM2.5 as a covariate: Y(su ) ∼ Bernoulli(π(su )) logit(π(su )) = β0 1179×1 + β1 X(su )

(10) (11)

PPQL for Spatial-Misalignment with Non-Normal Data

653

28 24

26

Latitude

30

32

PM 2.5 Monitors

−86

−84

−82

−80

Longitude

Figure 1. PM2.5 monitor locations in Florida. Because X(su ) is unobserved, we instead considered the model Y(su ) ∼ Bernoulli(π(su ))

(12)

logit(π(su )) = β0 1179×1 + β1 W + β1 x ,

(13)

x ∼ N(0179×1 , x ).

(14)

Using the estimates of ξ and τ x , the approximate conditional distribution in (5) was calculated. To estimate the parameters ˜ and in the regression model, we used the estimates of W, W, ˜ x , calculated using the estimated values of ξ and τ x in x ,  the pseudo-penalized quasi-likelihood algorithm to account for additional term β1 x and the estimation of ξ. In addition,

˜ as the true we considered the naive model that treats W covariate. For the naive approach, the estimate of β1 was 0.056 with a variance and standard error of 0.000308 and 0.018, respectively. When using penalized quasi-likelihood, the estimate of β1 was 0.058 with a variance and standard error of 0.000370 and 0.019, respectively. Accounting for the additional uncertainty due to kriging resulted in a 20% increase in the variance of the estimate of β1 , 0.000308–0.000370, and a 10% increase in the standard error, 0.018–0.019. The parameter estimate increased 3% from 0.056 to 0.058. A Wald test for the significance of the parameter resulted in a z-score of 3.19 in the naive case and 3.00 in the

654

Biometrics, September 2014

PPQL case, a 6% decrease. However, in both cases, the effect of PM2.5 was significant. Using the PPQL based estimate of β1 , during the time period under consideration, we estimated a unit increase in PM2.5 resulted in an estimated multiplicative effect on the odds of low birth weight of 1.059 with a 95% confidence interval of (1.02, 1.10). 3.2. Point-to-Areal Misalignment In light of the sudden onset of dire air quality conditions at the end of April and the increased public awareness of the wildfires by May 10, we hypothesized asthma hospital and emergency department (ED) admissions would be correlated with PM2.5 during the early stages of the fire. Because behavior likely changed shortly after the onset of the first extreme smoke event and the emergency declaration, we considered the total number of hospitalizations and ED cases between April 30, 2007 and May 6, 2007. The PM2.5 value defined as the maximum value of PM2.5 between April 30, 2007 and May 6, 2007 was observed at each of 32 Air Quality System (AQS) air quality monitors in Florida (U. S. Environmental Protection Agency (2007), Supplementary File 2). Florida’s Department of Health (FDOH) has a data-sharing agreement (DSA) with Florida’s Agency for Health Care Administration (AHCA), which provided access to two data sources: confidential hospitalization records and ED records. The health data consisted of all records where either the principal or other diagnosis codes were for asthma, ICD-9 codes 493.0-493.9 or 786.07 (National Center for Health Statistics, 1998). To account for demographic variation, the expected number of health events was calculated for each county adjusting for age, gender and ethnicity using the state of Florida as the standard population (Web Appendix A). Let Y (c) and E(c) denote the observed and expected number of cases in county c, respectively, c = 1, . . . , 67. Let X(s) denote the PM2.5 value at location s, s ∈ D ⊂R2 , where D denotes the state of Florida. As in the previous example, the number of observations from each station varies. The minimum number of samples during the specified time period was 1 and the maximum number of samples was 9. Of the 32 monitors, 12 had less than or equal to 3 samples, 7 had either 4 or 5 samples and the remaining 13 had at least 7 samples. For a complete listing of the number of samples from each site see Supplementary File 2. The lack of daily measurements from all monitors is a noted limitation of our study. Suppose X(s) ∼ Gaussian Process(ξ, KX (·)),

(15)

where KX (s, s ; τ x ) is an exponential covariance function defined for any two locations s and s parameterized by an unknown vector of 2 parameters τ x = [scale, range]. For each  county c, X(c) = |c|−1 s∈c X(s) ds is the average PM2.5 value in county c, where |c| is the area of county c. Let so denote the coordinates (latitude/longitude) of the n = 32 monitors. Let the vector X(so ) = [X(so,1 ), X(so,2 ), . . . , Y (su,32 )] be the vector of the PM2.5 values at the locations so . Note (15) implies X(so ) ∼ N(132×1 ξ, X,oo (τ x )). As in the previous example, in comparing three standard covariance structures, spherical, Gaussian, and exponential,

it was determined the exponential covariance structure performed best based on AIC and BIC. The parameters ξ and τ x were estimated based on restricted maximum likelihood and estimated generalized least squares as ˆξ = 27.89 and τˆx = [186.72, 1.42]. Let B correspond to the 67 counties in Florida. Let X(B) = [X(1), X(2), . . . , X(67)] . The block kriging distribution is X(B)|X(so ), ξ, τ x ∼ N(W = 167×1 ξ+B,so −1 X,oo (X(so )−1n ξ), x = B,B − B,so −1 X,oo so ,B ), where B,B is the covariance among the blocks, B,so is the point-to-block covariance and, X,oo is the covariance matrix among the points. To approximate W, first, a fine grid was laid over state. Let su denote the coordinates of the nu = 3,061 grid locations. Then the conditional mean, E(X(su )|X(so )) = 1nu ξ + X,uo −1 X,oo (X(so ) − 1n ξ), was estimated using the estimated values of ξ and τ x . Fiˆ by averaging the estimated nally, W was estimated as W conditional mean values within each county. The block kriging variance matrix was calculated in a similar fashion. See Schabenberger and Gotway (2005) for more details regarding block kriging. We modeled the number of cases, Y(B) = [Y (1), Y (2), . . . , Y (67)] using Poisson regression with the natural logarithm of the expected number of cases E(B) = [E(1), E(2), . . . , E(67)] as an offset and an effect due to PM2.5 . Because X(c) is not observed, we instead considered the model Y(B) ∼ Poisson(μ) log(μ) = log(E(B)) + β0 167×1 + β1 W + β1 x x ∼ N(067×1 , x ). To estimate the parameters in the model, we used the estiˆ and x ,  ˜ x , calculated using the estimated mates of W, W, values of ξ and τ x in the pseudo-penalized quasi-likelihood algorithm to account for the additional term β1 x and the estimation of ξ. In addition, we considered the naive model that ˆ as the true covariate. treats W For the naive approach, the estimate β1 was 0.019 with a variance and standard error of 5.3E−6 and 0.0023, respectively. When using pseudo-penalized quasi-likelihood, the estimate of β1 was 0.018 with a variance and standard error of 8.5E−6 and 0.0029, respectively. Accounting for the additional uncertainty due to kriging resulted in a 60% increase in the variance of the estimate of β1 , 5.3E−6 to 8.5E−6, and a 27% increase in the standard error, 0.0023 to 0.0029. The parameter estimate decreased 8% from 0.019 to 0.018. A Wald test for the significance of the parameter resulted in a z-score of 8.41 in the naive case and 6.09 in the PPQL case, a 27% decrease. Although the effect of PM2.5 was significant using the naive approach and the PPQL approach advocated here, the increased uncertainty associated

PPQL for Spatial-Misalignment with Non-Normal Data with PPQL estimator reflects the additional uncertainty associated with using the block kriging mean as the covariate. Based on the PPQL estimate of β1 , during the time period under consideration, a unit increase in PM2.5 results in an estimated multiplicative effect of 1.018 on the rate of asthma cases with a 95% confidence interval of (1.012, 1.024). 4. Simulation Study As noted earlier, when using a GLM with the canonical link function, the estimates of β1 are unbiased for Poisson distributed responses and biased for binomially distributed responses. In both cases, the standard errors are biased. To address whether these biases are of practical concern when interest lies in valid inference for β1 , a simulation study was conducted for the point-to-point misalignment problem. The spatial domain is [0,500] × [0,500] ⊂R2 . The response is observed at 100 locations and the covariate is observed at 100 different locations. The covariate surface was generated from X(s) ∼ Gaussian Process(ξ = 40, KX (·)), where KX (s, s ; τ x ) is an exponential covariance function defined for any two locations s and s parameterized by an unknown vector of 2 parameters τ x = [scale, range]. We consider two distributions for the response and use the canonical link: (1) Poisson (log link) and (2) Binomial (logit link). For Poisson and logistic regression, the usual GLM weights depend on the mean parameters β0 and β1 . The weights derived using the PPQL algorithm depend on β0 , β1 , the scale and the range. We vary β0 , β1 , the scale, the range and the Berkson parameters: • • • •

β1 : 4 levels: 0.02, 0.05, 0.07, 0.10. Scale: 2 levels: 35, 70. Range: 3 levels: 50, 200, 400. Berkson Parameters: 2 levels: known and unknown.

For the Poisson distribution, 14 levels of β0 were considered and for the Binomial distribution, 8 levels of β0 were considered: • Poisson: defined as β0 = log(M) − β1 × 40: 14 levels of M: 1, 2, 4, 8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48. log(M) − β1 × 40: 8 levels • Binomial: defined as β0 = 1 − log(M) of M: 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8. In addition, for the binomially distributed response we considered three levels of the number of trials at each location: 1, 10, 50. Crossing the factors leads to a total of 4 × 2 × 3 × 2 × 14 = 672 settings for the Poisson distribution, and 4 × 2 × 3 × 2 × 8 × 3 = 1,152 settings for the Binomial distribution. For each of the 1,824 simulation settings, 1,000 simulated datasets were generated. 4.1. Study Results Because the results are qualitatively similar for known and unknown Berkson parameters, all figures and tables correspond to the settings where the Berkson parameters were known. Although only a subset of results are presented here, the remaining simulation results support the conclusions drawn. In each ¯ table, the mean (βˆ1 ), variance (ˆ σˆβ2 ), bias, and mean squared 1

655

error associated with the 1,000 estimates of β1 are provided. ˆ1 (¯s2 ) In addition, the average of the estimated variances of β ˆ β1 and the proportion of the 95% confidence intervals covering the true value are provided.

4.1.1. Poisson. For both the naive method and the PPQL method, the estimates of the regression parameter β1 were virtually indistinguishable, regardless of the values of the other factors. For both methods, the biases were negligible. For the naive method, the variance estimator underestimated the true variability. The effect of each factor given all other factors fixed was explored. The magnitude of the bias increased as the expected value of Y at X = 40 increased, as the slope increased, as the range decreased, or as the scale increased (Table 1). In addition, the variability of the estimator was larger when the Berkson parameters were estimated. Because the naive estimate of the variability was downwardly biased, coverage for 95% confidence intervals was below 95% and followed the same decreasing relationship as the bias in the variance estimator (Figure 2). For the PPQL method, the variance estimator approximated the true variability well. The variability was larger when the Berkson parameters were estimated. Biases were negligible and coverage was near 95% in all scenarios (Figure 3 and Table 1).

4.1.2. Binomial. For 1 trial, the naive and PPQL estimates of the slope and its variance were not different. Biases in both the parameter estimate and its variance did not affect the performance of 95% confidence intervals. That is, coverage for 95% confidence intervals was near 95% for all settings (results not shown). For 10 trials, in almost all settings, the naive and PPQL estimates of the slope and its variance were not different and biases did not affect coverage (Table 2). The only exception was the case when the slope was equal to 0.1 and the scale was equal to 70. In this case, the naive variance estimate was downwardly biased and coverage for 95% intervals was between 80% and 90%. The PPQL method performed well in this case, and coverages were near 95% (Table 2). With 50 trials, for the naive method, the bias in estimates of β1 were negligible. However, the associated estimates of the variance underestimated the true variability. When the expected value of Y at X = 40 was equal to 0.5, the bias in the variability reached a maximum and decreased as the expected value of Y at X = 40 moved away from 0.5. This was expected based on the symmetry of the variance function around 0.5. With respect to the other variables, the same relationships discussed for the naive variance in the Poisson settings were observed. As a result, coverage for 95% confidence intervals is below 95% (Table 2). The PPQL-based estimates of β1 were essentially unbiased. Estimates of the variability were downwardly biased when the range was 50 and upwardly biased when the range was 200 or 400. The bias decreased as the expected value of Y at X = 40 increased. However, because the bias was relatively small, coverage for 95% confidence intervals was near 95% for all settings (Table 2).

656

Biometrics, September 2014 Table 1 Poisson response simulation results

Method Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL

β1

Range

E(Y|X=40)

¯ ˆ1 β

Absolute bias

MSE

σˆˆβ2

¯ ˆsˆβ2

Coverage

0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10

50 50 50 50 400 400 400 400 50 50 50 50 400 400 400 400

1 1 48 48 1 1 48 48 1 1 48 48 1 1 48 48

0.020 0.020 0.020 0.020 0.019 0.020 0.020 0.020 0.097 0.098 0.096 0.100 0.099 0.100 0.100 0.100

2.50E−06 3.54E−04 1.11E−04 1.61E−04 5.19E−04 3.44E−04 2.42E−05 1.54E−04 2.94E−03 1.81E−03 3.57E−03 9.27E−05 1.02E−03 3.86E−04 9.80E−05 1.35E−04

3.45E−04 3.49E−04 1.25E−05 1.14E−05 4.65E−04 4.65E−04 1.07E−05 1.05E−05 4.38E−04 3.99E−04 2.09E−04 9.33E−05 6.36E−04 6.38E−04 6.49E−05 4.29E−05

3.45E−04 3.49E−04 1.25E−05 1.14E−05 4.64E−04 4.65E−04 1.07E−05 1.05E−05 4.30E−04 3.95E−04 1.96E−04 9.33E−05 6.35E−04 6.37E−04 6.48E−05 4.28E−05

3.25E−04 3.36E−04 6.71E−06 1.14E−05 4.93E−04 4.97E−04 1.01E−05 1.17E−05 2.50E−04 4.15E−04 5.17E−06 9.57E−05 5.70E−04 6.23E−04 1.16E−05 4.01E−05

0.943 0.950 0.856 0.949 0.957 0.957 0.948 0.972 0.862 0.950 0.272 0.941 0.939 0.959 0.535 0.950

1

1

Known Berkson parameters: slope equal to 0.02 and 0.1; range equal to 50 and 400; scale equal to 70; E(Y |X = 40) equal to 1 and 48; ¯ ˆ1 and coverage. ˆˆ2 of the estimates of β1 ; average of the estimated variance ¯sˆ2 of β average βˆ1 , MSE, and variance σ β1

5. Discussion For linear models and generalized linear models with spatially misaligned data, if Berkson error is induced from using kriging to align the disparate datasets, then the REMLbased estimated generalized least squares estimator described in Lopiano et al. (2013) and the pseudo-penalized quasilikelihood estimator described here can be used to conduct inference. In practice, however, the parameters that define the kriging mean and variance, that is, the Berkson parameters, are unknown. A limitation of using estimated Berkson parameters is estimating λ introduces an additional source of uncertainty that can cause bias in parameter estimates and standard errors to be too small. If the Berkson parameters can be estimated with sufficient precision, then the approach presented here should be adequate. If not, bootstrap algorithms can be used in conjunction with the method presented in Section 2 in an attempt to account for the additional error. For complicated data-generating mechanisms, deriving a bootstrap procedure is not trivial. In the context of measurement error models, see Buonaccorsi (2009, Section 16.6) for a discussion of bootstrap methods. In our work, we only considered PPQL to account for the additional uncertainty in the generalized linear model. For GLMM, the PQL approach can perform poorly, particularly when estimating variance components (Agresti, 2002, Section 12.6.4). However, because the variance components λ are estimated using an external dataset and used in the PQL Fisher scoring algorithm as if they are known, the resulting pseudoPQL approximation for estimating the main effects β does not use variance components that are estimated via restricted penalized quasi-likelihood methods. As a result, as illustrated in our simulation results, the PPQL method performs well when conducting inference. Nonetheless, exploring the performance of other methods developed for GLMMs would be of interest. We showed, in the context of both point-to-point and pointto-areal spatial misalignment, because the unobserved covari-

β1

ate has a Berkson error relationship with the kriging mean, the methods can be applied to misaligned data in linear and generalized linear models. Although we accounted for estimating the mean parameters ξ, the additional uncertainty in estimating τ x was not accounted for here. To account for this in linear models, Szpiro et al. (2011) developed three bootstrap methodologies. Although bootstrap methods could be developed for non-normal responses, in the simulation study considered here, the additional uncertainty from estimating the Berkson parameters proved negligible. The simulation results also illustrate naive inference in Poisson regression models with Berkson error can perform poorly depending on the size of the effect, the size of the intercept, and the covariance matrix associated with the Berkson error. Consequently, incorrect inferences may be drawn when using Poisson regression in the presence of Berkson error. The result has implications for environmental association studies where relationships among misaligned data are routinely assessed using methods that induce Berkson error. For binomially distributed responses, the simulation results show as the number of trials increases, the effect of Berkson error becomes more pronounced. For the spatial misalignment problem, the performance of naive methods in Poisson and logistic regression depends on the parameters β0 , β1 , the range, and the scale. We illustrated coverage for 95% confidence intervals fell below the nominal level when, for all other values fixed, we increased the magnitude of β0 , increased the magnitude of β1 , decreased the range, and increased the scale. In practice, however, without knowing the true parameters, the true coverage of naive confidence intervals will be unknown and we recommend using an adjusted procedure to conduct inference. Our approach builds upon and provides a new approach in a recent line of work in the area of spatial misalignment (Madsen et al., 2008; Gryparis et al., 2009; Szpiro et al., 2011; Lopiano et al., 2013). First, we extend the basic ideas of these

PPQL for Spatial-Misalignment with Non-Normal Data

657

Known Berkson Parameters

1.0 0.8 0.6 0.4 10

Range=50 Range=200 Range=400

0.0

0.0

Range=50 Range=200 Range=400 20

30

40

0

10

20

30

40

E(Y|X=40)

Beta1=0.1, Scale=35

Beta1=0.1, Scale=70

0.8 0.6 0.4 0.2

0.2

0.4

0.6

0.8

1.0

E(Y|X=40)

1.0

0

Naive Coverage

Beta1=0.02, Scale=70

0.2

0.2

0.4

0.6

0.8

1.0

Beta1=0.02, Scale=35

0

10

Range=50 Range=200 Range=400

0.0

0.0

Range=50 Range=200 Range=400 20

30

E(Y|X=40)

40

0

10

20

30

40

E(Y|X=40)

Figure 2. Poisson response naive coverage: For known Berkson parameters, a plot of the coverage for 95% confidence intervals based on the naive method versus the expected value of Y at X = 40 with a separate color for each value of the range parameter. Each cell corresponds to a combination of the slope and scale parameter. The bold line is drawn at the nominal 95% coverage level.

658

Biometrics, September 2014

Known Berkson Parameters

1.0 0.8 0.6 0.4 10

Range=50 Range=200 Range=400

0.0

0.0

Range=50 Range=200 Range=400 20

30

40

0

10

20

30

40

E(Y|X=40)

Beta1=0.1, Scale=35

Beta1=0.1, Scale=70

0.8 0.6 0.4 0.2

0.2

0.4

0.6

0.8

1.0

E(Y|X=40)

1.0

0

PQL Coverage

Beta1=0.02, Scale=70

0.2

0.2

0.4

0.6

0.8

1.0

Beta1=0.02, Scale=35

0

10

Range=50 Range=200 Range=400

0.0

0.0

Range=50 Range=200 Range=400 20

30

E(Y|X=40)

40

0

10

20

30

40

E(Y|X=40)

Figure 3. Poisson response PPQL-based coverage: For known Berkson parameters, a plot of the coverage for 95% confidence intervals based on the PPQL method versus the expected value of Y at X = 40 with a separate color for each value of the range parameter. Each cell corresponds to a combination of the slope and scale parameter. The bold line is drawn at the nominal 95% coverage level.

PPQL for Spatial-Misalignment with Non-Normal Data

659

Table 2 Binomial response simulation results Method Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL Naive PPQL

Trials

β1

Range

E(Y|X=40)

¯ ˆ1 β

Absolute bias

MSE

σˆˆβ2

¯ ˆsˆβ2

Coverage

10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50

0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10

50 50 50 50 400 400 400 400 50 50 50 50 400 400 400 400 50 50 50 50 400 400 400 400 50 50 50 50 400 400 400 400

0.10 0.10 0.50 0.50 0.10 0.10 0.50 0.50 0.10 0.10 0.50 0.50 0.10 0.10 0.50 0.50 0.10 0.10 0.50 0.50 0.10 0.10 0.50 0.50 0.10 0.10 0.50 0.50 0.10 0.10 0.50 0.50

0.020 0.021 0.020 0.020 0.019 0.019 0.021 0.021 0.095 0.097 0.094 0.098 0.100 0.100 0.100 0.101 0.020 0.020 0.020 0.020 0.020 0.020 0.020 0.021 0.095 0.099 0.094 0.100 0.099 0.100 0.100 0.101

2.71E−04 6.55E−04 3.32E−05 3.77E−04 1.05E−03 9.08E−04 9.20E−04 1.09E−03 4.91E−03 2.72E−03 6.06E−03 1.99E−03 4.23E−04 4.12E−04 1.83E−04 1.05E−03 9.25E−05 2.79E−04 1.74E−04 1.54E−04 4.67E−04 3.08E−04 4.06E−04 5.58E−04 5.00E−03 1.42E−03 6.25E−03 4.89E−04 9.06E−04 1.69E−04 3.76E−04 5.86E−04

3.91E−04 3.97E−04 1.54E−04 1.54E−04 5.13E−04 5.13E−04 1.97E−04 1.98E−04 5.13E−04 5.22E−04 3.33E−04 3.13E−04 6.11E−04 6.13E−04 2.88E−04 2.90E−04 8.11E−05 8.32E−05 3.43E−05 3.42E−05 1.00E−04 1.00E−04 4.02E−05 4.05E−05 2.42E−04 2.02E−04 1.97E−04 1.29E−04 1.57E−04 1.50E−04 9.63E−05 8.79E−05

3.91E−04 3.96E−04 1.54E−04 1.54E−04 5.12E−04 5.13E−04 1.96E−04 1.97E−04 4.89E−04 5.14E−04 2.96E−04 3.09E−04 6.11E−04 6.13E−04 2.88E−04 2.89E−04 8.11E−05 8.31E−05 3.42E−05 3.42E−05 1.00E−04 1.00E−04 4.01E−05 4.02E−05 2.17E−04 2.00E−04 1.58E−04 1.29E−04 1.56E−04 1.50E−04 9.61E−05 8.76E−05

3.65E−04 3.77E−04 1.32E−04 1.40E−04 5.48E−04 5.52E−04 1.95E−04 1.98E−04 3.45E−04 4.98E−04 1.61E−04 2.84E−04 6.49E−04 6.99E−04 2.39E−04 2.80E−04 7.22E−05 7.89E−05 2.64E−05 3.20E−05 1.09E−04 1.11E−04 3.90E−05 4.08E−05 6.86E−05 1.88E−04 3.20E−05 1.31E−04 1.28E−04 1.68E−04 4.75E−05 8.13E−05

0.937 0.945 0.933 0.945 0.955 0.956 0.950 0.952 0.891 0.932 0.824 0.922 0.947 0.964 0.928 0.949 0.929 0.939 0.906 0.946 0.951 0.958 0.945 0.951 0.719 0.929 0.568 0.938 0.918 0.961 0.840 0.951

1

1

Known Berkson parameters: slope equal to 0.02 and 0.1; range equal to 50 and 400; scale equal to 70; E(Y |X = 40) equal to 0.1 and 0.5; ¯ ˆ1 and coverage. ˆˆ2 of the estimates of β1 ; average of the estimated variance ¯sˆ2 of β average βˆ1 , MSE, and variance σ β1

approaches to the case of generalized linear models. Second, although the uncertainty associated with estimating the mean parameter ξ is accounted for in all of this work, that associated with the estimation of τ x is not. Thus, in this article we investigate, through simulation, the effect of the additional uncertainty associated with estimating τ x . Third, the method presented here is fast and general enough to allow for low rank models of the covariance among the X values necessary to extend the method to large datasets (Cressie and Johannesson, 2008). This, plus its relative simplicity, make it a viable alternative to researchers and practitioners without expertise or time to develop complex Markov chain Monte Carlo (MCMC) algorithms necessary to implement fully Bayesian models in the spatial setting. Nonetheless, Bayesian hierarchical models are an appealing approach to account for Berkson error in spatial misalignment problems because the parameter estimation process would naturally propagate all types of uncertainty through the model. Such an approach can be computationally intensive and, in more complicated models, substantial effort may be required to set up algorithms to approximate posterior

β1

distributions. Although MCMC methods and more specifically the Gibbs sampler have been the standard approach for Bayesian models, integrated nested Laplace approximations (INLA) methods offer a computationally efficient alternative (Fong, Rue, and Wakefield, 2010). In the context of the Berkson error model considered here, such an approach could represent a harmony between the desirable properties of Bayesian inference and the computational efficiency of the proposed simpler approach. 6. Supplementary Materials Web Appendices and supplementary files referenced in Section 3 are available with this paper at the Biometrics website on Wiley Online Library. Software needed to replicate the simulation studies is provided. For confidentiality reasons, the applications presented in the text cannot be shared.

References Agresti, A. (2002). Categorical Data Analysis, 2nd edition. New York, NY: John Wiley and Sons.

660

Biometrics, September 2014

Banerjee, S. and Gelfand, A. E. (2002). Prediction, interpolation and regression for spatially misaligned data. Sankhya, Series A 64, 227–245. Breslow, N. E. and Clayton, D. G. (1993). Approximate inference in generalized linear mixed models. Journal of the American Statistical Association 88, 9–25. Buonaccorsi, J. P. (2009). Measurement Error: Models, Methods and Applications. Boca Raton, FL: Chapman & Hall/CRC. Cressie, N. and Johannesson, G. (2008). Fixed rank kriging for very large spatial data sets. Journal of the Royal Statistical Society, Series B 70, 209–226. Florida Department of Environmental Protection (2007). Health advisory: May 10, 2007. http://www.dep.state.fl.us/ mainpage/em/2007/wildfire/health/Advisory 051007.pdf. Fong, Y., Rue, H., and Wakefield, J. (2010). Bayesian inference for generalized linear mixed models. Biostatistics 11, 397–412. Gelfand, A. E., Zhu, L., and Carlin, B. P. (2001). On the change of support problem for spatio-temporal data. Biostatistics 2, 31–45. Georgia Forestry Commision (2007). The historic 2007 georgia wildfires: Learning from the past-planning for the future. http://www.wildfirelessons.net/documents/Historic 2007 GA Wildfires.pdf. Gong, G. and Samaniego, F. J. (1981). Pseudo maximum likelihood estimation: Theory and applications. The Annals of Statistics 9, 861–869. Gotway, C. A. and Young, L. J. (2002). Combining incompatible spatial data. Journal of the American Statistical Association 97, 632–648. Gryparis, A., Paciorek, C. J., Zeka, A., Schwartz, J., and Coull, B. A. (2009). Measurement error caused by spatial misalignment in environmental epidemiology. Biostatistics 10, 258– 274. Lopiano, K. K., Young, L. J., and Gotway, C. A. (2011). A comparison of errors in variables methods for use in regression models with spatially misaligned data. Statistical Methods in Medical Research 20, 29–47. Lopiano, K. K., Young, L. J., and Gotway, C. A. (2013). Estimated generalized least squares in spatially misaligned regression models with Berkson error. Biostatistics 14, 737–751.

Madsen, L., Ruppert, D., and Altman, N. S. (2008). Regression with spatially misaligned data. Environmetrics 19, 453–467. McCullagh, P. and Nelder, J. A. (1989). Monographs on Statistics and Applied Probability 37: Generalized Linear Models, 2nd edition. Boca Raton, FL: Chapman & Hall/CRC. Mugglin, A. S., Carlin, B. P., and Gelfand, A. E. (2000). Fully model-based approaches for spatially misaligned data. Journal of the American Statistical Association 95, 877–887. National Center for Health Statistics (1998). Classification of diseases and injuries. ftp://ftp.cdc.gov/pub/Health Statistics/NCHS/Publications/ICD-9/ucod.txt. Paciorek, C. J., Yanosky, J. D., Puett, R. C., Laden, F., and Suh, H. H. (2009). Practical large-scale spatio-temporal modeling of particulate matter concentrations. The Annals of Applied Statistics 3, 370–397. Sadler, T. W. (2011). Langman’s Medical Embryology, 12th edition. Baltimore, MD: Wolters Kluwer Health. Schabenberger, O. and Gotway, C. A. (2005). Statistical Methods for Spatial Data Analysis, 1st edition. Boca Raton, FL: Chapman & Hall/CRC. Szpiro, A. A., Sheppard, L., and Lumley, T. (2011). Efficient measurement error correction with spatially misaligned data. Biostatistics 12, 610–623. The State of Florida State Emergency Response Team (2007). Statewide disaster response morning briefing: May 9, 2007. http://www.dep.state.fl.us/mainpage/em/2007/wildfire/ response/20070509 briefing.pdf. Thurston, S. W., Spiegelman, D., and Ruppert, D. (2003). Equivalence of regression calibration methods in main study/external validation study designs. Journal of statistical planning and inference 113, 527–539. U. S. Environmental Protection Agency (2007). Air quality system (aqs) data for downloading. http://www.epa.gov/ ttn/airs/airsaqs/detaildata/downloadaqsdata.htm. Zhu, L., Carlin, B. P., and Gelfand, A. E. (2003). Hierarchical regression with misaligned spatial data: Relating ambient ozone and pediatric asthma ER visits in Atlanta. Environmetrics 14, 537–557.

Received May 2013. Revised February 2014. Accepted March 2014.