Spat Demogr DOI 10.1007/s40980-016-0018-4

Income Estimation Using Night Luminosity: A Continuous Spatial Model Ma´rcio Poletti Laurini1

Ó Springer International Publishing AG 2016

Abstract We propose an estimation methodology of per-capita incomes using satellite information on night luminosity (DMSP-OLS Nighttime Lights Time Series), using a structure of continuous spatial random effects and correction for measurement errors. This methodology allows the construction of income measures for disaggregated units and also building income estimates for years without census information. We apply this methodology using information for the Brazilian territory, and the results indicate that inclusion of spatial random effects is critical to the accuracy of the estimated income, and that these spatial effects are persistent in time. Keywords Income estimation  Nighttime lights  Continuous spacial models  Mate´rn covariance  Geo-referenced data

1 Introduction Measures of product and income are critical in economic and demographic analysis. Welfare analysis, quality of life and human development are dependent on the existence and quality of the used income measurements. These measurements are not only important at the national level, as aggregate measures of Gross Domestic Product, but are also equally important in smaller aggregations, such as states, municipalities or neighborhoods within the same municipality. The construction and evaluation of regional and urban development policies depend on good income measures. Similarly, disaggregated income measures are important to business investment decisions, location of shopping centers, real estate investment, & Ma´rcio Poletti Laurini [email protected] 1

Department of Economics, FEA-RP - University of Sa˜o Paulo, Av. dos Bandeirantes 3900, Ribeira˜o Preto 14040-950, Brazil

123

M. P. Laurini

marketing policies and thus a wide range of economic decisions depend on these disaggregated measurements. However, there is a fundamental limitation in the construction of income measures for smaller economic units. Normally we only have annual GDP measures for countries or states, with the lack of such measures at higher frequencies for municipalities and districts. For example, we have only the direct construction of per capita income measures for municipalities in Brazil, the dataset used in this study, through census information. The demographic Census is conducted in Brazil at an average frequency of 10 years, and the latest measures of per capita income of Brazilian municipalities were built with information from 1991, 2000 and 2010 demographic Census by the United Nations Development Program, Applied Economics Research Institute (IPEA) and the Joa˜o Pinheiro Foundation. Thus, we do not have direct observations of incomes for municipalities in years without census. Another fundamental problem is that the smallest available unit of measure for per-capita income is at municipality level and we have no observations of income for smaller units such as districts or neighborhoods. To construct measures in this scale must perform localized surveys, or use information derived from related measures and subject to severe measurement errors as taxation of real estate. There is no a simple method to get these measures at this level of disaggregation or years without census. In a groundbreaking work, Henderson et al. (2012), based on the works of Elvidge et al. (1997) and Doll et al. (2006), proposed using satellite observations about night luminosity as an income proxy. The fundamental idea is to explore the correlation between luminosity and economic development measures, such as access to electricity and public goods, using a night light measure created by the United States Air Force Defense Meteorological Satellite Program (DMSP), known as Version 4 DMSP-OLS Nighttime Lights Time Series. These measurements allow assessing night illumination on grids of 30 arcsec, covering the longitudes and latitudes of 180 -180 to -65 to 75, encompassing most of the populated areas of the globe. This measure is an integer between 0 (no illumination) and 63, and represents the intensity of night luminosity for a certain geospatial location. Henderson et al. (2012) use this proxy to improve the measurements of income in countries with low quality of official data, and shows that this method allows estimating an adequate measure of income, especially for developing countries. Henderson et al. (2012) also show that this measure of income allows assess important questions about economic development in sub-Saharan Africa. They discuss if regions on the coast grow faster than regions in the interior, and if regions more prone to malaria grow less than regions relatively free of this disease. The quality of this methodology in building income measures is shown by good approximation of economic growth observed in the Korean peninsula, the impact of the Asian crisis in 1990 years in Indonesia and the impact of the genocide in Rwanda in the level of income. Two other very important works analyzing the quality of indicators based on night luminosity data are Chen and Nordhaus (2011) and Nordhaus and Chen (2014). These works discuss the validity of using this measure in the construction of

123

Income Estimation Using Night Luminosity: A Continuous...

economic data, and they point out that these data are informative for countries with low-quality statistical systems, but low informative for countries with good official statistics, analyzing in detail the impact of measurement error problems in the official statistics on the construction of measures using nightime information. A complete description of the Nighttime Lights Time Series project can be found in Cauwels et al. (2014). In this work the authors also analyze the spatial distribution and the dynamic evolution of nighttime lights, and discuss some applications of these data in demographics and politics. They also calculate an interesting measure of center of light from the planet, comparing with the center of economic activity discussed in other studies. They also introduce a measure of spatial Gini index, using the spatial distribution of light as a proxy to calculate this measure of inequality. Other applications in demography and politics can be found in Elvidge et al. (2009), calculating a global poverty rate using this information, infanty mortality (Chen 2015), population modeling in China (Zhuo et al. 2009) and Amazon (Amaral et al. 2006), measurement of city limits (Imhoff et al. 1997) and the impact of military interventions (Agnew et al. 2008). A survey of applications using nightime lights can be found in Cauwels et al. (2014). In our work we start from the same idea of Henderson et al. (2012) on using night luminosity as a proxy for income, but focusing on estimating income for regions in the same country, combining available information of per capita income for municipalities with the night luminosity proxy. Our approach differs from that used in Henderson et al. (2012) by including continuous spatial random effects and corrections for measurement errors in the analysis using a Bayesian estimation methodology. The mechanism of spatial random effects allows capturing the spatial dependence in luminosity and income for each location in the spatial continuum, and this correction is especially important in the analysis due to strong spatial dependence on income within the same region. An important assumption in the methodology of Henderson et al. (2012) is that income estimation errors are not correlated. However, in the presence of spatial dependence this assumption may not be valid, especially if these measures are evaluated for shorter distances. The main focus of Henderson et al. (2012) is to use the night luminosity information to complement the information in official statistics for low quality regions, whereas our goal is to build a suitable methodology for estimation in the presence of spatial dependence and measurement errors, especially important in the estimation of income for disaggregated regions, i.e., neighboring regions. Our work addresses the problem of spatial correlation in errors in Henderson et al. (2012), and provides a methodology for continuous spatial estimation of income, combining the luminosity information with available information for per-capita income. To accomplish this we use a continuous geostatistical model based on equivalence between the weak solution of a stochastic partial differential equation and a Gaussian random field represented by a covariance function of Mate´rn class proposed in Lindgren et al. (2011). This representation allows to represent the income and luminosity processes in the continuous space, using a representation of finite elements using a triangulation of the analyzed space. This formulation allows to use a representation of Gaussian Markov Random Fields to the spatial

123

M. P. Laurini

dependence process, which enables the use of efficient computational algorithms for the representation and inference processes. Our estimation method is based on Bayesian inference, using the integrated nested Laplace approximations (INLA) proposed in Rue et al. (2009), which are accurate analytical approaches for the calculation of the posterior distribution of parameters, predictions and latent factors. This methodology does not need the use of Monte Carlo simulation methods, allowing the analysis of complex problems in reduced computational time. In this paper, we use the proposed methodology for the estimation of per capita income for the Brazilian space continuum, using the information of per capita incomes of Brazilian municipalities for the years 2000 and 2010. We discuss some different ways of implementing the methodology, with and without the inclusion of a measurement error mechanism in the analysis and static and dynamic implementations, using purely spatial and spatio-temporal models. From these methodologies we build predictions of per capita income for the entire Brazilian territorial space and also for states and municipalities, showing the results for the state of Sa˜o Paulo and the city of Sa˜o Paulo. The results indicate that the methodology proposed enables an adequated fitting for the observed incomes and also represents a useful spatial interpolation methodology for demographic and economic analyzes that require data for disaggregated regions. This paper has the following structure. In Sect. 2 we present an overview of night luminosity data used in our analysis. The spatial modeling methodology is described in Sect. 3. The Brazilian data used in our applications are presented in Sect. 4. The results of the purely spatial model are shown in Sect. 5, and the results of the spacetemporal dynamic model are discussed in Sect. 7. The final conclusions are in Sect. 8.

2 Night Luminosity Data The night luminosity data used in this study come from satellites mensurations from the United States Air Force Defense Meteorological Satellite Program (DMSP), through the Operational Linescan System (OLS) Sensor, which records the intensity of light from Earth using a system of satellites that circle the Earth 14 times a day. The satellites record the night luminosity between 20:30 and 22:00 in local time. The archiving if this data comprises the years after 1992, and are available through raster images. This data is filtered to remove contaminants such as forest fires, polar auroras and influence of clouds, and thus the final information mainly comprises light sources caused by human activity. The results reported are annual averages of all luminosity measures observed by the set of satellites that watched every location. See Cauwels et al. (2014) for a detailed description of this system. The data reported exclude high latitude regions, and the observations are reported for latitudes between -65 and 75 degrees due to lower measurement accuracy outside these regions. This exclusion is of little relevance, as only 0.0002 % of the human population lives in this locations. The values are reported with an accuracy of 30 arc-sec, which equates to an accuracy of 0.86 square kilometers at the equator. The reported illumination is an integer between 0 and 63, with 0 indicating absence

123

Income Estimation Using Night Luminosity: A Continuous...

of illumination in the region and 63 the maximum truncation observed. As discussed in Henderson et al. (2012) and Cauwels et al. (2014) these measurements are subject to various forms of measurement error, such as saturation of sensors and decay in the quality of measurement due to the age of each satellite, which are partially corrected by image processing algorithms. General descriptive statistics on these data are available in Henderson et al. (2012), and we will discuss in more detail only the years used in this study in Sect. 4. Figures 1 and 2 show a raster image with the night luminosity information for the years 2000 and 2010. We can note by viewing these figures that the data reflect the common understanding of income distribution in the world, and also the growth in income between 2000 and 2010. The figures clearly show the greater nocturnal intensity in the more developed regions, such as North America and Europe, and the smallest lighting reflects the lower developing regions such as Africa and parts of Latin America. In Sect. 4 we present descriptive statistics on income and luminosity information for the spatial grid used in this work, covering the Brazilian territorial space.

3 SPDE-Mate´rn Model The methodology used in this study allows estimating a structure of random effects defined continuously in space. The fundamental idea is to use computational formulation that allows representing the observations of income and night light in a continuous grid. This formulation complements the methodology used in Henderson et al. (2012), by including the spatial dependence in income and luminosity data,

Fig. 1 DMSP-OLS Nighttime Lights Time Series—2000

123

M. P. Laurini

Fig. 2 DMSP-OLS Nighttime Lights Time Series—2010

and making it possible to estimate and interpolate the income date for any disagreggation of interest. In Henderson et al. (2012) the construction mechanism of the income estimates using the light proxies consists of aggregating the lighting measures for the comparable region. For example, if the goal is to build a measure to a certain location, the methodology consists of summing the night lighting measures for all territorial space of this region and using this aggregation as an explanatory variable in a linear model with income observed in this region as the dependent variable. In this paper we use a finite element representation based on the equivalence between the solution of a stochastic partial differential equation and a spatial covariance function. This representation allows to project the observed values of income and luminosity in a triangulation of the analyzed region. This triangulation allows to represent continuously the income and brightness values in space and time. This mechanism allows incorporating the spatial dependence effects observed in income and night luminosity data, and especially allows for continuous predictions and forecasts of per capita income using night luminosity. One great advantage of this methodology is to avoid existing approximations in the use of models based on vertices (areal models) (e.g., Bivand et al. 2013) in the statistical analysis of spatial processes. Models based on vertices sometimes introduce artificial discontinuities in spatial analysis procedures, and in this case are especially inefficient because they are based on aggregations of light observations in a same area, in this case wasting the elevated resolution of nightime observations of DMSP/OLS program.

123

Income Estimation Using Night Luminosity: A Continuous...

3.1 Representation Using Gaussian Markov Random Fields This estimation/representation methodology used in this work is based on Gaussian Markov Random Fields (GMRF) structures, e.g., Rue and Held (2005) and Arbia (2006). This random field can be defined by a graph structure, where an undirected graph is defined by a tuple G ¼fV; Eg with V defining a set of nodes and E a set of fi; jg vertices, with i; j 2 V. In this representation a random vector x ¼ ðx1; x2 ; . . .; xk Þ> 2 Rn defines a Gaussian Markov Random Field with respect the graph G ¼fV; Eg, with expectation EðxÞ ¼ l and precision matrix Q if and only if their density is characterized by a Gaussian density pðxÞ ¼ ð2pÞn=2 jQj1=2 expð 12 ðx  lÞ> Qðx  uÞ, and Qij 6¼ 0 () fi; jg 2 G_i 6¼ j. The GRMF structure is usefull because we are interested in characterizing a continuous version of the Markov property, which requires each value observed in a certain location is conditionally independent given a local neighborhood. We can characterize this property in the continuous space using a spectral representation of spatial dependence function (Simpson et al. 2012) using a fundamental result obtained by Rozanov (1977) on the equivalence of certain random fields and covariance functions. A computationally efficient implementation of this result was obtained by Lindgren et al. (2011), using the results of Whittle (1954) and Rozanov (1977) and directly exploring the connection between stochastic partial differential equations and spatial covariance functions. The fundamental contribution of Lindgren et al. (2011) is to explore the relationship between a covariance function defined on a certain random field and the weak solution of the Stochastic Partial Differential Equation (SPDE): ðj  DÞa=2 xðuÞ ¼ WðuÞ;

u 2 Rd ;

a ¼ m þ d=2;

j [ 0;

m[0

ð1Þ

where ðj  DÞa=2 is a pseudo-differential operator, related to the Laplacian D given by: D¼

d X o2 ox2i i¼1

with W(u) is an innovation process corresponding to a spatial Gaussian white noise with unit variance, defining a covariance function of Mate´rn class. One advantage of the approach proposed in Lindgren et al. (2011) is that you can use a numeric representation of the solution of this SPDE through finite element methods (e.g., Brenner and Scott (2007) on a triangulation, where the random field can be written as: xðuÞ ¼

n X

wk ðuÞxk

ð2Þ

i¼1

where n is the number of vertices on triangulation, from a number of bases wk ðuÞ and weights xk , where these weights are distributed in a Gaussian form. The

123

M. P. Laurini

weights determine the value of the random field at each vertex, and within the triangle values are obtained by interpolation. This basis expansion allows for a continuous approach to discretely observed spatial data, and so avoids the usual distinction between areal and geostatistical models, e.g., Bivand et al. (2013). With this representation we need not aggregate the luminosity data to estimated the income measures, and thus exploring all disponible information to project expected incomes using the full resolution of DSP/ OLS nightime data. We use this covariance function through a hierarchical formulation of a continuously indexed random field Y(s): Y ðsÞ ¼ fyðsÞÞ : ðsÞ 2 D  R2  R where s defines an array of locations whose spatio-temporal covariance function ððsÞðs0 ÞÞ ¼ r2 CððsÞðs0 ÞÞ is set for each ðsÞ and ðs0 Þ  R2  R. We assume that this process is spatially stationary and thus the covariance depends only on the distance between locations, using a distance vector h ¼ ðs  s0 Þ. The hierarchical representation of spatial model is written as: yðsÞ ¼ zðsÞb þ nðsÞ þ eðsÞ nð s Þ ¼ x ð s Þ and z(s) denotes a set of covariates observed in locations s, n the latent component of spatial random effects, and e a component of idiosyncratic innovations unrelated to the spatial effects. In this specification, we have an observation equation: yðsi Þ ¼ zi b þ nðsi Þ þ eðsÞ with eðsÞ  Nð0; r2e Þ, and the variance r2e representing the so called (nugget effect). In this representation we have nðsÞ ¼ xðsÞ, and thus we use the Mate´rn covariance function for xðsÞ, characterizing the random field:    ð3Þ Cov xðsi Þ sj ¼ r2x CðhÞ with its covariance given by: CðhÞ ¼

1 ðkhÞm Km ðkhÞ CðmÞ2m1

with Km a modified Bessel function of the second kind. This model can be rewritten in compact notation as: y ¼ zb þ n þ e n¼x

ð4Þ

~ with yi ¼ ðyi ðs1 Þ; . . .; yðsd ÞÞ0 , with eðsÞ  Nð0; r2e Id Þ and xt  Nð0; R ¼ r2x RÞ, 0 0 ~ defines a zi ¼ ðzðs1 Þ; . . .; zðsd ÞÞ , ni ¼ ðnðs1 Þ; . . .; nðsd ÞÞ . In this formulation R  covariance matrix of d dimension, with elements C jjsi  sj jj .

123

Income Estimation Using Night Luminosity: A Continuous...

The to be estimated in this representation is given by the set  vector of parameters  h ¼ b; r2e ; r2x ; s; j , and their posterior distribution can be written as: pðh; njyÞ / pðyjn; hÞpðhÞ where we assume independent priors for pðhÞ. The yi are conditionally independent due to Markov property, and thus we can write this posterior distribution as: ! T Y pðh; njyÞ / pðyjn; hÞ ðpðnjhÞÞpðhÞ t¼1

using the Gaussian property this posterior is given by: pðh; njyÞ ¼

ðr2e ÞdT=2

! dimðhÞ T Y 1X 0 exp  2 ð y  zb  nÞ ð y  zb  nÞ pðhi Þ re t¼1 i¼1

The Gaussian Markov Random Field representation permits the use of the integrated nested Laplace approximations proposed in Rue et al. (2009), which allow the use of accurate analytical approximations to calculate the posterior distribution of parameters and latent factors in the general class of problems that can be written as GRMFs. In the estimation procedure we assume a set of independent priors with log-Gamma distributions for the parameters of the spatial covariance function, Gaussian distributions for the parameters of the conditional mean and autoregressive parameters and Gamma distribution parameters to accurately. The values used in priors and further details on the estimation procedure can be obtained with the author, and it is importante to note that the results are robust to values used in priors. The approximation of GMRF using the projection on the triangulated mesh, is a global approximation to the true random field. Using the results of Simpson et al. (2012), it’s possible to show the convergence of the approximated field with convergence order OðhÞ with h the length of maximum vertex in the approximation grid. The convergence properties of the approximated can be showed using a piecewise continuous function approximation in the Sobolev space H 1 ; defining R R functions f(s) in L2 space with jjf jj2H 1 ¼ j2 X f ðsÞ2 ds þ X rf ðsÞds\1: Assuming L ¼ ðj2  DÞ, for any f 2 H 1 , Simpson et al. (2012) show that R E X ð f ðsÞLð xðsÞ  xh ðsÞÞdsÞ2  ch2 jjf jj2H 1 , where c is a constant and h the length of maximum vertex.

4 Database Description To carry out our analysis, we will combine income information from municipalities in years with census information with the night luminosity (DMSP-OLS Nighttime Lights Time Series). Figure 3 shows an overlap of night luminosity in 2010 with the territorial division of Brazilian municipalities. The first step is to obtain a subset of

123

M. P. Laurini

Fig. 3 DMSP-OLS Nighttime Lights Time Series in 2010—overlap with Brazilian municipal division

light intensity data that only cover Brazilian territory. For this reason we held a cut in the raster file to obtain a rectangle covering the space in question. In Table 1 we present information related to raster files with the luminosity data in years 2000 and 2010. This file shows the size of the grid (dimensions), the resolution, geographic coverage (extent), the referencing system used, and the name of the original files from the DMSP-OLS Nighttime Lights Time Series used in the analysis (names). We use data from satellite F15 for the year 2000 to permit compatible projections for the years 2000–2007, since F14 satellite provides data for the period 1997–2003 only. The data from this satellite also have a higher luminosity compared to the F14 satellite data and a slightly larger number of observations with non-zero luminosity for the region used in this work. For year 2010 we used the information from F18 satellite. Figures 4 and 5 show the images for the years 2000 and 2010. In Table 2 we present some descriptive statistics for the light intensity data. We provide the full intensity measures for the analyzed grid (Total Intensity) and summary of intensities in the municipalities, which are the intensity measures near the geographic centroid of each municipality. We can note that there are a great number of points on the grid with zero values, which reflect the areas not occupied by human settlements in the Brazilian geographical space. You can also notice the growth in the intensity observed between the years 2000 and 2010 analyzed. The income data used include the average monthly income (per capita income) of individuals resident in the 5565 Brazilian municipalities, denominated in Brazilian Real of August 1, 2010, obtained from the Atlas of Human Development in Brazil 2013, which contains information on the Municipal Human Development Index IDHM, constructed using information extracted from the Brazilian demographic

123

Income Estimation Using Night Luminosity: A Continuous... Table 1 Raster files—2000 and 2010 year: 2000 class: RasterLayer dimensions: 4601, 4924, 22655324 (nrow, ncol, ncell) resolution: 0.008333333, 0.008333333 (x, y) extent: -73.4625, -32.42917, -33.65417, 4.6875 (xmin, xmax, ymin, ymax) coord. ref.: ?proj=longlat ?datum=WGS84 ?no_defs ?ellps=WGS84 ?towgs84=0,0,0 data source: C:\Users\laurini\AppData\Local\Temp\R_raster_laurini\raster_tmp_2014-0916_152248_54912_93128.grd names: F152000.v4b_web.stable_lights.avg_vis values: 0, 63 (min, max) year: 2010 class: RasterLayer dimensions: 4600, 4925, 22655000 (nrow, ncol, ncell) resolution: 0.008333333, 0.008333333 (x, y) extent: -73.46667, -32.425, -33.65, 4.683334 (xmin, xmax, ymin, ymax) coord. ref.: ?proj=longlat ?datum=WGS84 ?no_defs ?ellps=WGS84 ?towgs84=0,0,0 data source: C:\Users\laurini\AppData\Local\Temp\R_raster_laurini\raster_tmp_2014-0916_152414_54912_42700.grd names: F182010.v4d_web.stable_lights.avg_vis values: 0, 63 (min, max)

Fig. 4 DMSP-OLS Nighttime Lights Time Series—section analyzed—2000

123

M. P. Laurini

Fig. 5 DMSP-OLS Nighttime Lights Time Series—section analyzed—2010

Table 2 Descriptive statistics—nighttime lights intensity Min.

1st Qu.

Median

Mean

3rd Qu.

Max.

SE

Total intensity—2000

0

0

0

0.49950

0

63

25.67964

Total intensity—2000a

3

5

7

12.45000

12

63

22.81083

Total intensity—2010

0

0

0

0.77380

0

63

25.65833

Total intensity—2010a

4

6

8

13.98000

14

63

22.34726

Municipalities intensity—2000

0

0

0

5.94500

7

63

24.86793

Municipalities intensity—2000a

4

6

8

14.63000

15

63

22.28993

Municipalities intensity—2010

0

0

0

8.3860

10

63

24.6389

5

7

10.5

17.3400

20

63

21.6286

a

Municipalities intensity—2010 a

The measure with excluded zeros. Municipalities intensity is the measured luminosity near the geographic centroid of the municipality

Census of 2000 and 2010 by the consortium of United Nations Development Programme, Applied Economics Research Institute and the Joa˜o Pinheiro Foundation. Figures 6 and 7 show the distribution of this per capita income measure through the territorial division of municipalities in years 2000 and 2010. We can see in these figures the pattern of wealth concentration in the Center-West, Southeast and South regions, and the lowest relative incomes in municipalities of the North and Northeast regions. In Table 3 we show descriptive statistics of the average monthly per capita income in 2000 and 2010 years used in this work. These statistics reflect the concentration and income disparity among Brazilian municipalities.

123

Income Estimation Using Night Luminosity: A Continuous...

Fig. 6 Per capita income—Brazilian municipalities—2000

Fig. 7 Per capita income—Brazilian municipalities—2010

123

M. P. Laurini Table 3 Descriptive statistics—per capita income—Brazilian municipalities—2000 and 2010 Min.

1st Qu.

Median

Mean

3rd Qu.

Max.

SD

Per-capita income 2000

62.6500

173.5000

308.6000

338.5000

463.1000

1760.0000

602.2916

Per-capita income 2010

96.2500

281.1000

467.5000

493.6000

650.6000

2044.0000

686.4039

5 Static Spatial Model We start the analysis through a purely spatial model, based on a cross-section specification with the inclusion of spatial random effects using the methodology proposed in Sect. 3. We discussed first a model without correction for measurement errors, and in Sect. 6 we discuss how to include a classic measurement error mechanism in this spatial model. The specification without the presence of measurement error is given by: log yðsÞ ¼ a þ iðsÞb þ nðsÞ þ eðsÞ nð s Þ ¼ x ð s Þ where log yðsÞ denotes the logarithm of per-capita income in location s, iðsÞ denotes the average light intensity in the grid cell corresponding to the centroid of the municipality, and nðsÞ and eðsÞ has the interpretation of spatial random effects and idiosyncratic innovations, according to the methodology proposed in Sect. 3. The key component is the estimation of the continuous Mate´rn covariance matrix for spatial random effect. As discussed in Sect. 3 the spatial continuous formulation is based on a representation of finite elements in a triangulation for the analyzed region. To this we use the territorial division of municipalities to construct the triangulated mesh. As discussed in Sect. 3.1, the values inside each vertex in the grid are interpolated using the basis functions, and thus projecting the variables in the continuous neighborhood of each location. The continuous projection using the spde-Mate´rn covariance matrix performs a spatial interpolation process, thereby smoothing the values in the grid, conditionally on the designed model and the estimated parameters. Thus, the approximation given by the use of municipality centroid as location of income measure can be justified by the interpolation properties of the methodology used in the analysis. The triangulated mesh used in this work is shown in Fig. 8, and was built using a Delauney triangulation abased on fine grid of values required for proper representation of spatial process. Note that this triangulation there is a difference between internal and external triangulation of the analyzed domains. The external triangulation is necessary to avoid variance problems at the borders, as discussed in Lindgren et al. (2011). In this example we present only the results for the year 2010, but the results for the year 2000 are similar and can be obtained from the author. We estimated the posterior distribution of parameters using the Bayesian INLA methodology, and the

123

Income Estimation Using Night Luminosity: A Continuous... Fig. 8 Mesh for Brazilian municipalities

Table 4 Model without correction for measurement error Mean

SD

.025q

.5q

.975q

Mode

Intercept

5.9341

0.0740

5.7885

5.9337

6.0817

Intensity

0.0067

0.0003

0.0062

0.0067

0.0073

0.0067

Precision

24.5743

0.5815

23.4495

24.5683

25.7352

24.5574

log-s

0.2869

0.0676

0.1542

log-j

-0.5814

0.1077

-0.7926

Marginal Lik.

339.14

n obs

5561

5.9330

0.2869

0.4197

0.2870

-0.5814

-0.3697

-0.5814

results for the model without measurement error are shown in Table 4, with the results for the intercept a; the b parameter relating the observed intensity to the observed income, e the precision parameters for the idiossincratic errors and the s and j parameters of the spatial Mate´rn covariance, presented in the original log scale used the estimation. We can observe that the model estimates an average income of exp (5.9341) = 377.6999 for a zero night intensity, and for each additional unit of night lighting the average log monthly income grows 0.00067. These estimates indicate a symptom of possible measurement error in light intensity, since on these estimates the maximum expected value for per capita income, without the spatial random effects, of expð5:9341 þ 0:0067  63Þ ¼ 576:0532, which would lead to a severe underestimation of incomes for the locations with more intense luminosity. The importance of this spatial effect estimation is shown in Fig. 9, showing the spatial correlation as a function of distance, with the solid line indicating the expected pattern for a Mate´rn function with the estimated parameters reported in

123

M. P. Laurini

Fig. 9 Spatial correlation—SPDE-Mate´rn—without correction for measurement error

Table 4 and the points the observed correlations. We note that there is a strong pattern of spatial dependence in this data. The estimated random field with the continuous spatial random effects can be seen in Fig. 10, that captures the deviations in log scale for each location. We can observe that the spatial random effects recover part of the prediction error, but is not sufficient to recover the expected pattern of income prediction. This can be seen in Fig. 11, which shows the income predictions using using night luminosity over the random effect. We can observe that dispersion of predicted values can not capture the richest regions in the country, indicating that the adjustment mechanism without the inclusion of a correction for measurement errors is not appropriate.

6 Measurement Error To incorporate the presence of measurement error in the analysis we assume, similar to Henderson et al. (2012), a classic mechanism of measurement error, where the variable x can only be observed through a proxy: w¼xþu with u ¼ ðu1 ; . . .; un Þ0 denoting a measurement error vector. As discussed in Henderson et al. (2012), the light intensity presents some evident measurement problems due to the fact that this intensity is limited to 63 and is only observed in integer values, in addition to problems associated with the satellites used in the measurement. Another problem that can be noted is that the noise filtering tends to eliminate

123

Income Estimation Using Night Luminosity: A Continuous...

Fig. 10 Random effect—SPDE-Mate´rn—without correction for measurement error

Fig. 11 Income prediction—model without correction for measurement error—2010

123

M. P. Laurini

intensity to regions with low brightness, which is consistent with the lack of intensity values less than 4 in regions with positive luminosity. To address these problems we will use a classic measurement error mechanism, but using a Bayesian approach to this problem proposed in Muff et al. (2015). We will consider that this error mechanism is independent and a originating from a Gaussian distribution with zero mean and precision su , with the precision denoting the inverse of the variance. In this representation we assume that the measurement error is also independent of the explanatory and the dependent variables. This representation assumes that there is a true variable x, and repeated observations of the variable j from xi has the structure: wij  N ðxi ; su Þ An important result of this structure is that w is used as a proxy of x, the estimation of the conditional relationship of x will be biased. In a linear regression with measurement error the estimation of model y ¼ b0 þ bx w þ e instead of errorfree model y ¼ b0 þ bx x þ e will lead to an biased estimator jbx j\jbx j; which can be easily understood by the fact that the presence of a measurement error u with positive variance of the variance w is greater than the variance of x, and since the estimator bx is obtained as Cov(y, w) / Var(w), it will usually be less than Cov(y, x) / Var(x), since the variance of w is the sum of the variances of x and u. Because we are using a Bayesian estimation mechanism, we adopt the Bayesian approach proposed in Muff et al. (2015), which is also based on integrated nested Laplace approximation to a hierarchical representation of the measurement error process. We present the only results for a single explanatory variable measured with error in a regression problem. The general case with more than one variable may be found in Muff et al. (2015). This representation is based on a hierarchical model with three levels. The first is the level of observation yjv; h1 , which captures the dependence of response observed y in function of unobserved parameter v and hyperparameters h1 . The second level captures the latent factors in terms of a set of hyperparameters h2 . The third level defines the (hyper) priors for the hyperparameters h ¼ ðh1 ; h2 Þ. In this representation the posterior distribution of the parameters v and h is given by: pðv; hÞ / pðyjv; hÞpðvjhÞpðhÞ and inferences R R about the marginal distributions of vi can be obtained by pðvi jyÞ ¼ h vi pðv; hjyÞdvi dh, with vi denoting the vector v with the ith element excluded. As usual these integrals are not analytically known, and some numerical approximation method must be used. In the case where the variable y does not depend on any other independent variable beyond x, we can write the model for covariate x as: x  Nða0 ; sx Þ and thus the complete model is given by:

123

Income Estimation Using Night Luminosity: A Continuous...

EðyjxÞ ¼ bx x x ¼ a0 þ e x w¼xþu ex  Nð0; sx Þ eu  Nð0; su Þ In this representation a0 is a hyperparameter, the latent field v depends only on x and the parameter vector is given by h ¼ ðbx ; sx ; su ; a0 ; h1 Þ, where h1 defines a set of parameters unrelated to the measurement error mechanism. In this form the joint posterior distribution of x and h It is given by: pðx; hjyÞ / pðhÞpðxjhÞpðwjx; hÞpðyjx; hÞ / pðhÞpðxjw; hÞpðwjhÞpðyjx; hÞ as pðxjhÞpðwjx; hÞ ¼ pðxjw; hÞpðwjhÞ. By the assumed mechanism of measurement error, follow that:     1 1 wjh  N a0 ; s1 u þ sx and pðxjw; hÞ / pðx; hÞpðwjx; hÞ n s o su x / exp  ð x  a0 Þ0 ð x  a0 Þ  ð x  wÞ0 ð x  wÞ 2 2 and thus we have:   xjw; h  N ððsa0 þ su wÞðsx þ su Þ1 ; sx þ su

ð5Þ

which allows assessing the posterior distribution of interest to the expected value of x , given by the prediction using the Eq. 5. In our model we combined the representation of spatial random effects model based on the measurement error for the night luminosity variable, where the correction for measurement error is performed according to the representation: log yðsÞ ¼ a þ i ðsÞbme þ nðsÞ þ eðsÞ nð s Þ ¼ x ð s Þ with i denoting the night luminosity adjusted for measurement error. In Table 5 we present the results of SPDE-Mate´rn model with correction for measurement error. In this table the parameter bme captures the conditional effect of intensity on the log of income, and the parameters su and sx are the precisions in the representation of the classic measurement error model. We can see that the conditional effects are quite different from the model without correction for measurement error, especially in bme which is substantially higher than in the former model, as expected by the existing bias mechanism in the presence of a non-

123

M. P. Laurini Table 5 Model with correction for measurement error Mean

SD

.025q

.5q

.975q

Mode

Intercept

6.0725

0.0827

5.9095

6.0721

6.2377

6.0715

Precision

22.5028

0.5029

21.5142

22.5042

23.4913

22.5138

MEC bme

0.0172

0.0014

0.0142

0.0171

0.0204

0.0170

MEC su

0.0010

0.0001

0.0008

0.0010

0.0013

0.0009

MEC sx

0.0248

0.0043

0.0178

0.0243

0.0345

0.0232

log-s

0.3364

0.0671

0.2158

0.3319

0.4785

0.3183

log-j

-0.7269

0.1344

-1.0087

-0.7188

-0.4838

-0.6942

n obs

5561

Marginal Lik.

-1197.43

controlled measurement error. In this specification, we have that the average income for a zero intensity, given the exponential of intercept, would now assume the value of 433.7637, which is relatively close to the unconditional average for income for the year 2010 of 493.6. In this estimation for each additional unit of intensity, log income grows 0.0172, and he maximum prediction is given by exp(6.0725 ? 0.0172 9 63) = 1281.902, closer to the expected value for the richer regions and more intense night luminosity. In this specification the spatial random effects remain fundamental. This can be seen by the spatial correlation, correspondent to the estimated spatial covariance function parameters, shown in Fig. 12. We can observe a pattern of strong spatial dependence, similar to that observed in the model without control for measurement

Fig. 12 Spatial correlation—SPDE-Mate´rn—with correction for measurement error

123

Income Estimation Using Night Luminosity: A Continuous...

Fig. 13 Random effects—model with correction for measurement error—2010

error. The random field for the spatial random effects is shown in Fig. 13, and the corresponding predictions matching this field with the conditional mean function are shown in Fig. 14. We can observe that the model with the inclusion of correction for measurement error can explain more adequately the incomes for the richest regions, matching the highest light intensities. As a measure of adjustment of the model to the observed data, we show in Tables 6 and 7 some fit measures for the two models, with the correlations between the data for observed income in 2010 for the municipalities and the fitted values for the models with and without measurement errors. In these tables the notation is as follows—WE is the conditional adjustment (regression effects predicting income using ligth intensity) without random effects and without correction for measurement error; WE ? SRE the fit adding WE plus the spatial random effects. ME— conditional adjustment (regression effects) of the model with correction for measurement error, ME ? SRE—the ME plus the spatial random effects. We can observe that the spatial random effects are critical in the model, since the correlation of the models with and without random effects are quite different. For example, the correlation of ME fit without the spatial random effects with per capita income observed in 2010 is estimated as 0.40, while with the sum of the spatial random effects the correlation is about 0.89. Looking at the fit measures between the model fit and observed incomes presented in Table 7, with the mean error (ME), root mean square error (RMSE), mean absolute error (MAE), mean percentage error (MPE) and mean absolute percentage error (MAPE), we can see some important features. When we look at the

123

M. P. Laurini

Fig. 14 Income prediction—model with correction for measurement error—2010 Table 6 Estimated correlations

Year 2010 WE ? SRE WE ME ? SRE ME

y 2010

WE ? SRE

WE

ME ? SRE

ME

1

0.888

0.403

0.890

0.401

1

0.458

0.992

0.452

1

0.419

0.994

1

0.412 1

WE—mean effect (conditional regression effects) of model without correction for measurement errors. WE ? SRE—Mean effects plus the spatial random effects. ME—mean effect (conditional regression effects) of model with correction for measurement errors. ME ? SRE—mean effects plus the spatial random effects

models without random effects, we can see that the model without the correction for measurement error undervalues incomes with an average bias of -92.49, while the model with correction for measurement error overestimates incomes in 23.64. With the inclusion of spatial random effects both models have a bias near the value of -11. In general, the two models with the inclusion of spatial random effects obtain a proper fit, with an mean percentage error of about -1.9 % for the observed incomes. Note that in the model without control for measurement errors, we are correcting the possible measurement problem through the spatial random effects, which is inappropriate because the model with the correction allows to estimate correctly the conditional relationship between luminosity and income.

123

Income Estimation Using Night Luminosity: A Continuous... Table 7 Adjustment measures ME

RMSE

MAE

MPE

MAPE

WE ? SRE

-11.014

112.330

71.801

-1.909

14.875

WE

-92.497

247.647

194.582

-22.083

48.160

ME ? SRE

-11.588

111.584

71.923

-1.919

14.958

23.643

233.616

191.378

2.625

37.782

ME

WE—mean effect (conditional regression effects) of model without correction for measurement errors. WE ? SRE—mean effects plus the spatial random effects. ME—mean effect (conditional regression effects)of model with correction for measurement errors. ME ? SRE—mean effects plus the spatial random effects

7 Spatio-Temporal Model An alternative way to implement this approach is through a space-temporal model, combining the information from more than 1 year in the construction of the income proxy using night luminosity. For this we will use a space-temporal model where we introduce an autoregressive dynamic to the spatial random effects, similar the model proposed in Cameletti et al. (2013). This dynamic model is given by the following representation: yðsÞt ¼ zðsÞt b þ nðsÞt þ eðsÞt nðsÞt ¼ anðsÞt1 þ xðsÞt

ð6Þ

where now the parameter a is an autoregressive parameter for the random field that defines the spatial random effects. The posterior distribution of the parameters h, which now include the autoregressive parameter a, is given by: ! ! T T Y Y pðh; njyÞ / pðyjn; hÞ pðn1 jhÞ pðnt jnt1 ; hÞ pðhÞ t¼1

t¼2

Using the Gaussian property to this random field the posterior distribution can be written as: ! T 1X 0 2 dT=2 pðh; njyÞ ¼ ðre Þ exp  2 ðyt  zt b  nt Þ ðyt  zt b  nt Þ re t¼1  2 d=2  2 rx ~ 1 ~ 1=2 exp 1  a n0 Rn  j Rj 1  a2 2r2x 1   dðT1Þ=2 ðT1Þ=2 1 0~ ~ R n  r2x jRj exp ð n  an Þ ð  an Þ t1 t t1 2r2x t 

dim ðhÞ Y

pðhi Þ

i¼1

123

M. P. Laurini

We also modified the spatio-temporal model to include the classic measurement error mechanism as discussed in the previous section. For this, we unite the representation of spatio-temporal model with the measurement error correction, replacing zt for zt , adapting the notation presented in Sect. 6. Due to limitations of the existing data, we can only include in this specification the years 2000 and 2010, as we only have luminosity observations for 1992 forward, and the data of municipal income per capita for this period only cover the Census years of 2000 and 2010, representing a model with two periods. The estimated parameters for this model are shown in Table 8. We note that the parameters of the conditional mean are similar to those obtained by the model with correction for measurement error, and we can also see that the random effects are quite persistent, with an autoregressive parameter with posterior mean of of 0.994. To check the qualities of fit of this model, we presented in Tables 9 and 10 the measures of correlation and adjustment between observed income in 2000 and 2010, using the results of the models with and without the inclusion of spatial random effects. We can see that the results are quite similar to those obtained with the purely spatial model. This consistency results can be explained by the high persistence of random spatial effects. This characteristic can be seen in Figs. 15 and 16 showing the random fields for spatial effects in 2000 and 2010 adjusted by this model. The graphs show that the spatial random effect is virtually constant between years. This feature is especially important for the construction of income measures for years without census information. As the persistence of random field is close to one, it serves as a good estimate of expected spatial effects for years not observed, assuming a martingale structure. As an example, we estimate income predictions for the years 2001 and 2005 using this method, which are shown in Fig. 17. The income predictions for years with Census information of income (years 2000 and 2010) are showed in Figs. 18 and 19.

Table 8 Spatio-temporal model—2000 and 2010 Mean

SD

.025q

.5q

.975q

Mode

Intercept

5.962

0.8933

4.1432

5.9631

7.7774

5.9648

Precision

19.1988

0.2464

18.7595

19.1813

19.7226

19.1326

MEC bme

0.0200

0.0015

0.0167

0.0200

0.0233

0.0199

MEC su

0.0010

0.0001

0.0008

0.0009

0.0012

0.0009

MEC sx

0.0243

0.0030

0.0193

0.0239

0.0311

0.0231

log-s

0.1062

0.0398

0.0204

0.1096

0.1758

0.1196

log-j

-1.8381

0.1642

-2.1482

-1.8430

-1.5036

-1.8561

0.9944

0.0020

0.9896

0.9948

0.9972

0.9954

n obs

11122

AR(1)-a Marginal Lik.

123

-1759.78

Income Estimation Using Night Luminosity: A Continuous... Table 9 Correlations—spatio-temporal model Year 2000 Year 2000

ME ? SRE

1

ME ? SRE

ME

0.882

0.417

1

0.436

ME

1 Year 2010

Year 2010

ME ? SRE

1

ME ? SRE

ME

0.890

0.399

1

0.434

ME

1

ME—mean effect (regression effects). ME ? SRE—mean effects plus the spatial random effects

Table 10 Adjustment—spatio-temporal model ME

RMSE

MAE

MPE

MAPE

ME ? SRE 2000

-11.094

91.604

60.577

-2.796

18.820

ME 2000

114.853

223.663

179.822

24.139

39.632

ME ? SRE 2010

-9.859

111.152

71.889

-1.923

14.852

ME 2010

-4.125

243.197

194.679

-4.851

40.874

ME—mean effect (regression effects). ME ? SRE—mean effects plus the spatial random effects

7.1 Disaggregated Income Measures As applications of space-temporal model, we will present some disaggregated measures of income, using the estimation results reported in Table 8. These same disaggregated estimation are also constructed using the purely spatial model and are essentially equivalent, and for space reason will not be presented. We first present disaggregated income estimation results for state of Sa˜o Paulo in Brazil, and following the results for the municipality of Sa˜o Paulo. Figure 20 presents the spatial random effects estimated for the State of Sa˜o Paulo. We note that the random effects are consistent with the wealth patterns observed in the state, showing, for example, positive corrections to the regions of Campinas, and negative corrections for the poorest region of the state located in the south border with the Parana´ state. The income estimates obtained by combining the spatial random effects and the intensity effect for the years 2000 and 2010 are shown in Figs. 21 and 22. To facilitate the visualization of the growth pattern during this period, we show in Fig. 23 the income growth rate between 2000 and 2010 fitted by the spatio-temporal model. In Fig. 24 we show a state map with the major cities in the state to facilitate the interpretation of results. To a show a further disaggregated income measure, in Fig. 25 we show the estimation of income for the Municipality of Sa˜o Paulo, using the values projected

123

M. P. Laurini

Fig. 15 Random effects—spatio-temporal model with correction for measurement error—2000

Fig. 16 Random effects—spatio-temporal model with correction for measurement error—2010

123

Income Estimation Using Night Luminosity: A Continuous...

Fig. 17 Income prediction—years 2001 and 2005

Fig. 18 Income prediction—spatio-temporal model—2000

123

M. P. Laurini

Fig. 19 Income prediction—spatio-temporal model—2010

Fig. 20 Random effect—state of Sa˜o Paulo in 2000

123

Income Estimation Using Night Luminosity: A Continuous...

Fig. 21 Income prediction—state of Sa˜o Paulo in 2000

Fig. 22 Income prediction—state of Sa˜o Paulo in 2010

by the spatio-temporal model for the year 2010. In this figure we also show a graph with the division of subprefectures of this city. We can observe that the pattern of income predicted by the model is fairly consistent with the stylized facts about the

123

M. P. Laurini

Fig. 23 Growth—state of Sa˜o Paulo 2000–2010

Fig. 24 Major cities—state of Sa˜o Paulo

distribution of income in Sa˜o Paulo, with higher incomes in the central regions (Se´, Vila Mariana, Santana, and Mo´oca) and lower incomes in peripheral regions, such as Anhanguera, Pirituba and Itaquera, and especially for Parelheiros region. The

123

Income Estimation Using Night Luminosity: A Continuous...

Fig. 25 Income prediction—municipality of Sa˜o Paulo and subprefectures—2010

estimation of income for finer disaggregations and for years without census information shows some of the applications of the methodology proposed in this paper.

123

M. P. Laurini

8 Conclusions In this work we present an estimation methodology for income data using satellite images of night luminosity from the DMSP-OLS Nighttime Lights Time Series program, which explicitly include a continuous spatial component. This spatial component is especially important due to the strong pattern of spatial dependence in income observed in higher disaggregation levels, as municipalities or districts. In the presence of a spatial error component the residual independence assumption assumed in the methodology of Henderson et al. (2012) may be invalid, affecting the validity of income estimates using the night luminosity proxy. Our model is based on a continuous representation of the spatial dependence process, based on the equivalence of a spatial covariance matrix of Mate´rn class and the solution of a stochastic partial differential equation. The continuous representation allows projecting the estimated income using a finite elements projection in the spatial grid, allowing the direct construction of income estimates for any point or area of interest. Our methodology also allows to incorporate the presence of measurement error in the estimation, using a Bayesian formulation for a classic measurement error mechanism. The results indicate that without the incorporation of this correction the night luminosity information leads to an estimate with a negative bias for observed income, especially for regions with low luminosity. We also discussed an estimation method based on a spatio-temporal model, incorporating the information for multiple periods. The results of this model indicate a high persistence for spatial random effects, and allow to project incomes for years without census information using the spatial dependence in the data. We apply the proposed methodology to build income measures to sub-regions of the Brazilian territorial space and show that these estimates are consistent with those observed in incomes at aggregated levels. This method has important applications in demographic and regional economic analysis. It also has direct applications in urban development models and estate planning for areas without direct Census information, such as districts of the same municipality. Further developments of this methodology would be the combination of luminosity data with other satellite informations, such as population density, urban distribution or other processes related to income and economic development. Acknowledgments I thank the valuable comments received from two anonymous referees and the editor Jeremy Porter. I appreciate the support of CNPq and FAPESP. The analysis of this work were performed using the R software (R Core Team 2015), using the INLA, sp, raster, rasterVis, maptools, rgeos, dismo and rgdal packages, and Fig. 3 was created by the GRASS GIS software.

References Agnew, J., Gillespie, T., Gonzales, J., & Min, B. (2008). Baghdad nights: Evaluating the us military ‘surge’ using nighttime light signatures. Environment and Planning, 40, 2285–2295. Amaral, S., Monteiro, A. M. V., Caˆmara, G., & Quintanilha, J. A. (2006). DMSP/OLS night-time light imagery for urban population estimates in the Brazilian Amazon. International Journal of Remote Sensing, 5(10), 855–870.

123

Income Estimation Using Night Luminosity: A Continuous... Arbia, G. (2006). Spatial econometrics—Statistical foundations and applications to regional convergence. Berlin: Springer. Bivand, R. S., Pebesma, E., & Gomes-Rubio, V. (2013). Applied spatial data analysis with R. Berlin: Springer. Brenner, S. C., & Scott, R. (2007). The mathematical theory of finite element methods. Berlin: Springer. Cameletti, M., Lindgren, F., Simpson, D., & Rue, H. (2013). Spatio-temporal modeling of particulate matter concentration through the spde approach. AStA Advances in Statistical Analysis, 97, 109–131. Cauwels, P., Pestalozzi, N., & Sornette, D. (2014). Dynamic and spatial distibution of global nighttime lights. EPJ Data Science, 3, 1–26. Chen, X. (2015). Explaining infant mortality and poverty rates: What we can learn from night-time lights? Spatial Demography, 3, 27–53. Chen, X., & Nordhaus, W. (2011). Using luminosity data as proxy for economic statistics. Procedings of National Academy of Sciences, 108(21), 8589–8594. Doll, C. N., Muller, J.-P., & Morley, J. G. (2006). Mapping regional economic activity from night-time light satellite imagery. Ecological Economics, 57, 75–92. Elvidge, C. K., Baugh, E., Kihn, E., Kroehl, H., & Davis, E. (1997). Mapping of city lights using dmsp operational linescan system data. Photogrammetric Engineering and Remote Sensing, 63, 727–734. Elvidge, C. D., Sutton, P. C., Ghosh, T., Tuttle, B. T., Baugh, K. E., Bhaduri, B., et al. (2009). A global poverty map derived from satellite data. Computer and Geosciences, 35(8), 1652–1660. Henderson, J. V., Storeygard, A., & Weil, D. N. (2012). Measuring economic growth from outer space. American Economic Review, 102(2), 994–1028. Imhoff, M., Lawrence, W., Stutzer, D., & C.D., Elvidge. (1997). A technique for using composite DMSP/ OLS city lights satellite data to map urban areas. Remote Sensing of Environment, 61(4), 361–370. Lindgren, F., Rue, H., & Lindstrom, J. (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach. Journal of the Royal Statistical Society, Series B, 73(4), 423–498. Muff, S., Held, L., Riebler, A., Rue, H., & Saner, P. (2015). Bayesian analysis of measurement error models using integrated nested Laplace approximations. Journal of Royal Statistical Association, Series C, 64(2), 231–252. Nordhaus, W., & Chen, X. (2014). A sharper image: Estimates of the precision of nighttime lights as a proxy for economic statistics. Journal of Economic Geography, 15, 1–30. R Core Team. (2015). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Rozanov, J. A. (1977). Markov random fields and stochastic partial differential equations. Mathematics of the URSS Sbornik, 32, 515–534. Rue, H., & Held, L. (2005). Gaussian Markov random fields: Theory and applications. London: Chapman & Hall. Rue, H., Martino, S., & Chopin, N. (2009). Approximated Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society Series B, 71, 319–392. Simpson, D., Lindgren, F., & Rue, H. (2012). Think continuous: Markovian Gaussian models in spatial statistics. Spatial Statistics, 1, 16–29. Whittle, P. (1954). On stationary processes on the plane. Biometrika, 41(7), 434–449. Zhuo, L., Ichinose, T., Zheng, J., Chen, J., Shi, P. J., & Li, X. (2009). Modelling the population density of china at the pixel level based on DMSP/OLS non-radiance-calibrated nightime light images. International Journal of Remote Sensing, 30(4), 1003–1018.

123