SLX 模型 - 使用 splm 包和 slag 函数在 R 数据中使用面板的空间计量经济学
SLX Model - Spatial Econometrics with panel in R data using splm package and slag function
我需要单独估计空间滞后 X (SLX) 的空间计量模型，结合空间自回归模型 (SAR) 或空间误差模型 (SEM)。当它们组合在一起时，它们被称为空间 Durbin 模型 (SDM) 或空间 Durbin 误差模型 (SDEM)，遵循 Vega & Elhorst (2015) 的论文 "The SLX Model".
我打算使用 splm 包估计 R 中的所有空间面板模型，这也需要 spdep 函数。从这个意义上说，我从形状文件创建了 Queen 类型和 k = 4 的邻居列表：
> TCAL <- readOGR(dsn = ".", "Municipios_csv")
> coords <- coordinates(TCAL)
> contnbQueen <- poly2nb(TCAL, queen = TRUE)
> enter code herecontnbk4 <- knn2nb(knearneigh(coords, k = 4, RANN = FALSE))
然后我将这个邻居列表转换成权重矩阵：
> 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")
下一步，我为面板 SAR 和 SEM 模型创建了一个公式，该公式运行良好并生成了估计值：
> 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)
然后，我尝试通过创建协变量 X 的空间滞后来估计 SLX、SDM 和 SDEM 模型：
> vegX <- pdata.frame(veg, index = c("Mun","Year")); class(vegX)
[1] "pdata.frame" "data.frame"
然后我创建了 pseries 值：
> DistX <- vegX$Dist; class(DistX)
[1] "pseries" "numeric"
> LandX <- vegX$Land; class(LandX)
[1] "pseries" "numeric"
> CredAgX <- vegX$CredAg; class(CredAgX)
[1] "pseries" "numeric"
但是我在应用slag函数时出现错误:
DistX <- slag(agSPX$Dist, listw = W)
Error in lag.listw(listw, xt) : object lengths differ
我的面板数据有5年276个地区。所以，对象的特征是：
> length(DistX)
[1] 1380
> length(W)
[1] 3
> length(W$weights)
[1] 276
所以，我想知道，如果我可以在矩阵中转换 W$weights，例如 usaww 用作 slag 函数的示例，我可以应用函数 mat2listw然后用渣过X.
谁能告诉我哪里错了？
可能这不是最佳解决方案，但我使用这些 步骤:
计算了 SLX、SDM 和 SDEM 模型
1) 通过以下方式加载形状文件：
TCAL <- readOGR(dsn = ".", "Municipios_csv_BIO")
coords <- coordinates(TCAL)
2) 通过创建权重矩阵W:
contnbQueen <- poly2nb(TCAL, queen = TRUE)
contnbk4 <- knn2nb(knearneigh(coords, k = 4, RANN = FALSE))
3) 选择要应用的矩阵：
W <- nb2listw(contnbk4, glist = NULL, style = "W")
4) 在 pdata.frame 中转换 data.frame:
vegSPX <- pdata.frame(vegPainel, index = c("ID","Ano"))
5)创建特定的pdata.frame空间向量，例如：
vegIDDX <- vegSPX$IDD
vegSoilX <- vegSPX$Soil
vegQAIX <- vegSPX$QAI
6) 为这些模型指定公式：
fmSPvegX <- Area ~ IDD + Soil + QAI + slag(vegIDDX, listw = W) + slag(vegSoilX, listw = W) + slag(vegQAIX, listw = W)
7) 对原始数据应用 plm 和 splm 函数（data.table data.frame 对象）：
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)
希望对您有所帮助。
您还可以将滞后解释变量添加到数据中。
一个可重现的例子：
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
您可以验证空间解释是否正确。例如：
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
我需要单独估计空间滞后 X (SLX) 的空间计量模型，结合空间自回归模型 (SAR) 或空间误差模型 (SEM)。当它们组合在一起时，它们被称为空间 Durbin 模型 (SDM) 或空间 Durbin 误差模型 (SDEM)，遵循 Vega & Elhorst (2015) 的论文 "The SLX Model".
我打算使用 splm 包估计 R 中的所有空间面板模型，这也需要 spdep 函数。从这个意义上说，我从形状文件创建了 Queen 类型和 k = 4 的邻居列表：
> TCAL <- readOGR(dsn = ".", "Municipios_csv")
> coords <- coordinates(TCAL)
> contnbQueen <- poly2nb(TCAL, queen = TRUE)
> enter code herecontnbk4 <- knn2nb(knearneigh(coords, k = 4, RANN = FALSE))
然后我将这个邻居列表转换成权重矩阵：
> 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")
下一步，我为面板 SAR 和 SEM 模型创建了一个公式，该公式运行良好并生成了估计值：
> 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)
然后，我尝试通过创建协变量 X 的空间滞后来估计 SLX、SDM 和 SDEM 模型：
> vegX <- pdata.frame(veg, index = c("Mun","Year")); class(vegX)
[1] "pdata.frame" "data.frame"
然后我创建了 pseries 值：
> DistX <- vegX$Dist; class(DistX)
[1] "pseries" "numeric"
> LandX <- vegX$Land; class(LandX)
[1] "pseries" "numeric"
> CredAgX <- vegX$CredAg; class(CredAgX)
[1] "pseries" "numeric"
但是我在应用slag函数时出现错误:
DistX <- slag(agSPX$Dist, listw = W)
Error in lag.listw(listw, xt) : object lengths differ
我的面板数据有5年276个地区。所以，对象的特征是：
> length(DistX)
[1] 1380
> length(W)
[1] 3
> length(W$weights)
[1] 276
所以，我想知道，如果我可以在矩阵中转换 W$weights，例如 usaww 用作 slag 函数的示例，我可以应用函数 mat2listw然后用渣过X.
谁能告诉我哪里错了？
可能这不是最佳解决方案，但我使用这些 步骤:
计算了 SLX、SDM 和 SDEM 模型1) 通过以下方式加载形状文件：
TCAL <- readOGR(dsn = ".", "Municipios_csv_BIO")
coords <- coordinates(TCAL)
2) 通过创建权重矩阵W:
contnbQueen <- poly2nb(TCAL, queen = TRUE)
contnbk4 <- knn2nb(knearneigh(coords, k = 4, RANN = FALSE))
3) 选择要应用的矩阵：
W <- nb2listw(contnbk4, glist = NULL, style = "W")
4) 在 pdata.frame 中转换 data.frame:
vegSPX <- pdata.frame(vegPainel, index = c("ID","Ano"))
5)创建特定的pdata.frame空间向量，例如：
vegIDDX <- vegSPX$IDD
vegSoilX <- vegSPX$Soil
vegQAIX <- vegSPX$QAI
6) 为这些模型指定公式：
fmSPvegX <- Area ~ IDD + Soil + QAI + slag(vegIDDX, listw = W) + slag(vegSoilX, listw = W) + slag(vegQAIX, listw = W)
7) 对原始数据应用 plm 和 splm 函数（data.table data.frame 对象）：
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)
希望对您有所帮助。
您还可以将滞后解释变量添加到数据中。
一个可重现的例子：
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
您可以验证空间解释是否正确。例如：
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