Message added to fitted.sarlm() to direct users' attention to use of known response.

See https://github.com/rsbivand/iase/issues/5

The summary object from various spatial models has a Coef component:

> library(spatialreg)
> example(spautolm)
> summary(ecarIVaw)$Coef
             Estimate  Std. Error  z value     Pr(>|z|)
(Intercept) 1.4347053 0.225520449 6.361753 1.994644e-10
ft.NWBIR74  0.0409029 0.006298975 6.493580 8.382028e-11

The summary of an lm or glm object names this component coef with a lower-case "c":

> summary(lm(y~x,data.frame(x=runif(3),y=runif(3))))$coef
              Estimate Std. Error    t value  Pr(>|t|)
(Intercept)  0.5793891  0.4190393  1.3826605 0.3986232
x           -0.3553338  0.5276859 -0.6733813 0.6227157

Is it worth breaking anyone's existing code for a little bit of consistency across models?

consider moving expm from Imports: to Suggests: try to facilitate webr

Hello,
Regarding spatialreg::lmSLX version 1.2.3, I get the following error which I can not figure it out
Error in sum_lm_model$coefficients[((m2 + 2):m + LI), 1:2, drop = FALSE] :
subscript out of bounds
What could be the issue?
Thank you!

Hello there. This is such a great package - thank you for the work in creating it. I am curious about best practices for calculating average marginal effects for interaction terms while using errorsarlm model objects. Standard packages that do this with many other model objects do not work with Sarlm model objects created by the errorsarlm model prediction.

Consider the model Y ~ X1 + X2 + X1*X2, where Y is continuous, X1 is continuous, and X2 is categorical. I calculate this basic OLS interaction model within the errorsarlm command to consider the role of correlated errors in the OLS, and the roll of spillover from the model Xs. I use the impacts command to view the total effects (direct + indirect + lamda), total SE, and total p-value for each coefficient.

I am trying to calculate the beta coefficient and standard error/p-value for the interaction term at various levels of X1. Typically, one uses the margins package or ggeffects package to do this. However, they do not handle errorsarlm model objects. I have inquiries into ggeffects and margins about expanding thier model capabilities, but the issue of considering the spatail weights matrixes when calculating SE I think will prevent these packages from considering a solution.

I have "brute forced" the total effects beta coefficients from the impacts command into the nested lm regression model object in order to use ggeffects and margins, but the standard errors are wrong, of course. The lm model object does not use the weights matrix for accounting for Rho, lamda, etc.

I can also easily calculate point estimates at different levels of the continuous X variable, but I am unsure how to calculate correct SEs. ad different X1 values.

Are there best practices for considering marginal effects at various levels of the interaction model terms? Or is this completed within the impacts command and I am misunderstanding the output or flexibility of impacts? Alternatively, can I use the impacts output and model object to do the work manually? (stated another way: Is the impacts function flexible to predict beta, SE, and p at various levels?)

Seems like others have similar issues, and I have found no real workaround in the Stack communities. See here.

Attached is a short script with a toy example (nonsense model), but the issue exists for all errorsarlm model objects, I believe.

marginal_effects_help.txt

Thank you for you attention and help. Apologies if this question has been considered already, or if I am misunderstanding or underutilizing the impacts output.

https://stat.ethz.ch/pipermail/r-sig-geo/2019-April/027230.html

I am sorry to have to disturb users of the spdep package who have used it to
fit spatial regression models.

The package has become large and hard to maintain. Functions to create and
handle spatial neighbours, and to test for spatial autocorrelation remain in
spdep. Functions to fit spatial regression models are now in the new package
spatialreg (1.1-3 on CRAN). They are also in spdep 1.1-2 now on its way to
CRAN.

If you use spdep as you have done until now, the versions of the functions
in spdep will be used if you have not installed spatialreg, and you will see
warnings that the spdep functions are deprecated. If you have installed
spatialreg, you will see warnings too, but the estimation will be done using
the spatialreg namespace internally. If you install spatialreg, and attach
spdep followed by spatialreg, you will see that the spdep model fitting
functions are masked by their equivalents in spatialreg.

During May, I expect to make the spdep versions of the model fitting
functions defunct, gradually removing the masking problem. Most use of spdep
by packages importing its namespace or attaching it only use functions
creating and handling neighbour objects. The division of functionalities
will benefit that majority of packages, because spdep will load faster and
draw in fewer packages on which it depends. However, improving workflows for
users of those packages means that users of model fitting functions will
need to use spdep for neighbour objects and in addition spatialreg for model
estimation and evaluation.

This split was discussed on github:

r-spatial/spdep#31

and twitter: https://twitter.com/RogerBivand/status/1105023658341351424

Please take up issues on github:

https://github.com/r-spatial/spdep/issues
https://github.com/r-spatial/spatialreg/issues

on the list or by direct email.

I hope that this reconfiguration will improve the software going forward.

Linked spdep issue: r-spatial/spdep#31

Dear Roger,

I just went through your teaching materials of ECS530. Thanks a lot for this valuable resource, it is really incredibly helpful!

However, I was a bit confused by seeing the indirect impacts after running a Durbin error model in ECS530_h19/ECS530_VII.Rmd, line 118. If I'm not confusing something here, the impacts() function after the Durbin error and the SLX seems to take up the wrong lagged coefficients from the model if a lagged intercept is included. See the example (taken from your materials) below:

library(spatialreg)
library(spdep)
library(sf)

nc <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE)
st_crs(nc) <- "+proj=longlat +datum=NAD27"
row.names(nc) <- as.character(nc$FIPSNO)
nc$ft.SID74 <- sqrt(1000)*(sqrt(nc$SID74/nc$BIR74) + sqrt((nc$SID74+1)/nc$BIR74))
nc$ft.NWBIR74 <- sqrt(1000)*(sqrt(nc$NWBIR74/nc$BIR74) + sqrt((nc$NWBIR74+1)/nc$BIR74))
gal_file <- system.file("weights/ncCR85.gal", package="spData")[1]
ncCR85 <- read.gal(gal_file, region.id=nc$FIPSNO)
closeAllConnections()
nc_lw <- nb2listw(ncCR85, style="B")
nc_lw_w <- nb2listw(ncCR85, style="W") # row-standardized to test without intercept

# SDEM with unstandardized W and intercept in lagX
m1b <- spatialreg::errorsarlm(ft.SID74 ~ ft.NWBIR74, weights=BIR74, data=nc, listw=nc_lw, Durbin=TRUE)
summary(m1b)
summary(spatialreg::impacts(m1b))

# The indirect impact does not equal the lagged coef of ft.NWBIR74
all.equal(spatialreg::impacts(m1b)[[1]]$indirect, m1b$coefficients["lag.ft.NWBIR74"], 
          tolerance = 1e-5, check.attributes = FALSE)

# it equals the lagged intercept
all.equal(spatialreg::impacts(m1b)[[1]]$indirect, m1b$coefficients["lag.(Intercept)"], 
          tolerance = 1e-5, check.attributes = FALSE)


# Similarly, with more than one covariate, it starts with the intercept, and omits the last lag coefficient
m1c <- spatialreg::errorsarlm(ft.SID74 ~ ft.NWBIR74  + east, weights=BIR74, data=nc, listw=nc_lw, Durbin=TRUE)
summary(m1c)
summary(spatialreg::impacts(m1c))

# Indirect impact for "east" equals the lagged coefficient for "ft.NWBIR74"
all.equal(spatialreg::impacts(m1c)[[1]]$indirect["east"], m1c$coefficients["lag.ft.NWBIR74"], 
          tolerance = 1e-5, check.attributes = FALSE)


# Without intercept (row standardized W), everything works as expected
m1c <- spatialreg::errorsarlm(ft.SID74 ~ ft.NWBIR74  + east, weights=BIR74, data=nc, listw=nc_lw_w, Durbin=TRUE)
summary(m1c)
summary(spatialreg::impacts(m1c))



### Same problem occurs in impacts(lmSLX()) with lagged intercept
m2 <- spatialreg::lmSLX(ft.SID74 ~ ft.NWBIR74  + east, weights=BIR74, data=nc, listw=nc_lw, Durbin=TRUE)
summary(m2)
summary(spatialreg::impacts(m2))

# but everything works as expected without intercept
m2b <- spatialreg::lmSLX(ft.SID74 ~ ft.NWBIR74  + east, weights=BIR74, data=nc, listw=nc_lw_w, Durbin=TRUE)
summary(m2b)
summary(spatialreg::impacts(m2b))

To me, this looks like the lagged intercept is not taken into account when extracting the lagged coefficients for the impacts in SDEM and SLX (I did not test any other models so far)?

I'm happy to look into the code if you agree that there's something going wrong here.

Best wishes,
Tobias

Hey,
I stumbled upon some arguably unexpected behaviour with the lmSLX() function that exist on version 1.2-1.

1. It throws an error (probably due to a lagged intercept coefficient) if row-stochastic weights are provided naively.
2. If the weights are not row-stochastic a spatially lagged intercept is included by default. This is arguably unexpected behaviour -- what do you think?
3. Also, the impacts of SLX models are not computed in the correct manner (see Issue #23)

I can work on a PR fixing these, but I'll need some time. Below is a quick example.

library("spatialreg")
N <- 100L

# Build connectivity ---
xy <- cbind(runif(N), runif(N)) # Simulate locations
dist <- as.matrix(dist(xy)) # Compute distance
diag(dist) <- Inf # We want zeros on the diagonal
W <- dist^-2 # Compute inverse-distance with decay parameter two
W_rs <- W / rowSums(W) # Use row- and eigenvalue-scaling
W_ev <- W / max(eigen(W, symmetric = TRUE, only.values = TRUE)$values)

# Build DGP ---
X <- matrix(rnorm(N), N)
y <- X * 2 + W_ev %*% X * 1 + rnorm(N, sd = 0.01) # No need for noise

# Estimate using lmSLX() ---
l_rs <- spdep::mat2listw(W_rs)
l_ev <- spdep::mat2listw(W_ev)

lmSLX(y ~ X, listw = l_rs)
# Error in sum_lm_model$coefficients[((m2 + 2):m + LI), 1:2, drop = FALSE] : 
#  subscript out of bounds
lmSLX(y ~ 0 + X, listw = l_rs) # This works
lmSLX(y ~ X, Durbin = ~ 0 + X, listw = l_rs) # This works -- arguably expected
lmSLX(y ~ X, listw = spdep::mat2listw(W_rs, style = "W")) # This works too


lmSLX(y ~ X, listw = l_ev)
# Here, the intercept is lagged too, since W's scaling means it's not just 1
# This may qualifies as unexpected behaviour and can be annoying to avoid
lmSLX(y ~ 0 + X, listw = l_ev) # No intercept
lmSLX(y ~ X, Durbin = ~ 0 + X, listw = l_ev) # Arguably the expected behaviour

Thread https://stat.ethz.ch/pipermail/r-sig-geo/2019-October/027717.html

Dear Roger,

Perhaps you are already aware of this issue, but it seems that the gtools package is orphaned.

Because of this, all packages (like mine) that indirectly depends on spatialreg are now seeing build failures with the warning in R CMD check

Check: package dependencies
Result: WARN
    Requires orphaned package: ‘gtools’

How could I fix this problem?

Best regards,
Aritz

I was reproducing some results to compare spatialreg (version 1.2-1) and bsreg and stumbled upon an issue with the calculation of impacts for the SLX term (i.e. impacting SLX, SDM, and SDEM models).
For the SLX, average indirect impacts are often assumed to be reflected by θ directly. However, this is not generally the case. Partial effects of a univariate SLX are given by

∂y/∂x = Iβ + Wθ

The average indirect effect (i.e. observation i on j != i) is given by sum(W) / N * theta. For row-stochastic W we have sum(W) == N, so the impact is given by θ, but not for other scaling types. I am currently wrapping up a working paper, where I go into some more detail, that I'll link here later.

The issue is that this is not accounted for in the impacts() methods for models with the LX term, as I demonstrate below. I think this can be fixed in a straightforward manner and I can start work on a PR, but will need some time. A code example of what I am talking about is below.

library("spatialreg")
set.seed(42)
N <- 100L

# Build connectivity ---
xy <- cbind(runif(N), runif(N)) # Simulate locations
dist <- as.matrix(dist(xy)) # Compute distance
diag(dist) <- Inf # We want zeros on the diagonal
W <- dist^-2 # Compute inverse-distance with decay parameter two
W_rs <- W / rowSums(W) # Use row-, eigenvalue-, and norm-scaling
W_ev <- W / max(eigen(W, symmetric = TRUE, only.values = TRUE)$values)
W_fn <- W / norm(W, "F")

# Build DGP ---
X <- matrix(rnorm(N), N)
beta <- 2
theta <- 1
y <- X * beta + W_ev %*% X * theta + rnorm(N, sd = 0.01) # No need for noise

# Estimate SLX (using lm and spatialreg) ---

# Use lm and compare theta
(m1 <- lm(y ~ 0 + X + W_ev %*% X)) # True model -- correct
(m2 <- lm(y ~ 0 + X + W_rs %*% X)) # Appears way off
(m3 <- lm(y ~ 0 + X + W_fn %*% X))
# Results mirror spatialreg
(s1 <- lmSLX(y ~ 0 + X, listw = spdep::mat2listw(W_ev)))
(s2 <- lmSLX(y ~ 0 + X, listw = spdep::mat2listw(W_rs)))
(s3 <- lmSLX(y ~ 0 + X, listw = spdep::mat2listw(W_fn)))
# The thetas are very different, depending on W's scale

# Compute average effects
c("Direct" = coef(m1)[[1]], "Indirect" = sum(W_ev) / N * coef(m1)[[2]])
c("Direct" = coef(m2)[[1]], "Indirect" = sum(W_rs) / N * coef(m2)[[2]])
c("Direct" = coef(m3)[[1]], "Indirect" = sum(W_fn) / N * coef(m3)[[2]])
# W_ev and W_fn are equivalent since they only differ by scaling -- W_rs is off
# since it's not symmetric, but effects are closer than theta would imply.

# Compute average effects using spatialreg
impacts(s1)
impacts(s2)
impacts(s3)
# The package just uses the coefficients and ignores the strucutre of W

# Check against true model ---
c("Direct" = beta, "Indirect" = sum(W_ev) / N * theta,
  "Direct (empirical)" = mean(((X + 1) * beta + W_ev %*% X %*% theta) - y),
  "Indirect (empirical)" = mean((X * beta + W_ev %*% (X + 1) %*% theta) - y))
# Partial effects computed manually seem valid

In r-spatial/spdep#160 and r-spatial/spdep#161, changes are made to the use of n.comp.nb. It seems now more important to move checks of the possible existence of disjoint subgraphs to the creation and modification of nb neighbour objects, adding these to the nb object at creation or modification in spdep. This means that separate checking may be carried out on the "ncomp" attribute of an nb object from spdep version >= 1.3-6 (branch n_comp until r-spatial/spdep#161 is merged).

Seen this morning when trying to update spdep, which now entails installing spatialreg :

> install.packages("spatialreg")
Installation du package dans ‘/usr/local/lib/R/site-library’
(car ‘lib’ n'est pas spécifié)
essai de l'URL 'https://cloud.r-project.org/src/contrib/spatialreg_1.1-8.tar.gz'
Content type 'application/x-gzip' length 764534 bytes (746 KB)
==================================================
downloaded 746 KB

* installing *source* package ‘spatialreg’ ...
** package ‘spatialreg’ successfully unpacked and MD5 sums checked
** using staged installation
** libs
gcc -std=gnu99 -std=gnu11 -I"/usr/share/R/include" -DNDEBUG      -fpic  -g -O2 -ffile-prefix-map=/build/r-base-PT7Nxy/r-base-4.1.2=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g  -c eminmaxC.c -o eminmaxC.o
gcc -std=gnu99 -std=gnu11 -I"/usr/share/R/include" -DNDEBUG      -fpic  -g -O2 -ffile-prefix-map=/build/r-base-PT7Nxy/r-base-4.1.2=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g  -c init.c -o init.o
gcc -std=gnu99 -std=gnu11 -I"/usr/share/R/include" -DNDEBUG      -fpic  -g -O2 -ffile-prefix-map=/build/r-base-PT7Nxy/r-base-4.1.2=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g  -c listw2Matrix.c -o listw2Matrix.o
gcc -std=gnu99 -std=gnu11 -I"/usr/share/R/include" -DNDEBUG      -fpic  -g -O2 -ffile-prefix-map=/build/r-base-PT7Nxy/r-base-4.1.2=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g  -c ml_sse.c -o ml_sse.o
gcc -std=gnu99 -std=gnu11 -I"/usr/share/R/include" -DNDEBUG      -fpic  -g -O2 -ffile-prefix-map=/build/r-base-PT7Nxy/r-base-4.1.2=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g  -c mom_calc.c -o mom_calc.o
gcc -std=gnu99 -std=gnu11 -shared -L/usr/lib/R/lib -Wl,-z,relro -o spatialreg.so eminmaxC.o init.o listw2Matrix.o ml_sse.o mom_calc.o -llapack -lblas -lgfortran -lm -lquadmath -L/usr/lib/R/lib -lR
installing to /usr/local/lib/R/site-library/00LOCK-spatialreg/00new/spatialreg/libs
** R
** inst
** byte-compile and prepare package for lazy loading

 *** caught segfault ***
address 0x5568764bd5c1, cause 'memory not mapped'
An irrecoverable exception occurred. R is aborting now ...
Segmentation fault
ERROR: lazy loading failed for package ‘spatialreg’
* removing ‘/usr/local/lib/R/site-library/spatialreg’

The downloaded source packages are in
	‘/tmp/Rtmpil1u9K/downloaded_packages’
Message d'avis :
Dans install.packages("spatialreg") :
  installation of package ‘spatialreg’ had non-zero exit status

As a consequence, spdep is not upgraded...

Thank you for the great spatial tools you provide to the community.

Maybe I'm doing something wrong as I'm rather new to spatial regressions, but here is a behavior I find strange.

if I fit a spatial regression model the way it is suggested in the man page, it works and I get some fitted values in the object I get. For instance :
formula1<-paste0( "to_predict ~ predictor_1 + predictor_2") fit.sem1<-lagsarlm(formula1, data=df_b, listw = neighbs_weights_b)

will work. Then my fit.sem object has a $fitted.values and it seems everything works. df_b is an sf object, and neghbs_weights_b have been produced as intended.

But when I try to use this model to predict my dependent variable on another location, using something like :
predict.sem2.arras<-predict(object=fit.sem2, newdata = df_a, listw = neighbs_weights_a, pred.type = "TS", legacy.mixed=TRUE, power=TRUE)

I get an error message :

Error in x$terms (explo_spat_cor_models.R#210): $ operator is invalid for atomic vectors

It is a bit confusing, as I thought my data was wrong, like missing one of the predictors. I checked the way the weights are provided, including how the row.names need to be specified so that they are different in the training area and the testing area.

I found, thanks to the code of an article by Thibault Laurent (thank him for making it available), that the proper use seems to be:

1. Fit the non-spatial linear model on the "training data".
2. Feed the resulting linear model to lagsarlm in lieu of just the formula
3. Use the resulting spatial model to predict on the new area (providing that all the weights are properly computed and fed to the function, o course).

All in all, as it is working, maybe it is just that I'm not aware of practice that may seem obvious to the community, and if so, sorry for the inconvenience, but as an R user, it is rather puzzling that feeding the formula "almost" works and feeding the linear object "fully" works. Maybe just modifying the error message in predict.sarlm would be useful ?

If the issue is considered useful and the behavior can not be reproduced, I'll try and provide a small data set to reproduce it.

Is there a way to obtain confidence interval for a spautolm object?

While it is straightforward to get confidence intervals for objects from functions like lagsarlm by applying confint(mod) to, it didn’t work for spautolm object.
Could anyone please provide some guidance on how to get confidence interval from a spautolm object?

Instead, I get an error:

  no applicable method for 'vcov' applied to an object of class "Spautolm"```

Hello! The package "spatialreg" is great, and try it this morning and encounter some problems:
how to estimate SAR model with GMM, here I know the function "gstsls" is feasible for SARAR model, how about a SAR model without spatial error term?

Handling of formula Durbin objects does not work, get to line 233 in R/predict.sarlm.R.

We should document the following:

1. The mentioned method= argument to functions now named lagsarlm(), errorsarlm() and sacsarlm(), that is the spdep maximum likelihood estimation functions;

2. The recently updated and modified Durbin= argument to lagsarlm(), errorsarlm() and sacsarlm(), as well as lmSLX(), that is the functions in Halleck Vega and Elhorst (2015) named GNS, SDM, SDEM and SLX, the ones with WX included. In particular, the one-sided formula interface to the arguments permits subsets of X to be lagged. This came up in a discussion in an ERSA workshop in Vienna in 2016, that spatial lags of dummies may not make sense, so I implemented it to give us something to go on in the maximum likelihood setting.

3. Very recently, I've fitted in Bayesian estimation for SEM, SDEM, SLM, SDM, SAC and GNS, using the same Durbin= framework, and the internals of Virgilio Goméz-Rubio and Abhirup Mallik's GSoC project in 2011, translating the Spatial Econometrics toolbox code for sar_g, sdm_g, sem_g and sac_g. It's still preliminary, was done for spReg_lag() earlier, and spReg_sem() and spReg_sac() at the last CRAN release. It's still rough, but does have impacts. SEM/SLM have Griddy Gibbs and Metropolis/Hastings for the spatial coefficients, SAC only has Metropolis/Hastings. None have heteroskedastic variance in the code, because it was hard getting it to work, and in the lag case meant doing much more computation in each iteration. I have bits of a script running against Virgilio's SEMCMC implementations, and INLA's "slm" latent model. If you've read this far, might doing something together for the JSS Bayesian SI make sense?

That's enough work for the five years starting 2019, right?

Originally posted by @rsbivand in r-spatial/spdep#23 (comment)

Hello,

I have been working on an R function script that uses the Getis filter to remove the spatial component of variables, and I am wondering if it could be added to the spatialreg package or it is better suited for other packages such as spdep. I have never contributed to a package, but I would love to document the function and do what is needed for it to be used by other scientists/users.

Thanks in advance for your guidance.
Best regards,

Felipe

See r-spatial/spdep#37

Hi, I have this error in my project:

Error in dimnamesGets(x, value) : 
  invalid dimnames given for “dgRMatrix” object

I used some of the code from Spatial weights objects as sparse matrices and graphs first to get data and neighbor matrix:

library(spatialreg)
library(sf)
columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1])
library(spdep)
nb_q <- poly2nb(columbus)
nb_q

attr(nb_q, "region.id")[1:10]
col2 <- droplinks(nb_q, 21)
View(col2)

nb_B <- nb2listw(col2, style="B", zero.policy=TRUE)

After this, it will be my code for reproducing the error.

First, I want to make sure it works with original data, so I extract some of them and run the lag model, it works.

columbus_extract <- columbus[,7:18] %>% st_set_geometry(NULL)

test <- lagsarlm(HOVAL~., data=columbus_extract, listw = nb_B, zero.policy = TRUE)
summary(test)

Next, I did a similar processing to what I have done in my project.
I separated the original data by 2 time variables, which are season and time periods of one day(like 0-3am, 3-6am...) and make the the number of lines has become 32(4 seasons * 8 periods) times of the original.
Since there is no time-related data in columbus, I just simply copied the columbus_extract data 32 times, just for keeping the same structure.

columbus_new <- data.frame(matrix(data=NA,nrow = 1568,ncol = 12))

for(i in 1:12){
  columbus_new[,i] <- rep(columbus_extract[,i],32)
}
colnames(columbus_new) <- colnames(columbus_extract)
View(columbus_new)

For the neighbor matrix, I use this code to copy. Now the neighbor matrix is like a combined 32 separated neighbor matrix of the original one.

col2_new <- col2
for(k in 2:32) {
  for(l in 1:49) {
    a <- col2_new[l]
    b <- lapply(a,function(x) x + (k-1L)*49)
    col2_new[[l+(k-1)*49]] <- b[[1]]
  }
}
View(col2_new)

for (i in 1:length(col2_new)){
  col2_new[[i]] <- as.integer(col2_new[[i]])
}

nb_B_new <- nb2listw(col2_new, style="W", zero.policy=TRUE)
View(nb_B_new)

Then when I run the lag model, whether I have added method="LU" or not, I got an error. And in the errorsarlm model I got the same.

test2 <- lagsarlm(HOVAL~., data=columbus_new, listw = nb_B_new, zero.policy = TRUE)

test3 <- lagsarlm(HOVAL~., data=columbus_new, listw = nb_B_new, zero.policy = TRUE,
                  method = "LU")

Error:

Error in dimnamesGets(x, value) : 
  invalid dimnames given for “dgRMatrix” object

But in the SLX model, it's working fine.

test4 <- lmSLX(HOVAL~., data=columbus_new, listw = nb_B_new, zero.policy = TRUE)
summary(test4)

In another case, when I try to assign the new neighbor matrix to "CsparseMatrix" as you have shown in Spatial weights objects as sparse matrices and graphs, I got the same error.

# another case
B <- as(nb_B_new, "CsparseMatrix")

In my case, is there any way to solve this error?
And will you recommend using spatial regression in this way？
Thanks!

When observation weights are used the reported AIC is the standard AIC. When observation weights are used a dAIC measure might be more appropriate, see http://www.isr.umich.edu/src/smp/asda/J%20Surv%20Stat%20Methodol-2015-Lumley-jssam_smu021.pdf. Moved from r-spatial/spdep#13

@edzer how may I set up Travis on this repo, as it is part of r-spatial, and r-spatial already has Travis authorization?

spData is migrating from ESRI Shapefile to GPKG, so use of ESRI Shapefiles from spData in this package also needs to migrate before or by the end of August 2024, when ESRI Shapefiles will be dropped from spData. Needs 2.3.2 because of missing NY8_utm18.gpkg in 2.3.1; the spData declaration in DESCRIPTION should be changed to spData (>= 2.3.2).

I have some data here available at: https://github.com/zhenchun/data

require(spdep)
require(spatialreg)
require(sf)

nc = st_read("nc.shp", package="sf")

nb <- poly2nb(nc, queen=TRUE)

wt<-nb2listw(nb,zero.policy=TRUE, style="W")

fit<-glm.nb(n~pnhw13_17+offset(log(tt13_17)), nc, na.action=na.omit)
> summary(fit)

#Call:
#glm.nb(formula = n ~ pnhw13_17 + offset(log(tt13_17)), data = nc, 
#    na.action = na.omit, init.theta = 0.6835949402, link = log)

#Deviance Residuals: 
#    Min       1Q   Median       3Q      Max  
#-2.4175  -1.1024  -0.5685   0.1666   4.1108  

#Coefficients:
 #           Estimate Std. Error z value Pr(>|z|)    
#(Intercept)  -8.2594     0.1569  -52.63   <2e-16 ***
#pnhw13_17     3.0596     0.2121   14.43   <2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#(Dispersion parameter for Negative Binomial(0.6836) family taken to be 1)

#    Null deviance: 1132.94  on 800  degrees of freedom
#Residual deviance:  910.99  on 799  degrees of freedom
#AIC: 5665.7

#Number of Fisher Scoring iterations: 1


#              Theta:  0.6836 
#          Std. Err.:  0.0390 

# 2 x log-likelihood:  -5659.6900 

me.fit<-spatialreg::ME(n~pnhw13_17,offset=offset(log(tt13_17)), data = nc, family=negative.binomial(0.68236, link=log),zero.policy=TRUE, listw = wt, verbose=T,alpha=.05,  na.action=na.omit)


#Error in glm.fit(x = iX, y = Y, weights = weights, offset = offset, family = family) : 
 # NA/NaN/Inf in 'x'

Anybody knows how this error come from? Thanks a lot for your help.

I also learnt from these links https://www.mail-archive.com/[email protected]/msg03967.html
https://rpubs.com/corey_sparks/111362

Thanks for making this package.

When porting over sphet::stslshac, please consider either enabling a user to run OLS regressions with SHAC standard errors or provide a separate function; i.e. have both shac and stslshac.

This is probably an issue more related to the documentation.
I asked this question on stack overflow but I think it is helpful to bring it here.

I am fitting a Spatial Error Model using the errorsarlm() function in the spdep library.
The Breusch-Pagan test for spatial models, calculated using the bptest.sarlm() function, suggest the presence of heteroskedasticity.

A natural next step would be to get the robust standard error estimates and update the p-values. In the documentation of the bptest.sarlm() function says the following:

"It is also technically possible to make heteroskedasticity corrections to standard error estimates by using the “lm.target” component of sarlm objects - using functions in the lmtest and sandwich packages."

and the following code (as reference) is presented:

lm.target <- lm(error.col$tary ~ error.col$tarX - 1)
if (require(lmtest) && require(sandwich)) {
  print(coeftest(lm.target, vcov=vcovHC(lm.target, type="HC0"), df=Inf))} 

where error.col is the spatial error model estimated.

Now, I can easily adapt the code to my problem and get the robust standard errors.
Nevertheless, I was wondering:

• What exactly is the “lm.target” component of sarlm objects? I can not find any mention to it in the spdep documentation.
• What exactly are $tary and $tarX? Again, it does not seem to be mentioned on the documentation.
• Why documentation says it is "technically possible to make heteroskedasticity corrections"? Does it mean that proposed approach is not really recommended to overcome issues of heteroskedasticity?

See #10:

data(oldcol, package="spdep")
listw <- spdep::nb2listw(COL.nb, style="W")
ev <- eigenw(listw)
anyNA(COL.OLD$INC)
COL.OLD1 <- COL.OLD
is.na(COL.OLD1$INC) <- 10L
anyNA(COL.OLD1$INC)
res <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD1, listw, control=list(pre_eig=ev))
# Error in eigen_pre_setup(env, pre_eig = pre_eig, which = which) : 
#  length(pre_eig) == get("n", envir = env) is not TRUE

Hello,

I am using the newest version of spatialreg package 1.3.2 in R Console as well as in RStudio, in both cases getting the Error: 'predict.Sarlm' is not an exported object from 'namespace:spatialreg'.

I cannot find the function neither in listing , which I paste bellow, nor in NAMESPACE file of the package. Is it supposed to be public function? How can I use the estimated errorsarlm model to make predictions?

Thank you and appologize, if the bug is on my side.

ls("package:spatialreg")
[1] "aple" "aple.mc" "aple.plot"
[4] "as.spam.listw" "as_dgRMatrix_listw" "as_dsCMatrix_I"
[7] "as_dsCMatrix_IrW" "as_dsTMatrix_listw" "bptest.Sarlm"
[10] "can.be.simmed" "coerce" "create_WX"
[13] "do_ldet" "eigen_pre_setup" "eigen_setup"
[16] "eigenw" "errorsarlm" "get.ClusterOption"
[19] "get.coresOption" "get.mcOption" "get.VerboseOption"
[22] "get.ZeroPolicyOption" "GMargminImage" "GMerrorsar"
[25] "griffith_sone" "gstsls" "Hausman.test"
[28] "cheb_setup" "impacts" "intImpacts"
[31] "invIrM" "invIrW" "Jacobian_W"
[34] "jacobianSetup" "l_max" "lagmess"
[37] "lagsarlm" "lextrB" "lextrS"
[40] "lextrW" "listw2U_Matrix" "listw2U_spam"
[43] "lmSLX" "localAple" "LR.Sarlm"
[46] "LR1.Lagmess" "LR1.Sarlm" "LR1.Spautolm"
[49] "LU_prepermutate_setup" "LU_setup" "Matrix_J_setup"
[52] "Matrix_setup" "mcdet_setup" "MCMCsamp"
[55] "ME" "mom_calc" "mom_calc_int2"
[58] "moments_setup" "powerWeights" "sacsarlm"
[61] "SE_classic_setup" "SE_interp_setup" "SE_whichMin_setup"
[64] "set.ClusterOption" "set.coresOption" "set.mcOption"
[67] "set.VerboseOption" "set.ZeroPolicyOption" "similar.listw"
[70] "spam_setup" "spam_update_setup" "SpatialFiltering"
[73] "spautolm" "spBreg_err" "spBreg_lag"
[76] "spBreg_sac" "stsls" "subgraph_eigenw"
[79] "trW" "Wald1.Sarlm"

As explained by @rsbivand here, adding a Durbin= argument to splm fitting functions would be helpful for estimating the impacts from panel spatial Durbin models (SDM) as well as drawing inference for those impacts.

This would be possible if spatialreg exposes a single portable function handling such an argument.

Hi,

I have done a spatial regression using spatialreg::lagsarlm and I would like to predict using spatialreg::predict.sarlm however it shows an error that the function is not found.

I have previously used this function on another computer but I cannot find the function anymore.

Has the function been moved? Would you know how can I run predict.sarlm?

Thanks in advance,
Daniel

See this screenshot,

Current predict methods do not look for case weights.

Hi
I am looking into a way to validate a spatial error model I built.

From one hand I would like to perform a cross validation but I can't predict on a dataset with different rows id using the predict function .. and I guess it is because of the neighborhood matrix.

Is there a way to validate the model on a different dataset?

Also is the Rsquared computed as a good measure for the model fit in case of a spatial error ?
Or do you suggest using other measures

Thank you
Angela

Hi,

I am working on some spatial regression models using the errorsarlm() function. Unfortunatly it seems that this function uses all cores on a computer even when it is previously limited to a subset using the parallel package. I was hoping there is a way to either set the number of cores in the funciton itself or reoslve this.

Many thanks!

Justin

library(spatialreg)
library(spdep)
data(oldcol, package = "spdep")
W <- spdep::nb2listw(COL.nb, style = "W")
eq <- CRIME ~ 0 + INC + HOVAL
sdem <- errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed", method = "eigen")

in https://stat.ethz.ch/pipermail/r-sig-geo/2019-August/027557.html.

> sdem <- errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed",
+                    method = "eigen")
Error in estimable.default(lm.target, cm) :
   `cm' argument must be of type vector, list, or matrix.
In addition: Warning messages:
1: Function errorsarlm moved to the spatialreg package
2: In spatialreg::errorsarlm(formula = formula, data = data, listw =
listw,  :
   model configuration issue: no total impacts
> traceback()
6: stop("`cm' argument must be of type vector, list, or matrix.")
5: estimable.default(lm.target, cm)
4: estimable(lm.target, cm)
3: as.matrix(estimable(lm.target, cm)[, 1:2, drop = FALSE])
2: spatialreg::errorsarlm(formula = formula, data = data, listw = listw,
        na.action = na.action, Durbin = Durbin, etype = etype, method =
method,
        quiet = quiet, zero.policy = zero.policy, interval = interval,
        tol.solve = tol.solve, trs = trs, control = control)
1: errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed", method =
"eigen")
> dem <- lmSLX(eq, data = COL.OLD, listw = W)
Error in lmSLX(eq, data = COL.OLD, listw = W) : 
  object 'dirImps' not found
In addition: Warning message:
In lmSLX(eq, data = COL.OLD, listw = W) :
  model configuration issue: no total impacts
> traceback()
1: lmSLX(eq, data = COL.OLD, listw = W)

draft tests for SDM

If pre_eig= is used, and subsetting is done because of incomplete cases, protect against invalid pre-computed values (force re-computing).