Calculating spatial lags per year in R

Question

At the moment I am having some difficulties with calculating a spatial lag in R. I know how to calculate the lag in space-wide format but am unable to do it in long form, i.e. have repeated observations for the unit of analysis.

Below is some mock data to illustrate what I am trying to do. Let's start by generating some observations of events that I'm interested in.

# Create observations
pts<-cbind(set.seed(2014),x=runif(30,1,5),y=runif(30,1,5),
       time=sample(1:5,30,replace=T))
require(sp)
pts<-SpatialPoints(pts)

x and y are the coordinates while time represents the time period in which the event takes place. The events need to be aggregated to polygons which is the unit of analysis. In this example the polygons are grid-cells and for simplicity the boundaries are fixed over time.

# Observations take place in different areas; create polygons for areas
X<-c(rep(1,5),rep(2,5),rep(3,5),rep(4,5),rep(5,5)) 
Y<-c(rep(seq(1,5,1),5))
df<-data.frame(X,Y)
df$cell<-1:nrow(df) # Grid-cell identifier
require(raster)
coordinates(df)<-~X+Y 
rast<-raster(extent(df),ncol=5,nrow=5)
grid<-rasterize(df,rast,df$cell,FUN=max)
grid<-rasterToPolygons(grid) # Create polygons 

We can plot the data just to get an overview of the distribution:

For space-wide format I would calculate the spatial lag the following way:

pointsincell=over(SpatialPolygons(grid@polygons),SpatialPoints(pts),
              returnList=TRUE)
grid$totalcount<-unlist(lapply(pointsincell,length))
require(spdep)
neigh<-poly2nb(grid) # Create neighbour list
weights<-nb2listw(neigh,style="B",zero.policy=TRUE) # Create weights (binary)
grid$spatial.lag<-lag.listw(weights,
                            grid$totalcount,zero.policy=TRUE) # Add to raster

However, as you can see doing it this way doesn't take into account the fact that the events happen at different moments at time. It just simply aggregates everything to the polygon level. Now I want to calculate this spatial-lag taking into account this temporal dimension so aggregating the data in this case to the polygon-time level.

I wonder whether anyone has a useful suggestion on how this could be accomplished? What is the most convenient way of calculating spatial lags in long format?

I had a look at the spacetime package but was unsuccessful in applying it.

Answer 1

I think the most easy way to accomplish this is to use loops, and create the lag.listw() for your count variable for each year.

Something like this?

spatlag <- data.frame(id=NULL, time=NULL, lag=NULL)
for (y in sort(unique(data$time))){
  print(y)

Then inside the for loop you subset both the points and polygons, and execute the overlay. Then you summarize the number of points for each time point and bind them to the spatlag dataframe, one time point at the time.

pointsincell=over(SpatialPolygons(grid@polygons),SpatialPoints(pts),
              returnList=TRUE)
grid$totalcount<-unlist(lapply(pointsincell,length))
require(spdep)
neigh<-poly2nb(grid) # Create neighbour list
weights<-nb2listw(neigh,style="B",zero.policy=TRUE) # Create weights (binary)
grid$spatial.lag<-lag.listw(weights,grid$totalcount,zero.policy=TRUE) # Add to raster
rbind(spatlag, grid)
}

The code above is just for exemplification. So: 1. Create empty data frame for storing the lags 2. For loop for each time point 3. Create subset for the points where time equals time in for loop run 4. Overlay the points on the grid/polygon 5. Sum the number of points in each polygon overlay (could use dplyr to aggregate) 6. Bind the summed number of points to the empty data frame.

Answer 2

This would be much easier using the slag function of the splm package.

Tell R your data.frame is a panel data frame, then work with the pseries.

Please, note this will work only with balanced panel. Just to give you an example:

library(plm)
library(splm)
library(spdep)

data("EmplUK", package = "plm")

names(EmplUK)
table(EmplUK$year)
#there should be 140 observations /year, but this is not the case, so tomake it balanced

library(dplyr)
balanced_p<-filter(EmplUK, year>1977 & year<1983)
table (balanced_p$year)
#now it is balanced

firm<-unique(balanced_p$firm)
#I'm using the coordinates (randomly generated) of the firms, but it works also if you use the polygons as you did in your question
coords <- cbind(runif(length(firm),-180,+180), runif(length(firm),-90,+90))
pts_firms<-SpatialPoints(coords)

#now tell R that this is a panel, making sure that the firm id and years are the first two columns of the df
p_data<-pdata.frame(balanced_p)
firm_nb<-knn2nb(knearneigh(pts_firms))
firm_nbwghts<-nb2listw(firm_nb, style="W", zero.policy=T)

#now you can easily create your spatial lag 
#I'm assuming here that the dependent variable is wage! 
p_data$splag<-slag(p_data$wage,firm_nbwghts)

p_data$wage is of class pseries, while firm_nbwghts a listw

Answer 3

So I think I've found a method to do this. The output data will come in the form of a normal data frame. It's a bit clumsy but it works.

# Start by creating a panel (CSTS) data frame
grid$cc<-1:nrow(grid)
tiempo<-1:5
polygon<-as.vector(unique(unlist(grid$cc,use.names=FALSE)))

# Loop to create panel data frame
timeCol<-rep(tiempo,length(polygon))
timeCol<-timeCol[order(timeCol)]

polCol <- character()
for(i in tiempo){ 
 row <- polygon
 polCol <- c(polCol, row)
}

df<-data.frame(time=timeCol,nrow=polCol)
df$nrow<-as.numeric(df$nrow)
df<-df[order(df$time,df$nrow),] # Order data frame 

# Assign each point to its corresponding polygon
pts<-SpatialPointsDataFrame(pts,data.frame(pts$time)) # This is a bit clumsy
pts$nrow=over(SpatialPoints(pts),SpatialPolygons(grid@polygons),
              returnlist=TRUE) 

# Aggregate the data per polygon
pts$level<-1
pts.a<-aggregate(level~nrow+time,pts,FUN=sum) # No NA's

# Merge the data
df2<-merge(df,pts.a,all.x=T)
df2[is.na(df2$level),]$level<-0 # Set NA's to 0

# Create contiguity matrix
k<-poly2nb(grid,queen=TRUE) # Create neighbour list
W<-nb2listw(k,style="B",zero.policy=TRUE) # Spatial weights; binary
con<-as.matrix(listw2mat(W)) # Contiguity matrix

# Calculate spatial lag using Kronecker product
N<-length(unique(df2$nrow))
T<-length(unique(df2$time))
ident<-diag(1,nrow=T)
df2$SpatLag<-(ident%x%con)%*%df2$level # Done