O pacote plm é responsável pelos painéis convencionais (não espaciais) que usaremos para comparação. O pacote splm é responsável pelos painéis espaciais. Os autores do pacote lançaram um artigo sobre ele neste link.
library(spatialreg)
Loading required package: spData
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
(status 2 uses the sf package in place of rgdal)
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Loading required package: Matrix
Loading required package: sf
Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
library(plm)library(splm)
Attaching package: 'splm'
The following object is masked from 'package:spatialreg':
impacts
6.2 Shapefile
# Pacoteslibrary(sf)library(sp)# Abra o arquivonatregimes_montana.shp <-st_read("data/natregimes_montana.gml")
Reading layer `natregimes_montana' from data source
`/home/raphael/projects/ecoespacial/data/natregimes_montana.gml'
using driver `GML'
Simple feature collection with 55 features and 74 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -116.0625 ymin: 44.35373 xmax: -104.0426 ymax: 49
Geodetic CRS: WGS 84
natregimes_montana.shp <-st_make_valid(natregimes_montana.shp)natregimes_montana.shp <-as_Spatial(natregimes_montana.shp)# Plotar o mapaplot(natregimes_montana.shp)
Vamos separar algumas variáveis para usarmos no modelo
dados <- natregimes_montana.shp@datadados <-subset(dados, select=c("POLY_ID", "HR90", "HR80", "RD90", "RD80","UE90", "UE80"))
6.4 Matriz de vizinhança
Para rodar os paineis espaciais, vamos precisar de uma matriz de vizinhança.
library(spdep)
Attaching package: 'spdep'
The following objects are masked from 'package:spatialreg':
get.ClusterOption, get.coresOption, get.mcOption,
get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
set.coresOption, set.mcOption, set.VerboseOption,
set.ZeroPolicyOption
Characteristics of weights list object:
Neighbour list object:
Number of regions: 55
Number of nonzero links: 272
Percentage nonzero weights: 8.991736
Average number of links: 4.945455
Link number distribution:
2 3 4 5 6 7 8 9
3 8 11 16 6 7 3 1
3 least connected regions:
1 7 28 with 2 links
1 most connected region:
2 with 9 links
Weights style: W
Weights constants summary:
n nn S0 S1 S2
W 55 3025 55 23.70729 228.2789
fe <-plm(esp, data = dadosp)re <-plm(esp, data = dadosp, model ="random")ph <-phtest(fe, re) # H0: efeitos aleatóriosprint(ph)
Hausman Test
data: esp
chisq = 18.503, df = 2, p-value = 9.598e-05
alternative hypothesis: one model is inconsistent
6.7.2 Teste espacial de Hausman
error_type <-"b"sar_random <-spml(esp, data = dadosp, listw = w1, lag =TRUE, model ="random", effect ="individual", spatial.error ="none")sar_fixed <-spml(esp, data = dadosp, listw = w1, lag =TRUE, model ="within", effect ="individual", spatial.error ="none")sem_random <-spml(esp, data = dadosp, listw = w1, lag =FALSE, model ="random", effect ="individual", spatial.error = error_type)sem_fixed <-spml(esp, data = dadosp, listw = w1, lag =FALSE, model ="within", effect ="individual", spatial.error = error_type)sac_random <-spml(esp, data = dadosp, listw = w1, lag =TRUE, model ="random", effect ="individual", spatial.error = error_type)sac_fixed <-spml(esp, data = dadosp, listw = w1, lag =TRUE, model ="within", effect ="individual", spatial.error = error_type)test_sar <-sphtest(sar_random, sar_fixed)test_sem <-sphtest(sem_random, sem_fixed)test_sac <-sphtest(sac_random, sac_fixed)res <-cbind(c(test_sar$statistic, test_sar$p.value),c(test_sem$statistic, test_sem$p.value),c(test_sac$statistic, test_sac$p.value))dimnames(res) <-list(c("test", "p-value"), c("SAR","SEM","SAC"))round(x = res, digits =5)
SAR SEM SAC
test 57.60107 59.82911 73.12357
p-value 0.00000 0.00000 0.00000
6.7.3 Teste Pesaran CD (cross-section dependence)
cd <-pcdtest(esp, data = dadosp) # H0: ausência de dependência CS
Warning: Insufficient number of observations in time to estimate heterogeneous
model: using within residuals
print(cd)
Pesaran CD test for cross-sectional dependence in panels
data: HR ~ RD + UE
z = 0.47708, p-value = 0.6333
alternative hypothesis: cross-sectional dependence
6.7.4 Teste BSK
res_lmh <-bsktest(esp, data = dadosp, listw = w1, test ="LMH")res_lm1 <-bsktest(esp, data = dadosp, listw = w1, test ="LM1")res_lm2 <-bsktest(esp, data = dadosp, listw = w1, test ="LM2")res_clm_mu <-bsktest(esp, data = dadosp, listw = w1, test ="CLMmu")res_clm_lambda <-bsktest(esp, data = dadosp, listw = w1, test ="CLMlambda")res <-cbind(c(res_lmh$statistic, res_lmh$p.value),c(res_lm1$statistic, res_lm1$p.value),c(res_lm2$statistic, res_lm2$p.value),c(res_clm_mu$statistic, res_clm_mu$p.value),c(res_clm_lambda$statistic, res_clm_lambda$p.value))dimnames(res) <-list(c("test", "p-value"), c("LM joint","LM mu","LM lambda", "CLM mu", "CLM lambda"))round(x = res, digits =5)
LM joint LM mu LM lambda CLM mu CLM lambda
test 0.00067 -0.16962 0.02596 0.16929 0.01471
p-value 0.73956 1.13469 0.97929 0.86557 0.98826
6.7.5 Teste BSJK
# res_c1 <- bsjktest(esp, data = dadosp, listw = w1, test = "C.1")# res_c2 <- bsjktest(esp, data = dadosp, listw = w1, test = "C.2")# res_c3 <- bsjktest(esp, data = dadosp, listw = w1, test = "C.3")# res_j <- bsjktest(esp, data = dadosp, listw = w1, test = "J")# # res <- cbind(# c(res_c1$statistic, res_c1$p.value),# c(res_c2$statistic, res_c2$p.value),# c(res_c3$statistic, res_c3$p.value),# c(res_j$statistic, res_j$p.value)# )# # dimnames(res) <- list(c("test", "p-value"), c("C.1","C.2","C.3", "J"))# round(x = res, digits = 5)
6.8 Modelos
6.8.1 OLS
modOLS <-plm(esp, data=dadosp)summary(modOLS)
Oneway (individual) effect Within Model
Call:
plm(formula = esp, data = dadosp)
Balanced Panel: n = 55, T = 2, N = 110
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-1.6018e+01 -1.4836e+00 -9.7145e-17 1.4836e+00 1.6018e+01
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
RD -1.70968 2.03753 -0.8391 0.40518
UE -0.60668 0.30327 -2.0005 0.05058 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 1945.6
Residual Sum of Squares: 1778.9
R-Squared: 0.085641
Adj. R-Squared: -0.88048
F-statistic: 2.48204 on 2 and 53 DF, p-value: 0.093239
6.8.2 SAR (ML)
modSAR <-spml(esp, data = dadosp, listw = w1, lag =TRUE, model ="within", effect ="individual", spatial.error ="none")summary(modSAR)
