{"id":544,"date":"2014-02-28T13:17:55","date_gmt":"2014-02-28T18:17:55","guid":{"rendered":"http:\/\/homepages.uc.edu\/~yaozo\/wordpress\/?p=544"},"modified":"2014-02-28T13:17:55","modified_gmt":"2014-02-28T18:17:55","slug":"simple-spatial-correlograms-for-cross-country-analysis-in-r","status":"publish","type":"post","link":"https:\/\/zhuoyao.net\/index.php\/2014\/02\/28\/simple-spatial-correlograms-for-cross-country-analysis-in-r\/","title":{"rendered":"Simple Spatial Correlograms for Cross-Country Analysis in R"},"content":{"rendered":"

Accounting for temporal dependence in econometric analysis is important, as the presence of temporal dependence violates the assumption that observations are independent units. Historically, much less attention has been paid to correcting for spatial dependence, which, if present, also violates this independence assumption.<\/p>\n

The comparability of temporal and spatial dependence is useful for illustrating why spatial dependence may matter. For example, a time-series analysis of a county\u2019s GDP would examine how related this year\u2019s figure was to preceding years, whereas a spatial analysis would examine how related one country\u2019s figure was to that of their neighbours.<\/p>\n

However, these temporal and spatial issues are not completely analogous. This difference is highlighted with the creation of lagged variables. In time-series, the lag structure is intuitively easy to comprehend. A variable y_{t} has lags of y_{t-1}, \u2026, y_{t-k}.<\/p>\n

Lagging a spatial variable introduces more ambiguity. A spatially lagged variable of y_{i} can be calculated as \\rho(y_{i}) where \\rho is a spatial weighting matrix. Interested readers should consult further sources on this, but needless to say there are a number of ways (nearest neighbour, distance etc.) in which we could define \\rho.<\/p>\n

The difference between spatial and temporal processes in the creation of lagged variables is somewhat compounded when we consider that it is not possible to create \u2018lags of lags\u2019 in the same manner as in time-series (so \\rho^2 is not the valid). Your\u00a0neighbour\u2019s\u00a0neighbour might be your neighbour.<\/p>\n

Given the difficulty in creating spatial lags of spatial lags, is there an easy way of recreating a simple time-series correlogram for a spatial variable? Yes, indeed it is. The solution I present here is not an exclusive way to analyze a spatially dependent variable.<\/p>\n

For this analysis I use different countries spatial coordinates and create a correlogram of population density in 1900 (these data are taken from thispaper<\/a>). The first step involves loading these data into the R workspace such that the data frame contains the latitude and longitude coordinates alongside the potentially spatially autocorrelated variable.<\/p>\n

 <\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n<\/td>\n
\n
\n
> pop <- <\/code>read.csv<\/code>(<\/code>\"dat2.csv\"<\/code>)<\/code><\/div>\n
> <\/code>head<\/code>(pop)<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>lon\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 lat\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 country\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 dens<\/code><\/div>\n
1\u00a0 66.02554\u00a0 33.83319\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 Afghanistan\u00a0 79.2543748<\/code><\/div>\n
2\u00a0 17.55091 -12.29862\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 Angola\u00a0 24.2375278<\/code><\/div>\n
3\u00a0 20.07006\u00a0 41.14267\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 Albania 282.7854408<\/code><\/div>\n
4\u00a0 54.33089\u00a0 23.91317 United Arab Emirates\u00a0\u00a0 5.0438742<\/code><\/div>\n
5\u00a0 44.93819\u00a0 40.29368\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 Armenia\u00a0\u00a0 0.2136146<\/code><\/div>\n
6 134.48679 -25.73251\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 Australia\u00a0\u00a0 4.8896695<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

 <\/p>\n

With the data frame, I calculated a distance matrix with the \u2018rdist.earth\u2019 function from the package \u2018fields<\/a>\u2019. This function measures the distance between each set of data points based on their coordinates. This function recognizes that the world is not flat, and as such calculates what are known as great-circle distances.<\/p>\n

 <\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n<\/td>\n
\n
\n
# calculate great cirlce distances with fields package<\/code><\/div>\n
library<\/code>(fields)<\/code><\/div>\n
# popdists matrix<\/code><\/div>\n
popdists <- <\/code>as.matrix<\/code>(<\/code>rdist.earth<\/code>(<\/code>cbind<\/code>(pop$lon, pop$lat), miles = F, R = <\/code>NULL<\/code>))<\/code><\/div>\n
diag<\/code>(popdists) <- 0<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

 <\/p>\n

With this distance based matrix it is possible to calculate a correlogram. This is performed through the function \u2018autocorr\u2019. There are three arguments. The first is w is the distance-based weights matrix. The second, x is the variable of interest (note: the order of x and w should be the same as in the original data frame), and the third are the distance bands.<\/p>\n

 <\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n
13<\/div>\n
14<\/div>\n
15<\/div>\n
16<\/div>\n
17<\/div>\n
18<\/div>\n
19<\/div>\n<\/td>\n
\n
\n
# autocorr function<\/code><\/div>\n
autocorr <- <\/code>function<\/code>(w,x,dist=500){<\/code><\/div>\n
\u00a0\u00a0<\/code>aa <- <\/code>ceiling<\/code>(<\/code>max<\/code>(w)\/dist)<\/code><\/div>\n
\u00a0\u00a0<\/code>dists <- <\/code>seq<\/code>(0,aa*dist,dist)<\/code><\/div>\n
\u00a0\u00a0<\/code>cors <- <\/code>NULL<\/code><\/div>\n
\u00a0\u00a0<\/code>for<\/code>(i <\/code>in<\/code> 1:aa){<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0<\/code>w1 <- <\/code>ifelse<\/code>(w > dists[i] & w <= dists[i+1], 1, 0) <\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0<\/code>w2 <- w1<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0<\/code>for<\/code>(j <\/code>in<\/code> 1:<\/code>dim<\/code>(w1)[1]){<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>nu <- <\/code>sum<\/code>(w1[j,])<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>if<\/code>(nu>0){<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>w2[j,] <- w1[j,]\/nu<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>}\u00a0 <\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0<\/code>}<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0<\/code>lag <- w2 %*% x<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0<\/code>cors <- <\/code>c<\/code>(cors,<\/code>cor<\/code>(x,lag))<\/code><\/div>\n
\u00a0\u00a0<\/code>}<\/code><\/div>\n
\u00a0\u00a0<\/code>return<\/code>(cors)<\/code><\/div>\n
}<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

 <\/p>\n

The function works in the following manner. For each step there is a distance band. Here, I have specified that the distance bands are 500km. For the first-step, the function calculates the correlation of the variable of interest (x) and its spatially lagged neighbours \u2013 who are defined as having a centre-point within a 500km radius. The next step calculates this correlation for neighbours falling into the 500km to 1,000 km distance band, so there will be a new spatial weighting matrix \\rho. This is an iterative procedure, and continues up until the bands exceed the furthest neighbours in the data frame.<\/p>\n

 <\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n
13<\/div>\n
14<\/div>\n
15<\/div>\n
16<\/div>\n
17<\/div>\n
18<\/div>\n
19<\/div>\n
20<\/div>\n
21<\/div>\n
22<\/div>\n
23<\/div>\n
24<\/div>\n
25<\/div>\n
26<\/div>\n<\/td>\n
\n
\n
# at 500km<\/code><\/div>\n
ac1 <- <\/code>autocorr<\/code>(w=popdists,x=pop$dens,dist=500)<\/code><\/div>\n
# 1000 km<\/code><\/div>\n
ac2 <- <\/code>autocorr<\/code>(w=popdists,x=pop$dens,dist=1000)<\/code><\/div>\n
<\/div>\n
# MC analysis<\/code><\/div>\n
it <- 1000<\/code><\/div>\n
mc <- <\/code>matrix<\/code>(<\/code>NA<\/code>,nrow=it,ncol=<\/code>length<\/code>(ac1))<\/code><\/div>\n
for<\/code>(i <\/code>in<\/code> 1:it){<\/code><\/div>\n
\u00a0\u00a0<\/code>pop$rand <- <\/code>sample<\/code>(pop$dens,<\/code>length<\/code>(pop$dens),replace=F)<\/code><\/div>\n
\u00a0\u00a0<\/code>mc[i,] <- <\/code>autocorr<\/code>(w=popdists,x=pop$rand,dist=500)<\/code><\/div>\n
}<\/code><\/div>\n
<\/div>\n
library<\/code>(ggplot2)<\/code><\/div>\n
ac1 <- <\/code>data.frame<\/code>(<\/code>cbind<\/code>(ac1,<\/code>seq<\/code>(500,20000,500)))<\/code><\/div>\n
ac1 <- <\/code>cbind<\/code>(ac1,<\/code>t<\/code>(<\/code>apply<\/code>(mc,2,quantile, probs = <\/code>c<\/code>(0.025,0.975))))<\/code><\/div>\n
names<\/code>(ac1) <- <\/code>c<\/code>(<\/code>\"ac\"<\/code>,<\/code>\"dist\"<\/code>,<\/code>\"lci\"<\/code>,<\/code>\"uci\"<\/code>)<\/code><\/div>\n
<\/div>\n
ggplot<\/code>(ac1, <\/code>aes<\/code>(dist, ac)) +<\/code><\/div>\n
\u00a0\u00a0<\/code>geom_point<\/code>(colour = <\/code>\"darkblue\"<\/code>, size = 3) +<\/code><\/div>\n
\u00a0\u00a0<\/code>geom_line<\/code>(colour = <\/code>\"red\"<\/code>) +<\/code><\/div>\n
\u00a0\u00a0<\/code>scale_x_continuous<\/code>(<\/code>'Great Circle Distance (KM)'<\/code>,limits=<\/code>c<\/code>(0,20000)) + <\/code><\/div>\n
\u00a0\u00a0<\/code>scale_y_continuous<\/code>(<\/code>'Autocorrelation'<\/code>,limits=<\/code>c<\/code>(-0.4,0.7)) +<\/code><\/div>\n
\u00a0\u00a0<\/code>theme_bw<\/code>() + <\/code><\/div>\n
\u00a0\u00a0<\/code>geom_hline<\/code>(yintercept=0) +\u00a0\u00a0 <\/code><\/div>\n
\u00a0\u00a0<\/code>geom_smooth<\/code>(<\/code>aes<\/code>(ymin = lci, ymax = uci), stat=<\/code>\"identity\"<\/code>,fill=<\/code>\"blue\"<\/code>,colour=<\/code>\"darkblue\"<\/code>) <\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

 <\/p>\n

I also calculated 95% confidence intervals using a simple Monte Carlo technique. This code, and that to generate the nice ggplot2 figure is shown below.<\/p>\n

The analysis in the above provides a simple example of how R can be used to generate descriptive spatial statistics. Those interested in implementing spatial methodologies should examine the wide variety of packages on CRAN which can be used for spatial analysis. The \u2018spdep<\/a>\u2018 package provides most of the essentials. The recent release of the \u2018splm<\/a>\u2018\u00a0package, containing functions for panel regression analysis, is also noteworthy, and potentially very important, contribution.<\/p>\n

\"\"<\/a><\/p>\n

 <\/p>\n","protected":false},"excerpt":{"rendered":"

Accounting for temporal dependence in econometric analysis is important, as the presence of temporal dependence violates the assumption that observations are independent units. Historically, much… <\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[17],"tags":[],"class_list":["post-544","post","type-post","status-publish","format-standard","hentry","category-phddissertation"],"_links":{"self":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/544","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/comments?post=544"}],"version-history":[{"count":0,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/544\/revisions"}],"wp:attachment":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/media?parent=544"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/categories?post=544"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/tags?post=544"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}