SLX Model - Spatial Econometrics with panel in R data using splm package and slag function

Question

I need to estimate spatial econometric models with spatial lags of X (SLX) either alone, combined with spatial autoregressive models (SAR) or with spatial error model (SEM). When they are combined, they are called Spatial Durbin Model (SDM) or Spatial Durbin Error Model (SDEM), following Vega & Elhorst's (2015) paper "The SLX Model".

I intend to estimate all spatial panel models in R using splm package, which also requires spdep functions. In this sense, I created neighbors lists type Queen and k = 4 from a shape file:

> TCAL <- readOGR(dsn = ".", "Municipios_csv")
> coords <- coordinates(TCAL)

> contnbQueen <- poly2nb(TCAL, queen = TRUE)
> enter code herecontnbk4    <- knn2nb(knearneigh(coords, k = 4, RANN = FALSE))

Then I converted this neighbor list in a Weight matrix:

> W <- nb2listw(contnbk4, glist = NULL, style = "W")

attributes(W)
$names
[1] "style"      "neighbours" "weights"   

$class
[1] "listw" "nb"   

$region.id
 [1] "1"   "2"   "3"   "4"   "5"   "6"   "7"   "8"   "9"   "10"  "11"  "12"  "13"  "14"  "15"  "16"  "17"  "18"  "19"  "20" 
 [21] "21"  "22"  "23"  "24"  "25"  "26"  "27"  "28"  "29"  "30"  "31"  "32"  "33"  "34"  "35"  "36"  "37"  "38"  "39"  "40" 
 [41] "41"  "42"  "43"  "44"  "45"  "46"  "47"  "48"  "49"  "50"  "51"  "52"  "53"  "54"  "55"  "56"  "57"  "58"  "59"  "60" 
 [61] "61"  "62"  "63"  "64"  "65"  "66"  "67"  "68"  "69"  "70"  "71"  "72"  "73"  "74"  "75"  "76"  "77"  "78"  "79"  "80" 
 [81] "81"  "82"  "83"  "84"  "85"  "86"  "87"  "88"  "89"  "90"  "91"  "92"  "93"  "94"  "95"  "96"  "97"  "98"  "99"  "100"
[101] "101" "102" "103" "104" "105" "106" "107" "108" "109" "110" "111" "112" "113" "114" "115" "116" "117" "118" "119" "120"
[121] "121" "122" "123" "124" "125" "126" "127" "128" "129" "130" "131" "132" "133" "134" "135" "136" "137" "138" "139" "140"
[141] "141" "142" "143" "144" "145" "146" "147" "148" "149" "150" "151" "152" "153" "154" "155" "156" "157" "158" "159" "160"
[161] "161" "162" "163" "164" "165" "166" "167" "168" "169" "170" "171" "172" "173" "174" "175" "176" "177" "178" "179" "180"
[181] "181" "182" "183" "184" "185" "186" "187" "188" "189" "190" "191" "192" "193" "194" "195" "196" "197" "198" "199" "200"
[201] "201" "202" "203" "204" "205" "206" "207" "208" "209" "210" "211" "212" "213" "214" "215" "216" "217" "218" "219" "220"
[221] "221" "222" "223" "224" "225" "226" "227" "228" "229" "230" "231" "232" "233" "234" "235" "236" "237" "238" "239" "240"
[241] "241" "242" "243" "244" "245" "246" "247" "248" "249" "250" "251" "252" "253" "254" "255" "256" "257" "258" "259" "260"
[261] "261" "262" "263" "264" "265" "266" "267" "268" "269" "270" "271" "272" "273" "274" "275" "276"

    $call
    nb2listw(neighbours = contnbk7, glist = NULL, style = "W")

Next step, I created a formula for panel SAR and SEM models, which worked well and produced estimates:

> fmPanel <- Area ~ Dist + Land + CredAg
> vegSAR <- spml(fmPanel, data = veg, index = c("Mun","Year"), listw = W, model = "within", effect = "twoways", spatial.error = "none", lag = TRUE)
> vegSEM <- spml(fmPanel, data = veg, index = c("Mun","Year"), listw = W, model = "within", effect = "twoways", spatial.error = "b", lag = FALSE)

Then, I tried to estimate SLX, SDM and SDEM models by creating spatial lags of the covariates X:

> vegX <- pdata.frame(veg, index = c("Mun","Year")); class(vegX)

    [1] "pdata.frame" "data.frame"

I then created pseries values:

> DistX <- vegX$Dist; class(DistX)
    [1] "pseries" "numeric"
> LandX <- vegX$Land; class(LandX)
    [1] "pseries" "numeric"
> CredAgX <- vegX$CredAg; class(CredAgX)
    [1] "pseries" "numeric"

But the error happened when I applied slag function:

DistX <- slag(agSPX$Dist, listw = W)

    Error in lag.listw(listw, xt) : object lengths differ

My panel data has 5 years and 276 regions. So, objects characteristics are:

> length(DistX)
[1] 1380

> length(W)
[1] 3

> length(W$weights)
[1] 276

So, I was wondering that, if I could transform W$weights in a matrix such as usaww used as example of slag function, I could apply function mat2listw and then use slag over X.

Could someone tell me where I am wrong?

Answer 1

You can also add the lagged explanatory variable to the data.

A reproducible example:

library(plm)
library(spatialreg)
library(splm)

# load data
data(Produc, package = "plm")
data(usaww, package = "splm")

d <- pdata.frame(Produc, index = c("state","year"), drop.index = FALSE)

# create a spatial explanatory variable
d$unemp_l <- slag(d$unemp, usaww)

# run model
m <- splm::spml(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp + unemp_l,
                  data = d, listw = mat2listw(usaww) , model="within")
summary(m)
Spatial panel fixed effects error model


Call:
splm::spml(formula = log(gsp) ~ log(pcap) + log(pc) + log(emp) + 
    unemp + unemp_l, data = d, listw = mat2listw(usaww), model = "within")

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.1211492 -0.0234013 -0.0040218  0.0167919  0.1787587 

Spatial error parameter:
    Estimate Std. Error t-value  Pr(>|t|)    
rho 0.542254   0.033772  16.056 < 2.2e-16 ***

Coefficients:
            Estimate Std. Error t-value Pr(>|t|)    
log(pcap)  0.0090575  0.0251036  0.3608  0.71824    
log(pc)    0.2152367  0.0234077  9.1951  < 2e-16 ***
log(emp)   0.7833003  0.0277672 28.2096  < 2e-16 ***
unemp     -0.0014795  0.0011443 -1.2930  0.19603    
unemp_l   -0.0031210  0.0015790 -1.9766  0.04808 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

You can verify that the spatial explanatory is correct. For instance:

a <- usaww["ALABAMA",]
a <- a[a!=0]
a
    FLORIDA     GEORGIA MISSISSIPPI    TENNESSE 
       0.25        0.25        0.25        0.25 
mean(d[d$year=="1970" & d$state %in% names(a) , "unemp"])
[1] 4.525

d[d$state=="ALABAMA" & d$year=="1970", "unemp_l"]
ALABAMA-1970 
   4.525 

Answer 2

Probably it's not the best solution, but I calculated SLX, SDM and SDEM models with these steps:

1) Loading the Shape file by:

TCAL <- readOGR(dsn = ".", "Municipios_csv_BIO")
coords <- coordinates(TCAL)

2) Creating weighting matrix W through:

contnbQueen  <- poly2nb(TCAL, queen = TRUE)
contnbk4     <- knn2nb(knearneigh(coords, k = 4,  RANN = FALSE))

3) Choosing which matrix to apply:

W <- nb2listw(contnbk4, glist = NULL, style = "W")

4) Converting data.frame in a pdata.frame:

vegSPX  <- pdata.frame(vegPainel, index = c("ID","Ano"))

5) Creating specific pdata.frame spatial vectors, for instance:

vegIDDX  <- vegSPX$IDD
vegSoilX <- vegSPX$Soil
vegQAIX  <- vegSPX$QAI

6) Specifying a formula for these models:

fmSPvegX <- Area ~ IDD + Soil + QAI + slag(vegIDDX, listw = W) + slag(vegSoilX, listw = W) + slag(vegQAIX, listw = W)

7) Applying plm and splm functions on the original data (data.table data.frame object):

vegSLX<-plm(fmSPvegX,data=vegPainel,listw=W,index=c("ID","Ano"),model="within",effect="twoways",spatial.error="none",lag=F)
vegSARX<-spml(fmSPvegX,data=vegPainel,listw=W,index=c("ID","Ano"),model="within",effect="twoways",spatial.error="none",lag=T)
vegSEMX<-spml(fmSPvegX,data=vegPainel,listw=W,index=c("ID","Ano"),model="within",effect="twoways",spatial.error="b",lag=F)

I hope it may help.