SpatialEconometrics2

Tải dữ liệu

library(ggplot2)
setwd("/Users/thuphan/Desktop/hocplotly")
#library(rio)
#dulieu <- import("SpatialEconometrics.dta")
#head(dulieu)
#class(dulieu)
#plot(dulieu)
library(spdep)

## Loading required package: sp

## 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'.

library(ggthemes)
setwd("/Users/thuphan/Desktop/hocplotly/VNM_adm")
bando <- readShapePoly('VNM_adm1.shp')

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

mtay <- bando[c(3,7,13,17,18,20,31,33,39,51,58,59,61),]
head(mtay,13)

## class       : SpatialPolygonsDataFrame 
## features    : 13 
## extent      : 103.4582, 106.7919, 8.563332, 11.03249  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 9
## names       : ID_0, ISO,  NAME_0, ID_1,    NAME_1,                    TYPE_1,                   ENGTYPE_1, NL_NAME_1, VARNAME_1 
## min values  :  250, VNM, Vietnam,    3,  An Giang, Thành phố trực thuộc tỉnh, City|Municipality|Thanh Pho,        NA,  An Giang 
## max values  :  250, VNM, Vietnam,   61, Vĩnh Long,                      Tỉnh,                    Province,        NA, Vinh Long

class(mtay)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

ggplot() + geom_polygon(data=mtay,aes(x=long,y=lat,group=group,fill=id),color = "red", show.legend=F) +
           labs(x = NULL, y =  NULL) +
           theme_hc() ->m

## Regions defined for each Polygons

m

library(readxl)
setwd("/Users/thuphan/Desktop/hocplotly")
dulieu <-read_excel("dmtay.xlsx")
mtay2 <- SpatialPolygonsDataFrame(mtay,data=data.frame(dulieu, row.names=row.names(mtay)))


library(sf)
mtay2sf <- sf::st_as_sf(mtay2)
ggplot(mtay2sf) +
  geom_sf(aes(fill = Gdp), color = "red", show.legend=F) +
  scale_fill_gradient(low="green", high="blue")+
  labs(x = NULL, y =  NULL) +
  theme(legend.position = "none")+
  theme_wsj()

3. Vẽ đồ thị với tmap

library(tmap) # vẽ bản đồ bằng hàm qtm()
qtm(mtay2,"Gdp")

## Warning: Currect projection of shape mtay2 unknown. Long-lat (WGS84) is assumed.

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(mtay2) +
    tm_polygons("Gdp") +
    tm_borders("white", lwd = .5) +
    tm_symbols(col = "red", size = "Lab", scale = .5) +
    tm_text("id", size = "Gdp", col="blue", ymod=0.7) +
    tm_legend(show = FALSE)

## Warning: Currect projection of shape mtay2 unknown. Long-lat (WGS84) is assumed.

## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).

## Text size will be constant in view mode. Set tm_view(text.size.variable = TRUE) to enable variable text sizes.

## Legend for symbol sizes not available in view mode.

2 Chuẩn bị dữ liệu ma trận W

danhsachqueen<-poly2nb(mtay2, queen=TRUE)

W<-nb2listw(danhsachqueen, 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(mtay2)
plot(W,coordinates(mtay2), add=T, col="blue")

nb2listw(danhsachqueen, style="B", zero.policy=TRUE) ->B
plot(B,coordinates(mtay2))
     
nb2listw(danhsachqueen, style="C", zero.policy=TRUE) ->C
plot(C,coordinates(mtay2))

nb2listw(danhsachqueen, style="U", zero.policy=TRUE) ->U
plot(U,coordinates(mtay2))
     
nb2listw(danhsachqueen, style="S", zero.policy=TRUE) -> S
plot(S,coordinates(mtay2))
     
nb2listw(danhsachqueen, style="minmax", zero.policy=TRUE)  ->MM
plot(MM,coordinates(mtay2))

3 Spatial Stochastic Frontier Analysis

library(ssfa)

## Loading required package: Matrix

## Loading required package: maxLik

## Loading required package: miscTools

## 
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
## 
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/

## Loading required package: spatialreg

## 
## 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

congthuc <- Gdp ~Lab +  Cap  + Inv +  Mat