1.Láy dữ liệu map online
1.1 lệnh getData chuẩn bị loại bỏ
library(raster)
## Loading required package: sp
vietnam <-getData("GADM", country="Vietnam", level=1)
## Warning in getData("GADM", country = "Vietnam", level = 1): getData will be removed in a future version of raster
## . You can use access these data, and more, with functions from the geodata package instead
#head(vietnam,63)
library(ggplot2)
ggplot() + geom_polygon(data=vietnam,aes(x=long,y=lat,group=group,fill=id), show.legend=F)
## Regions defined for each Polygons
viet <-fortify(vietnam)
## Regions defined for each Polygons
head(viet)
## long lat order hole piece id group
## 1 105.3745 10.24604 1 FALSE 1 1 1.1
## 2 105.3362 10.23442 2 FALSE 1 1 1.1
## 3 105.3154 10.26842 3 FALSE 1 1 1.1
## 4 105.3120 10.27393 4 FALSE 1 1 1.1
## 5 105.3065 10.26894 5 FALSE 1 1 1.1
## 6 105.3013 10.27698 6 FALSE 1 1 1.1
ggplot() + geom_polygon(data=viet,aes(x=long,y=lat,group=group,fill=factor(group)), show.legend=F)
chon = c("Vietnam", "Thailand", "Cambodia", "Malaysia", "Laos")
indo = map_data("world", region = chon)
ggplot(indo, aes(x=long, y=lat, group=group, fill=factor(group))) + geom_polygon(col="white") + theme(legend.position="none")
1.2 Vẽ theo lệnh mới
#library(geodata)
#vn <- gadm(country="VNM", level=1, path=tempdir())
#head(vn)
#plot(vn)
#head(vn)
#vn <-fortify(vn)
#tail(vn)
#ggplot() + geom_polygon(data=vn,aes(x=lon,y=lat,group=group,fill=id), show.legend=F)
#ggplot(data=vn, aes(x=long, y=lat, group=group)) + geom_polygon(aes(fill=id), col="grey30", show.legend=F)
#ggplot() + geom_sf(data=vn,aes(x=long,y=lat,group=group,fill=id), show.legend=F)
2. Lấy dữ liệu map offline
library(spdep)
## Loading required package: spData
## 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: sf
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(maptools)
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
setwd("/Users/thuphan/Desktop/hocplotly/VNM_adm")
bando <- readShapePoly('VNM_adm1.shp')
## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read
#head(bando, 63)
#plot(bando)
mtay <- bando[c(3,7,13,17,18,20,31,33,39,51,58,59,61),]
head(mtay,13)
## ID_0 ISO NAME_0 ID_1 NAME_1 TYPE_1
## 2 250 VNM Vietnam 3 Đồng Tháp Tỉnh
## 6 250 VNM Vietnam 7 An Giang Tỉnh
## 12 250 VNM Vietnam 13 Bạc Liêu Tỉnh
## 16 250 VNM Vietnam 17 Bến Tre Tỉnh
## 17 250 VNM Vietnam 18 Cà Mau Tỉnh
## 19 250 VNM Vietnam 20 Cần Thơ Thành phố trực thuộc tỉnh
## 30 250 VNM Vietnam 31 Hậu Giang Tỉnh
## 32 250 VNM Vietnam 33 Kiên Giang Tỉnh
## 38 250 VNM Vietnam 39 Long An Tỉnh
## 50 250 VNM Vietnam 51 Sóc Trăng Tỉnh
## 57 250 VNM Vietnam 58 Tiền Giang Tỉnh
## 58 250 VNM Vietnam 59 Trà Vinh Tỉnh
## 60 250 VNM Vietnam 61 Vĩnh Long Tỉnh
## ENGTYPE_1 NL_NAME_1 VARNAME_1
## 2 Province <NA> Dong Thap
## 6 Province <NA> An Giang
## 12 Province <NA> Bac Lieu
## 16 Province <NA> Ben Tre
## 17 Province <NA> Ca Mau
## 19 City|Municipality|Thanh Pho <NA> Can Tho
## 30 Province <NA> Hau Giang
## 32 Province <NA> Kien Giang
## 38 Province <NA> Long An
## 50 Province <NA> Soc Trang
## 57 Province <NA> Tien Giang
## 58 Province <NA> Tra Vinh
## 60 Province <NA> Vinh Long
plot(mtay)
class(mtay)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
ggplot() + geom_polygon(data=mtay,aes(x=long,y=lat,group=group,fill=id),color = "white", show.legend=F)
## Regions defined for each Polygons
3. Tạo dữ liệu cho mô hình không gian
library(readxl)
setwd("/Users/thuphan/Desktop/hocplotly")
dulieu <-read_excel("dmtay.xlsx")
head(dulieu)
## # A tibble: 6 × 9
## ID_1 id la lo Gdp Lab Cap Inv Mat
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 2 10.7 106. 9.46 9.31 8.86 8.50 5.23
## 2 7 6 10.5 105. 9.27 8.10 7.84 6.77 6.99
## 3 13 12 9.27 106. 11.1 9.91 8.51 8.89 7.65
## 4 17 16 10.2 106. 9.48 10.6 7.48 8.96 8.73
## 5 18 17 9.17 105. 12.0 10.8 10.9 8.16 7.71
## 6 20 19 10.1 106. 15.6 14.1 11.7 8.67 6.61
SpatialPolygonsDataFrame(mtay,data=data.frame(dulieu, row.names=row.names(mtay))) ->mtay2
class(mtay2)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
mtay2
## class : SpatialPolygonsDataFrame
## features : 13
## extent : 103.4582, 106.7919, 8.563332, 11.03249 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 9
## names : ID_1, id, la, lo, Gdp, Lab, Cap, Inv, Mat
## min values : 3, 2, 9.16722337493423, 105.069781224349, 9.10438714828342, 8.09985558751294, 7.48194800830551, 6.2188316421091, 5.23218109901466
## max values : 61, 60, 10.7273877574565, 106.376752166402, 15.6282472137827, 14.0767291363272, 11.6761749519313, 9.64245341832953, 8.73663711323625
plot(mtay2)
4. Tạo dữ liệu cho vẽ đồ thị
#Giải mã polygon thành data frame
mtay <-fortify(mtay)
## Regions defined for each Polygons
class(mtay)
## [1] "data.frame"
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
mtay %>% select(id)%>%table()
## id
## 12 16 17 19 2 30 32 38 50 57 58 6 60
## 4658 2659 7254 837 1828 1870 14156 3675 3656 2227 2734 1378 1210
mtaydata <- merge(mtay, dulieu, by=c("id"), all=TRUE)
ggplot(data=mtaydata) +
geom_polygon(aes(x=long,y=lat,group=group,fill=id),color = "white", show.legend=F) +
geom_point(aes(x=lo, y=la), col="#8C3F4D")
5.Test cho Mô hình không gian
congthuc <- Gdp ~Lab + Cap + Inv + Mat
hoiquyols <-lm(data=mtay2,congthuc)
summary(hoiquyols)
##
## Call:
## lm(formula = congthuc, data = mtay2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0517 -0.7393 -0.2186 0.9411 1.2826
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0306 3.7664 0.539 0.60448
## Lab 1.1186 0.2578 4.338 0.00248 **
## Cap 0.4826 0.3398 1.420 0.19330
## Inv -0.8909 0.3390 -2.628 0.03029 *
## Mat 0.1023 0.2944 0.348 0.73718
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.029 on 8 degrees of freedom
## Multiple R-squared: 0.8616, Adjusted R-squared: 0.7924
## F-statistic: 12.45 on 4 and 8 DF, p-value: 0.001631
5.1 Ma trận trọng số
ds.queen<-poly2nb(mtay2, queen=TRUE)
W<-nb2listw(ds.queen, style="W", zero.policy=TRUE)
W
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 13
## Number of nonzero links: 52
## Percentage nonzero weights: 30.76923
## Average number of links: 4
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 13 169 13 6.898175 53.9846
plot(W,coordinates(mtay2))
5.2 Moran’s I Test
moran.lm<-lm.morantest(hoiquyols, W, alternative="two.sided")
print(moran.lm)
##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## Moran I statistic standard deviate = -0.83429, p-value = 0.4041
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I Expectation Variance
## -0.23438639 -0.08279977 0.03301290
5.3 Lagrange Multiplier Test
LM<-lm.LMtests(hoiquyols, W, test="all")
print(LM)
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## LMerr = 1.3459, df = 1, p-value = 0.246
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## LMlag = 1.1341, df = 1, p-value = 0.2869
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## RLMerr = 0.3042, df = 1, p-value = 0.5813
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## RLMlag = 0.092372, df = 1, p-value = 0.7612
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## SARMA = 1.4383, df = 2, p-value = 0.4872
6. Chạy mô hình không gian
6.1 SAR model
library(spatialreg )
## Loading required package: Matrix
##
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption
sar<- lagsarlm(congthuc, data=mtay2, W)
summary(sar)
##
## Call:lagsarlm(formula = congthuc, data = mtay2, listw = W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.93915 -0.67135 -0.15593 0.67169 1.39147
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.271461 5.472697 1.3287 0.1839537
## Lab 1.111048 0.189793 5.8540 4.799e-09
## Cap 0.451822 0.256932 1.7585 0.0786578
## Inv -0.945038 0.251771 -3.7536 0.0001743
## Mat 0.057258 0.222602 0.2572 0.7970070
##
## Rho: -0.34059, LR test value: 1.3419, p-value: 0.24669
## Asymptotic standard error: 0.28095
## z-value: -1.2123, p-value: 0.22541
## Wald statistic: 1.4696, p-value: 0.22541
##
## Log likelihood: -14.98611 for lag model
## ML residual variance (sigma squared): 0.57145, (sigma: 0.75594)
## Number of observations: 13
## Number of parameters estimated: 7
## AIC: 43.972, (AIC for lm: 43.314)
## LM test for residual autocorrelation
## test value: 0.83572, p-value: 0.36062
# Độ nhạy biên
impacts(sar, listw=W)
## Impact measures (lag, exact):
## Direct Indirect Total
## Lab 1.14096099 -0.31218486 0.82877613
## Cap 0.46398636 -0.12695396 0.33703240
## Inv -0.97048123 0.26553892 -0.70494231
## Mat 0.05879989 -0.01608857 0.04271132
# Vẽ đồ thị
spplot(mtay2)
6.2 SEM model
sem<-errorsarlm(congthuc, data=mtay2, W)
summary(sem)
##
## Call:errorsarlm(formula = congthuc, data = mtay2, listw = W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.948490 -0.431011 0.019093 0.357872 0.870181
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.32718 2.00122 -0.1635 0.87013
## Lab 0.85519 0.17248 4.9581 7.120e-07
## Cap 1.14364 0.29207 3.9156 9.018e-05
## Inv -1.19339 0.22728 -5.2508 1.514e-07
## Mat 0.29671 0.16159 1.8362 0.06632
##
## Lambda: -1.2872, LR test value: 4.1108, p-value: 0.042611
## Asymptotic standard error: 0.23651
## z-value: -5.4424, p-value: 5.2556e-08
## Wald statistic: 29.62, p-value: 5.2556e-08
##
## Log likelihood: -13.60168 for error model
## ML residual variance (sigma squared): 0.30512, (sigma: 0.55238)
## Number of observations: 13
## Number of parameters estimated: 7
## AIC: 41.203, (AIC for lm: 43.314)
6.3 SARMA model
sarma <-sacsarlm(congthuc,data=mtay2,W)
summary(sarma)
##
## Call:sacsarlm(formula = congthuc, data = mtay2, listw = W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.936553 -0.424806 0.061752 0.368121 0.827903
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.94139 3.78841 -0.2485 0.80375
## Lab 0.84778 0.17795 4.7641 1.897e-06
## Cap 1.14719 0.29319 3.9128 9.123e-05
## Inv -1.17676 0.23856 -4.9328 8.107e-07
## Mat 0.29916 0.16154 1.8520 0.06403
##
## Rho: 0.042029
## Asymptotic standard error: 0.22317
## z-value: 0.18833, p-value: 0.85062
## Lambda: -1.2978
## Asymptotic standard error: 0.26035
## z-value: -4.9847, p-value: 6.2046e-07
##
## LR test value: 4.1465, p-value: 0.12578
##
## Log likelihood: -13.58382 for sac model
## ML residual variance (sigma squared): 0.30133, (sigma: 0.54893)
## Number of observations: 13
## Number of parameters estimated: 8
## AIC: 43.168, (AIC for lm: 43.314)
6.4 Spatial Durbin Model SDM
sdm <-lagsarlm(congthuc, data=mtay2,W, type="mixed")
summary(sdm)
##
## Call:lagsarlm(formula = congthuc, data = mtay2, listw = W, type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02524310 -0.19577033 0.00041074 0.34140114 0.49353990
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.87190 5.55658 -1.5966 0.1103442
## Lab 0.66325 0.15550 4.2652 1.998e-05
## Cap 1.59985 0.27505 5.8166 6.007e-09
## Inv -1.55376 0.21914 -7.0901 1.340e-12
## Mat 0.50058 0.14794 3.3836 0.0007155
## lag.Lab 0.16973 0.54199 0.3132 0.7541592
## lag.Cap 3.39334 0.65603 5.1726 2.309e-07
## lag.Inv -2.31815 0.56785 -4.0823 4.459e-05
## lag.Mat 0.85619 0.37033 2.3120 0.0207803
##
## Rho: -1.1247, LR test value: 6.4265, p-value: 0.011243
## Asymptotic standard error: 0.274
## z-value: -4.1046, p-value: 4.0497e-05
## Wald statistic: 16.848, p-value: 4.0497e-05
##
## Log likelihood: -8.799813 for mixed model
## ML residual variance (sigma squared): 0.16487, (sigma: 0.40604)
## Number of observations: 13
## Number of parameters estimated: 11
## AIC: 39.6, (AIC for lm: 44.026)
## LM test for residual autocorrelation
## test value: 3.5705, p-value: 0.058816
6.5 SLX model
slx <-lmSLX(congthuc,data=mtay2,W)
summary(slx)
##
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
## data = as.data.frame(x), weights = weights)
##
## Residuals:
## 2 6 12 16 17 19 30 32
## -0.524782 0.150500 0.658789 -0.170473 -0.053515 -0.386364 0.622427 -0.409230
## 38 50 57 58 60
## 0.701966 -1.495272 -0.007559 0.126713 0.786799
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.8744 14.9869 -0.792 0.4725
## Lab 0.7193 0.4143 1.736 0.1576
## Cap 1.4256 0.7363 1.936 0.1249
## Inv -1.5248 0.5932 -2.571 0.0619 .
## Mat 0.3918 0.4002 0.979 0.3830
## lag.Lab -1.0998 0.9524 -1.155 0.3125
## lag.Cap 2.8043 1.7695 1.585 0.1882
## lag.Inv -0.8402 1.2731 -0.660 0.5454
## lag.Mat 0.5946 0.9925 0.599 0.5814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.099 on 4 degrees of freedom
## Multiple R-squared: 0.921, Adjusted R-squared: 0.763
## F-statistic: 5.828 on 8 and 4 DF, p-value: 0.05313
6.6 Spatial Durbin Error Model (SDEM)
sdem <- errorsarlm(congthuc, data=mtay2, W, etype="mixed")
summary(sdem)
##
## Call:errorsarlm(formula = congthuc, data = mtay2, listw = W, etype = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.1624601 -0.1198215 0.0043052 0.0680737 0.2564317
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -16.150663 2.748600 -5.8760 4.204e-09
## Lab 0.917963 0.052363 17.5309 < 2.2e-16
## Cap 1.396654 0.099193 14.0802 < 2.2e-16
## Inv -1.806936 0.134720 -13.4125 < 2.2e-16
## Mat 0.461031 0.053570 8.6061 < 2.2e-16
## lag.Lab -2.892750 0.263370 -10.9836 < 2.2e-16
## lag.Cap 5.182253 0.494689 10.4758 < 2.2e-16
## lag.Inv -1.255192 0.156134 -8.0392 8.882e-16
## lag.Mat 1.159281 0.180924 6.4075 1.479e-10
##
## Lambda: -1.7239, LR test value: 26.147, p-value: 3.1633e-07
## Asymptotic standard error: 0.022676
## z-value: -76.023, p-value: < 2.22e-16
## Wald statistic: 5779.5, p-value: < 2.22e-16
##
## Log likelihood: 1.060602 for error model
## ML residual variance (sigma squared): 0.014171, (sigma: 0.11904)
## Number of observations: 13
## Number of parameters estimated: 11
## AIC: 19.879, (AIC for lm: 44.026)
7 lựa chọn model
# Lựa chọn ols với sar
lm.morantest(hoiquyols,W)
##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## Moran I statistic standard deviate = -0.83429, p-value = 0.7979
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## -0.23438639 -0.08279977 0.03301290
# Lựa chọn sar với sem
lm.LMtests(hoiquyols,W,test="LMerr")
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## LMErr = 1.3459, df = 1, p-value = 0.246
lm.LMtests(hoiquyols,W,test="LMlag")
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## LMlag = 1.1341, df = 1, p-value = 0.2869
lm.LMtests(hoiquyols,W,test="RLMerr")
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## RLMerr = 0.3042, df = 1, p-value = 0.5813
lm.LMtests(hoiquyols,W,test="RLMlag")
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = congthuc, data = mtay2)
## weights: W
##
## RLMlag = 0.092372, df = 1, p-value = 0.7612
Hausman.test(sem)
##
## Spatial Hausman test (asymptotic)
##
## data: NULL
## Hausman test = 4.5458, df = 5, p-value = 0.4738
# Lựa chọn sem với sdem
#Likelihood ratio test (LR.Sarm(x,y))
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
LR.Sarlm(sem,sdem)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = -29.325, df = 4, p-value = 6.716e-06
## sample estimates:
## Log likelihood of sem Log likelihood of sdem
## -13.601677 1.060602
lrtest(sem,sdem)
## Likelihood ratio test
##
## Model 1: Gdp ~ Lab + Cap + Inv + Mat
## Model 2: Gdp ~ Lab + Cap + Inv + Mat
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 7 -13.6017
## 2 11 1.0606 4 29.325 6.716e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# lựa chọn sar với sdm
LR.Sarlm(sar,sdm)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = -12.373, df = 4, p-value = 0.01479
## sample estimates:
## Log likelihood of sar Log likelihood of sdm
## -14.986106 -8.799813
lrtest(sar,sdm)
## Likelihood ratio test
##
## Model 1: Gdp ~ Lab + Cap + Inv + Mat
## Model 2: Gdp ~ Lab + Cap + Inv + Mat
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 7 -14.9861
## 2 11 -8.7998 4 12.373 0.01479 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1