地理加權回歸模型(GWR, Geographically Weighted Regression)

  GWR是局部模型，公式如下所示。回歸係數隨位置發生變化，針對每個給定的位置分別進行求解。求解時確定以下三件事：確定帶寬（即以任意一點為中心，鄰域的範圍）；鄰域對該點影響的大小即權重，通過核函數確定。遵循原則為越近的數據點所待估計點的權重越高；距離類型，如歐式距離。

常用的核函數如下所示，其中參數b為核函數的帶寬（bandwidth）。

帶寬選擇：使偏差和方差達到均衡最優

‍交叉驗證（cross-validation,CV):用已有觀測值位置的估算y值減去觀測值，對差值求平方和。找到差值平方和最小的b。    

赤池信息量準則（Akaike information criterion，AIC）。AIC越小越好。‍

library(GWmodel)

#加載數據包

data(DubVoter)

#線性回歸模型

lm.global<- lm(GenEl2004 ~ DiffAdd + LARent + SC1 + Unempl +

               LowEduc +Age18_24 + Age25_44 +Age45_64,

                data = Dub.voter)

summary(lm.global)

#地理加權回歸帶寬確定

bw.gwr.1<- bw.gwr(GenEl2004 ~ DiffAdd + LARent + SC1 + Unempl +

                     LowEduc +Age18_24 +Age25_44 + Age45_64,

                   data = Dub.voter, approach ="AICc",

                   kernel ="bisquare", adaptive = TRUE)

bw.gwr.1

#執行地理加權回歸

gwr.res<- gwr.basic(GenEl2004 ~ DiffAdd + LARent + SC1 + Unempl +

                   LowEduc + Age18_24 +Age25_44 + Age45_64,

                   data = Dub.voter,bw =bw.gwr.1, kernel = "bisquare",

            adaptive = TRUE, F123.test = TRUE)

#gwr.res為結果輸出

#其中每條記錄的截距和係數保存在SDF

#SDF為sp對象

gwr.res

#-以LowEduc係數可視化為例

names(gwr.res$SDF)

library(RColorBrewer)

mypalette.6<- brewer.pal(6, "Spectral")

spplot(gwr.res$SDF,"LowEduc", key.space = "right",

       col.regions = mypalette.6, at = c(-8,-6, -4, -2, 0, 2, 4),

        main = "Basic GW regressioncoefficient estimates for LowEduc")

#---根據箱線圖原理去除異常值-

rgwr.res<- gwr.robust(GenEl2004 ~ DiffAdd + LARent + SC1 + Unempl

                 +LowEduc + Age18_24 + Age25_44 +Age45_64,

               data = Dub.voter, bw = bw.gwr.1,

               kernel = "bisquare", adaptive =TRUE, F123.test = TRUE)

spplot(rgwr.res$SDF,"LowEduc", key.space = "right",

      col.regions = mypalette.6, at = c(-8, -6,-4, -2, 0, 2, 4),

       main = "Robust GW regressioncoefficient estimates for LowEduc")

#- 共線性----

DeVar<- "GenEl2004"

InDeVars<- c("DiffAdd"," LARent", "SC1","Unempl",

             "LowEduc","Age18_24", "Age25_44",

              "Age45_64")

#過多的自變量可能導致局部共線性問題

#逐步回歸，前向選擇

model.sel<- model.selection.gwr(DeVar ,InDeVars,

                                 data = Dub.voter,

                                kernel ="bisquare",

                           adaptive = TRUE, bw= 80)

#對AICc排序

sorted.models<- model.sort.gwr(model.sel, 

numVars = length(InDeVars), ruler.vector = model.sel[[2]][,2])，

model.list<- sorted.models[[1]]

#對GWR的逐步回歸過程可視化

model.view.gwr(DeVar,InDeVars, model.list = model.list)

plot(sorted.models[[2]][,2],col = "black", pch = 20, lty = 5,

      main = "Alternative view of GWRmodel selection procedure", 

      ylab = "AICc",

    xlab = "Model number", type ="b")

#獲取行數

nobs<- dim(Dub.voter)[1]

#bw為行數即為普通線性回歸模型

#條件數

lcrm1<- gwr.lcr(GenEl2004 ~ DiffAdd + LARent + SC1 + Unempl + LowEduc +Age18_24 + Age25_44 +Age45_64, 

data = Dub.voter, bw = nobs,

   kernel = "boxcar",adaptive=TRUE)

summary(lcrm1$SDF$Local_CN)

#局部條件數----

#----cv方法確定帶寬與bw.gwr計算結果相同----

lcrm2.bw<- bw.gwr.lcr(GenEl2004 ~ DiffAdd + LARent + SC1 +    Unempl +LowEduc + Age18_24 +Age25_44 + Age45_64,

data = Dub.voter, kernel ="bisquare", adaptive=TRUE)

lcrm2.bw

lcrm2<- gwr.lcr(GenEl2004 ~ DiffAdd + LARent + SC1 + Unempl +    LowEduc +Age18_24 + Age25_44 +Age45_64,

            data = Dub.voter, bw = lcrm2.bw,

             kernel ="bisquare", adaptive = TRUE)

summary(lcrm2$SDF$Local_CN)

mypalette.7<- brewer.pal(8, "Reds")

spplot(lcrm2$SDF,"Local_CN", key.space = "right",

        col.regions = mypalette.7, cuts=7,

       main="Local condition numbers frombasic GW regression")

#---對條件數>30的地區，應用嶺回歸.--

#查找帶寬。

lcrm3.bw<- bw.gwr.lcr(GenEl2004 ~ DiffAdd + LARent + SC1 +  Unempl + LowEduc + Age18_24 +Age25_44 + Age45_64,

      data = Dub.voter, kernel= "bisquare",

              adaptive = TRUE, lambda.adjust =TRUE, cn.thresh = 30)

lcrm3.bw

lcrm3<- gwr.lcr(GenEl2004 ~ DiffAdd + LARent + SC1+ Unempl + LowEduc +Age18_24 + Age25_44 + Age45_64,

data=Dub.voter, bw = lcrm3.bw,kernel ="bisquare",

 adaptive = TRUE, lambda.adjust= TRUE, cn.thresh = 30)

summary(lcrm3$SDF$Local_CN)

summary(lcrm3$SDF$Local_Lambda)

spplot(lcrm3$SDF,"Local_CN", key.space = "right",

       col.regions = mypalette.7, cuts = 7,

       main = "Local condition numbersbefore adjustment")

spplot(lcrm3$SDF,"Local_Lambda", key.space = "right",

       col.regions = mypalette.7,cuts = 7,

      main = "Local ridge terms for LCR GW regression")