# 空間自己回帰モデル
# 
# 
：目的変数，
:説明変数，
:空間重み行列
# 
：回帰係数，
:空間ラグ（自己相関）パラメータ
# コロンバス市の犯罪データの例（空間従属性の検定）
# oldcolにCOL.OLD(データフレーム)とCOl.nb(nbクラス)が格納
# CRIME : 1,000世帯辺りの犯罪率
# HOVAL : 家の評価額(単位は$1,000)
# INC : 世帯収入
library(spdep)
data(oldcol)
names(COL.OLD)
# 重回帰モデル
col.ols <- lm(CRIME~INC+HOVAL,data=COL.OLD)
summary(col.ols)
# spdepにあるlm.LMtests()を使って，残差に関する空間従属性の検定を行う
# LMlag : 
の検定
# LMerr : 
の検定
W <- nb2listw(COL.nb,style="W") # 区域ごとに標準化した空間重み行列
lm.LMtests(col.ols,listw=W,test=c("LMlag","LMerr"))
# 空間ARモデル
# 空間自己回帰モデルで，
とした場合，つまり，
# 
# 隣接地域の値でを説明する1次の自己回帰モデル
# Rでは，spdepパッケージのlagsarlm()を使う
# コロンバス市の犯罪データの例（空間ARモデル）
col.ar <- lagsarlm(CRIME~1,data=COL.OLD,listw=W,method="eigen")
summary(col.ar)
# 空間同時自己回帰（SAR）モデル
# 空間自己回帰モデルで，
とした場合，つまり，
# 
# Rでは，spdepパッケージのlagsarlm()を使う
# コロンバス市の犯罪データの例（SARモデル）
col.sar <- lagsarlm(CRIME~INC+HOVAL,data=COL.OLD,listw=W,method="eigen")
summary(col.sar)
# 空間誤差モデル（SEM; spatial error model）
# 空間自己回帰モデルで，
とした場合，つまり，
# 
# 重回帰モデルの誤差項に空間自己相関がある形となる
# Rでは，spdepパッケージのerrorsarlm()を使う
# コロンバス市の犯罪データの例（SEM）
col.sem <- errorsarlm(CRIME~INC+HOVAL,data=COL.OLD,listw=W,method="eigen")
summary(col.sem)
# 空間Durbinモデル(SDM)
# 
# 目的変数と説明変数の空間自己相関を混合したモデル（SAR混合モデルとも呼ばれる）
# spdepパッケージのlagsarlm(...,type="mixed")を使う
# コロンバス市の犯罪データの例（SDM）
col.sdm <- lagsarlm(CRIME~INC+HOVAL,data=COL.OLD,listw=W,method="eigen",type="mixed")
summary(col.sdm)
# anova()によるsarlmクラスの比較
anova(col.ar,col.sar,col.sem,col.sdm,col.ols)
# 条件付き自己回帰（CAR）モデル
# 空間誤差が条件付き分布
に基づいている
# Rでは，spdepパッケージのspautolm(...,family="CAR")を使う
# spautolm(...,family="SAR")とすればSARモデルの推定もできる
# コロンバス市の犯罪データの例（CARモデル）
col.car <- spautolm(CRIME~INC+HOVAL,data=COL.OLD,family="CAR",listw=W)
summary(col.car)
# 地理的加重回帰（GWR）モデル
# グローバル検定
# spgwrパッケージに実装されるグローバル検定
# LMZ.F1GWR.test()
# LMZ2.F1GWR.test()
# LMZ3.F1GWR.test() : 説明変数ごとの有意性を調べる
# BFC99.gwr.test()
# BFC02.gwr.test()
# コロンバス市の犯罪データの例（GWRモデル）
library(spgwr)
data(georgia)
f <- formula(PctBach~TotPop90+PctRural+PctEld+PctFB+PctPov+PctBlack)
g.gauss <- gwr.sel(f,gweight=gwr.Gauss,data=gSRDF) # バンド幅の選択
res <- gwr(f,gweight=gwr.Gauss,data=gSRDF,bandwidth=g.gauss) # GWRの実行
str(res,max.level=1) # resに含まれるオブジェクトの表示(第一階層まで)
class(res$SDF) # 回帰係数はSpatialPolygonsDataFrameクラスの属性値として保存されている
# 黒人人口割合(PctBlack)の回帰係数を地図上に視覚化
brks <- c(-0.025,0,0.01,0.025,0.075)
cols <- grey(5:2/6)
plot(res$SDF,col=cols[findInterval(res$SDF$PctBlack,brks,all.inside=T)])
legend("bottomleft",leglabs(brks),fill=cols,bty="n")
# GWRのグローバル検定
res <- gwr(f,gweight=gwr.Gauss,data=gSRDF,bandwidth=g.gauss,hatmatrix=T)
LMZ.F1GWR.test(res)
LMZ.F2GWR.test(res)
LMZ.F3GWR.test(res)
BFC99.gwr.test(res)
BFC02.gwr.test(res)