GWR在R里面可以用好几个不同的包来实现，其中步骤最简单的是spgwr。思路就两步：建立窗口、用窗口扫全局。这其实就是GWR本质上的两步。比如我要在全美国范围内统计某两个（或多个）变量之间的回归关系，我可以做一个全局回归（global regression），但因为这些变量在空间分布上或许会有异质性（heterogeneity），表现在统计结果上就是空间不稳定性（nonstationarity），因此只看全局的统计，可能看不出什么结果来。举个不完全恰当但是很容易领会精神的例子，你比如说，我要分析亚洲范围内，经济发展程度与牛肉销量之间的关系，经济越发达的地方，人们就越吃得起牛肉。可是等我统计到印度的时候，坏了，印度大部分人不吃牛肉，这不是经济状况导致的，这一下就影响了全局统计的参考价值，那怎么办呢？我们可以建立一个窗口（正规说法是带宽窗口，bandwidth window），每次只统计窗口范围内的经济与牛肉销量的关系，然后用这个窗口去扫过全局的范围。等统计到印度的时候，印度内部的各地和印度自己比，吃牛肉的人的比例就不会突然减少，这样就能减少这种空间不稳定性对全局统计的影响。
所以，第一步就是要建立这样一个『窗口』。当然了，首先要安装包，我们要用到的R包有：
library(spgwr)
library(rgdal)
library(sf)
library(spData)
library(sp)
library(lattice)
library(ggplot2)
library(ggthemes)
其中，spgwr是做GWR的包，rgdal是用来读取矢量要素的，sf，sp和spData都是用来处理矢量数据的，别的基本都是画图用。
以下默认你会R和GWR的基本操作。并且，以下只展现方法，不要纠结我的数据和结果，我随便找的数据，这个数据本身没有什么意义，所以做出的统计看起来很『壮观』。
我们先导入数据。这里我用的是美国本土48州各个县（county，也有翻译成郡的）的人口普查数据和农业数据，来源是ESRI Online数据库。为啥用这个数据呢？因为...我电脑里面就存了这么个可以用来做GWR的数据...
我们用rgdal读取数据，然后把它画出来看看
require(rgdal)
usa_agri <- readOGR(dsn = "~/Documents/Spatial", layer = "usa_counties")
plot(usa_agri)
会得到这个东西：
readOGR里面，dsn后面加储存shp的路径（加到文件夹为止），layer后面写shp的文件名（不加.shp）。不喜欢rgdal的同学可以不用，用maptools或者spData等别的处理shp的R包代替。不过如果用maptools，要注意处理一下参考系。
我们看一下这个shp里面的列联表都有什么：
可见，shp里面有3108个县的数据，数据有61种。然后再看data下面有什么：
总之就是各种人口普查的数据，后面截不完图，还有经济、房地产和农业之类的数据。那我们就随便选两个来当变量。我就随便挑了，因变量选AVESIZE12，即2012年各个县农场的平均占地面积。自变量选POP_SQMI，也就是人口密度（每平方英里的人口）。
现在正式建立窗口，调用的是spgwr里面的gwr.sel函数：
bw <- gwr.sel( AVE_SIZE12 ~ POP_SQMI, data = usa_agri, gweight = gwr.Gauss,
verbose = FALSE, method = "cv")
其中~前后分别是因变量和自变量。GWR里因变量只能有1个，但自变量可以选多个，如果需要多个自变量的话，就在代码POP_SQMI之后用+号连接就行。gweight是你的空间加权的函数（随空间距离增大而不断衰减的函数，衰减率由下面要提到的带宽控制），这里用的是比较常用的高斯函数，其余的还有gwr.bisquare等函数可以调用。verbose决定是否汇报制定窗口的过程。method是决定构建带宽窗口模型的方法，这里用的cv指的是cross validation，即交叉验证法，也是最常用的方法，简单说就是把数据分成不同的组，分别用不同的方法来做回归计算，计算完了之后记录下结果，然后打乱重新分组，再回归计算，再看结果，周而复始，最后看哪种计算方法的结果最靠谱，这种方法就是最优解。还有一种很常见的选择最佳拟合模型的方法是AIC optimisation法，把method后面的cv改成aic就可以用。具体AIC optimisation是什么：AIC（赤池信息准则）_百度百科。总之，空间加权函数和带宽窗口构建方法的选择是GWR里面十分重要的步骤。
以上便是固定带宽窗口的示意图。比如我在对佐治亚做GWR，这一轮的regression target是红色的这个县，根据做出来的窗口，圆圈以内的县都要被算为红色县的邻县，其权重根据高斯函数等空间权重函数来赋值，而圆圈以外的县，空间权重都赋为0。
不喜欢固定带宽窗口的同学也可以不用它，而是用符合Tobler地理学第一定律的非固定带宽邻域统计，操作方法是在gwr.sel里面加一个命令adapt = TRUE，这样的情况下，根据你设置的k邻居数，每一轮统计的时候，和本轮对象在k以内相邻的多边形的权重参数会被赋值为0到1之间的一个数，比如下图:
我在对佐治亚做GWR，这一轮的regression target是红色的这个县，那么图上标为1的县就是红色县的1阶邻县，标为2的是2阶（邻县的邻县），标为3的是3阶（邻县的邻县的邻县）。如果用非固定带宽邻域统计，k为3，那么1、2、3都被定义为红色县的邻县，它们的权重从3到1依次增加，会按比例被赋上0和1之间的值，而其它没有标注的县，权重为0。
下一步就是用前一步做出的窗口去扫过全局区域：
gwr_result <- gwr(AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwr.Gauss, hatmatrix = TRUE)
这一步如果数据量大，可能会要跑一阵，跑完之后我们看看结果里面有什么：
Call:
gwr(formula = AVE_SIZE12 ~ POP_SQMI, data = usa_agri, bandwidth = bw,
gweight = gwr.Gauss, hatmatrix = TRUE)
Kernel function: gwr.Gauss
Fixed bandwidth: 205880.3
Summary of GWR coefficient estimates at data points:
Min. 1st Qu. Median 3rd Qu.Max. Global
X.Intercept. 7.3883e+01 2.1081e+02 3.2802e+02 6.6691e+02 8.5705e+03 625.5656
POP_SQMI -8.0085e+01 -4.5983e-01 -1.4704e-01 -7.3703e-02 -2.1859e-03 -0.0426
Number of data points: 3108
Effective number of parameters (residual: 2traceS - traceS'S): 119.6193
Effective degrees of freedom (residual: 2traceS - traceS'S): 2988.381
Sigma (residual: 2traceS - traceS'S): 1048.78
Effective number of parameters (model: traceS): 84.90185
Effective degrees of freedom (model: traceS): 3023.098
Sigma (model: traceS): 1042.741
Sigma (ML): 1028.4
AICc (GWR p. 61, eq 2.33p. 96, eq. 4.21): 52109.55
AIC (GWR p. 96, eq. 4.22): 52017.7
Residual sum of squares: 3287040139
Quasi-global R2: 0.4829366
基本上你做GWR该需要的结果这里都有了。比如窗口大小（Fixed bandwidth）是205880.3，意思是前一步构建的带宽窗口是半径205.88千米的圆。Effective number of parameters显示的是你带宽窗口的大小合不合适。Sigma是残差的标准差，这个值要尽量小。Residual sum of squares（RSS）也是对拟合程度的一个评估值。最重要的是最后那个R2，越靠近1说明统计的拟合度越好。我这里面Sigma很大，R2也不是很大，因为我这里只是呈现方法，用的数据本来就是互不相干、没什么太大意义的，所以不用太纠结。如果你是真正的统计数据要来做GWR，就需要注意这些值了。
然后，我们就可以把每个县的R2画在地图上。首先，前面报告里的这些数据，比如R2，要先自己去生成的GWR结果里面去找，然后自己再算一下每个县的local R2，并把它们赋值到shp里面去：
globalR2 <- (1 - (gwr_result$results$rss/gwr_result$gTSS))
sp <- gwr_result$SDF
sf <- st_as_sf(sp)
然后就可以各显神通画地图了，我一般还是用ggplot外加ggthemes来画这种地图，毕竟最方便：
ggplot() + geom_sf(data = sf, aes(fill=localR2)) +
coord_sf() +
theme(panel.grid.major = element_line(color = "black")) +
ggtitle(paste("Local R2")) +
labs(subtitle = paste("Global R2:", round(globalR2, 6) ) )
看效果：
如果你想plot别的东西，比如gwr.e等等，就从gwr_result$SDF提取或者计算，然后改sp和sf就行。
以上是一种情况，即你拥有的是带有不同变量的多边形矢量数据。那还有一种情况，就是你现在没有多边形，而是一堆点状数据，想要用R来做GWR，应该怎么办。比如我现在导入一个储存为csv的点状数据。数据还是和上面一样的，只不过每一个县都变成了一个点。
usa_points <-
read.csv("~/Documents/Spatial/usa_counties.csv",
stringsAsFactors = FALSE)
str(usa_points) #查看它的数据结构
查看数据结构，这个csv里面有经纬度数据，分别是lat和long。
然后我们把AVE_SALE12的数据画在地图上看看：
attach(usa_points)#让head里面的lat和long可读
map = SpatialPointsDataFrame(data=usa_points, coords=cbind(long, lat))
colours = c("dark blue", "blue", "red", "dark red")
spplot(map, "AVE_SALE12", cuts=quantile(AVE_SALE12), col.regions=colours, cex=1, main = "AVE_SALE12")
这里是给分了4个组，用了四种颜色，cuts后面直接用的quantile命令。如果要自己手动分组，cuts后面可以自己改成seq( )。效果：
再画一个POP_SQMI的数据：
然后plot一个全局的看看：
plot(AVE_SALE12 ~ POP_SQMI, data=map)
好吧这个数据简直了，不过没关系，只展示过程，数据不重要233
然后尝试一下线性回归（电脑你不要揍我233）
line_fit1 = lm(AVE_SALE12 ~ POP_SQMI, data=map)
abline(line_fit1, col= "red")
summary(line_fit1)
得到：
以及：
Call:
lm(formula = AVE_SALE12 ~ POP_SQMI, data = map)
Residuals:
Min 1Q Median 3Q Max
-197406 -153151 -80701 53722 5203063
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 197413.019 5308.363 37.189 <2e-16 ***
POP_SQMI-8.763 2.867 -3.057 0.00225 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 292800 on 3106 degrees of freedom
Multiple R-squared: 0.003, Adjusted R-squared: 0.002679
F-statistic: 9.346 on 1 and 3106 DF, p-value: 0.002254
总之这个数据很扯，但是这上面该统计的基本都有了。
然后我们用地图plot一下它的residual：
resids = residuals(line_fit1)
colours = c("dark green", "green", "yellow", "orange")
map.resids = SpatialPointsDataFrame(data=data.frame(resids),coords=cbind(long, lat))
spplot(map.resids, cuts=quantile(resids), col.regions=colours,cex=1, main = "Residuals")
以上是全局的。然后我们如法炮制，构建窗口，再用窗口扫描全局。这次我们试一下adaptive bandwidth：
bw2 = gwr.sel(AVE_SALE12 ~ POP_SQMI, data=map, adapt=TRUE)
gwr.points = gwr(AVE_SALE12 ~ POP_SQMI, data=map, adapt=bw2, hatmatrix=T, se.fit=T)
得到：
Call:
gwr(formula = AVE_SALE12 ~ POP_SQMI, data = map, adapt = bw2,
hatmatrix = T, se.fit = T)
Kernel function: gwr.Gauss
Adaptive quantile: 0.001270918 (about 3 of 3108 data points)
Summary of GWR coefficient estimates at data points:
Min. 1st Qu. Median 3rd Qu.Max. Global
X.Intercept. 3.3258e+03 7.9850e+04 1.4648e+05 2.7318e+05 2.2266e+06 197413.0190
POP_SQMI -1.4418e+04 -2.1415e+02 -5.1995e+01 -3.3218e+00 6.0489e+04 -8.7635
Number of data points: 3108
Effective number of parameters (residual: 2traceS - traceS'S): 734.2409
Effective degrees of freedom (residual: 2traceS - traceS'S): 2373.759
Sigma (residual: 2traceS - traceS'S): 199028.2
Effective number of parameters (model: traceS): 521.6917
Effective degrees of freedom (model: traceS): 2586.308
Sigma (model: traceS): 190674.6
Sigma (ML): 173937.2
AICc (GWR p. 61, eq 2.33p. 96, eq. 4.21): 85082.39
AIC (GWR p. 96, eq. 4.22): 84346.86
Residual sum of squares: 9.402988e+13
Quasi-global R2: 0.6478699
看起来adaptive带宽比前面的固定带宽要好一点。
然后给它画出来：
head(gwr.points$SDF)
spplot(gwr.points$SDF, "POP_SQMI", cuts=quantile(gwr.points$SDF$POP_SQMI), col.regions=colours, cex=1, main="GWR")
然后我们再多做一步吧，查看一下哪些地区的农场面积和人口密度之间的关系是显著的，这就要用t检验，比如假设我们根据这个例子，查了双侧t分布检验，置信区间95%，发现他落在-4到4之间说明关系显著：
t = gwr.points$SDF$POP_SQMI / gwr.points$SDF$POP_SQMI_se
sig.map = SpatialPointsDataFrame(map, data.frame(t))
ramp=c("pink","light blue","pink")
breaks=c(min(t),-4,4,max(t))
spplot(sig.map, cuts=breaks, col.regions=ramp, cex=c(0.5, 0.3, 0.5), main = "t - Value")
图里面，显著的地方用浅蓝表示，具体的student t分布的理论内容可以参考http://www.zhihu.com/question/30753175 。
根据人口密度来估计的农场规模（可以和真实的农场规模地图对比）：
spplot(gwr.points $SDF, "pred",col.regions=colours, main = "Predicted Value")
标准误差（好大的误差哈哈哈）：
spplot(gwr.points $SDF, "pred.se", col.regions=colours, main = "Standard Error")
其实点状数据GWR做到这里已经差不多了，可是有个搞事的。ESRI的网上说，点状数据是不建议直接做GWR的，其实它的意思是建议要想办法把点状数据转换为多边形，就可以做了。那怎么搞呢？有两种常见思路。
第一种是做希森多边形（Shan Ye：哪些新事物的出现，引发了看似毫无关联的领域的变革？）：
先添加几个R包：
library(spatstat)
library(maptools)
library(tmap)
library(raster)
我们先找一个美国本土单个的shp多边形要素，然后把它和点状要素花在一张图上：
require(rgdal)
usa <- readOGR(dsn = "~/Documents/Spatial", layer = "usa_boundary")
tm_shape(usa) + tm_polygons() + tm_shape(map) + tm_dots()
然后我们发现因为数据来源不一样，有的点状数据落在了美国轮廓的外面。为了做希森多边形，我们必须cut一下点状数据，把落在外面的暂时删除掉。
map@bbox <- usa@bbox
然后建立希森多边形：
th <- as(dirichlet(as.ppp(map)), "SpatialPolygons")
th <- as(th ,"SpatialPolygons")
之后把点状数据列联表里的东西全都赋值到希森多边形里，再沿着美国的边缘裁剪一下：
th.z <- over(th, map, fn=mean)
th.spdf <- SpatialPolygonsDataFrame(th, th.z)
th.crop <-crop(th, usa)
搞定，然后就可以按照前面说的多边形GWR的方法做下去了。
第二种思路是建立空间网格（和这个思路一样：Shan Ye：中国历史上发生有史料记载的战争最多的区域是哪里？）
画格子有两种常见操作：
第一种，首先查一下你的数据的空间范围
bbox(usa)
bbox(map)
然后根据范围（from和to），以及你想要的网格大小（by），建立网格的x和y坐标，并组成网格。
x <- seq(from = -13883452, to = 2898563, by = 15000)
y <- seq(from = -7454985, to = 6338173, by = 15000)
xy <- expand.grid(x = x, y = y)
class(xy)
str(xy)
然后把网格转换成spatial data frame：
grid.pts<-SpatialPointsDataFrame(coords= xy, data=xy)
plot(grid.pts)
gridded(grid.pts)
class(grid.pts)
gridded(grid.pts) <- TRUE
gridded(grid.pts)
str(grid.pts)
plot(grid.pts)
然后把它和你的点数据对应，转换成spatial data frame polygon（并且一起添加id之类的）
grid <- as(grid.pts, "SpatialPolygons")
plot(grid)
str(grid)
class(grid)
summary(grid)
gridspdf <- SpatialPolygonsDataFrame(grid, data=data.frame(id=row.names(grid), row.names=row.names(grid)))
names.grd<-sapply(gridspdf@polygons, function(x) slot(x,"ID"))
text(coordinates(gridspdf), labels=sapply(slot(gridspdf, "polygons"), function(i) slot(i, "ID")), cex=0.3)
points(map)
str(gridspdf@polygons)
第二种，是直接在你的点状数据上画
extent <- extent(bbox(map)) # 格子覆盖点数据的范围
raster <- raster(extent) #画格子
dim(raster) <- c(50, 100) # x和y方向上各自要多少个格子
projection(raster) <- CRS(proj4string(map)) # 设定投影，和点数据一样
grid <- as(raster, 'SpatialPolygonsDataFrame') # 转化为多边形
tm_shape(grid) + tm_polygons() + tm_shape(map) + tm_dots(col="red", size=0.3)
然后，我们把点状要素的值给赋到格子里去。这里可以用另一个包：
library(GWmodel)
Dist_matrix <- gw.dist(dp.locat=coordinates(map),rp.locat=coordinates(grid))
用这个包建立一个距离矩阵，然后这个包里面也有可以做GWR的功能：
gwr.res <- gwr.basic(AVE_SIZE12 ~ POP_SQMI, data=map, regression.points=grid, bw=500000, dMat=Dist_matrix, kernel='gaussian')
这里的带宽窗口我手动给设成了500公里。空间权重还是用的高斯。然后结果是：
***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2018-12-13 02:50:16
Call:
gwr.basic(formula = AVE_SIZE12 ~ POP_SQMI, data = map, regression.points = g,
bw = 5e+05, kernel = "gaussian", dMat = DM)
Dependent (y) variable: AVE_SIZE12
Independent variables: POP_SQMI
Number of data points: 3108
***********************************************************************
*Results of Global Regression *
***********************************************************************
Call:
lm(formula = formula, data = data)
Residuals:
Min 1Q Median 3QMax
-702 -464 -366 -110 37326
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 625.56560 25.90080 24.152 <2e-16 ***
POP_SQMI -0.042650.01399 -3.049 0.00231 **
---Significance stars
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1429 on 3106 degrees of freedom
Multiple R-squared: 0.4984
Adjusted R-squared: 0.4663
F-statistic: 9.297 on 1 and 3106 DF, p-value: 0.002314
***Extra Diagnostic information
Residual sum of squares: 6338159253
Sigma(hat): 1428.502
AIC: 53979.52
AICc: 53979.53
***********************************************************************
* Results of Geographically Weighted Regression *
***********************************************************************
*********************Model calibration information*********************
Kernel function: gaussian
Fixed bandwidth: 5e+05
Regression points: A seperate set of regression points is used.
Distance metric: A distance matrix is specified for this model calibration.
****************Summary of GWR coefficient estimates:******************
Min.1st Qu. Median3rd Qu. Max.
Intercept 625.565604 625.565604 625.565604 625.565604 625.5656
POP_SQMI -0.042649 -0.042649
