Calculating spatial lags per year in R

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.