On Wed, 10 Jan 2018, Sima Usvyatsov wrote:
Hello,
I am running a negative binomial model (MASS) on count data collected on a
grid. The dataset is large - ~4,000 points, with many predictors. Being
counts, there are a lot of zeroes. All data are collected on a grid with 20
points, with high spatial autocorrelation.
I would like to filter out the spatial autocorrelation. My question is:
since I have very limited spatial info (only 20 distinct spatial
locations), is it possible to simplify ME() so that I don't have to run it
on the whole dataset? When I try to run ME() on a 100-point subset of the
data, I get error in glm.fit: NA/NaN/Inf in 'x'. When I run it on a single
instance of the grid, I "get away" with a warning ("algorithm did not
converge").
Here's a fake dataset. It was grinding for a while but not throwing errors
(like my original data would). Regardless, it demonstrates the repeated
sampling at the same points and the large number of zeroes.
The data set has 1000 values in Lon, so is probably bigger than you
intended, and when 100 is used is not autocorrelated. You seem to have a
hierarchical model, with repeated measurements at the locations, so a
multi-level treatment of some kind may be sensible. If you want to stay
with ME-based spatial filtering, maybe look at the literature on spatial
panel (repeated measurements are in time) with ME/SF, and on network
autocorrelation (dyadic relationships with autocorrelation among origins
and/or destinations). Both these cases use Kronecker products on the
selected eigenvectors, I think.
Alternatively, use a standard GLMM with a grouped iid random effect and/or
a spatially structured random effect at the 20 location level. If the
groups are repeated observations in time, you should model the whole
(non-)separable space-time process.
Hope this helps,
Roger
Any advice would be most welcome.
library(spdep)
library(MASS)
df <- data.frame(Loc = as.factor(rep(1:20, each = 5)), Lat = rnorm(100, 30,
0.1), Lon = rnorm(1000, -75, 1), x = rnegbin(100, 1, 1))
coordinates(df) <- ~Lon + Lat
proj4string(df) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
nb <- dnearneigh(x=coordinates(df), d1=0, d2=200,longlat = TRUE)
dists <- nbdists(nb, coordinates(df), longlat=TRUE)
glist <- lapply(dists, function(x) 1/x)
lw <- nb2listw(nb, glist, style="W")
me <- ME(x ~ 1, data = df, family = "quasipoisson", listw = lw, alpha = 0.5)
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
