第12回演習マニュアル

演習の目的

  • R言語を使って,空間計量経済モデルを構築する.

空間計量経済パッケージを読み込む.

  • 以下の二つのライブラリを使って空間計量経済モデルを推定します.
    • 空間的自己相関モデル spdep
    • 地理的加重回帰モデル spgwr
> library(spdep)
要求されたパッケージ tripack をロード中です
要求されたパッケージ maptools をロード中です
要求されたパッケージ foreign をロード中です
要求されたパッケージ sp をロード中です
要求されたパッケージ SparseM をロード中です
[1] "SparseM library loaded"
要求されたパッケージ boot をロード中です
> library(spgwr)

誤差項の空間的自己相関(SEM)モデル(errorsarlm)

  • 第5回目の演習問題で用いたコロンバス市の犯罪データと同じデータを使います.但しファイル名は”oldcol”で,これを読み込んで属性データCOL.OLDと隣接関係COL.nbを用いて分析を行います.
  • 空間計量経済学のLuc Anselin教授(イリノイ大)によるAn Introduction to Spatial Regression Analysis in Rも参考にしてください.
> data(oldcol)
> summary(COL.OLD)
  • 属性の内容についてはspdep.pdfのoldcolの欄を参照してください.
  • 誤差項の空間的自己相関(SEM)モデルの推定にはerrorsarlmコマンドを使います.
> COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, nb2listw
(COL.nb, style="W"), method="eigen", quiet=FALSE)
  • ここでCRIME ~ INC + HOVALは,各地区の犯罪発生件数CRIMEを被説明変数,世帯所得水準INCと住宅価格HOVALを説明変数とした回帰モデルであることを意味します.
  • dataはモデルを適用するデータ,methodは対角行列(eigenvalue)かスパース行列(sparse)かを選べます.
  • quiet=FALSEにすると,最尤法の収束計算過程の一部を出力します.空間重み付け行列はnb2listwコマンドを使って定義していますが,その内容は第7回目の演習を参照してください.
  • 結果は,以下のように表示されます.空間的自己相関項に対するスケールパラメータλは0.562でそのZ値は4.197であることから,統計的に有意であることがわかります.
  • また,モデル全体の説明力は赤池の情報量基準AICで示されていますが,通常最小二乗法lmの結果と比較してAICが小さいため,より当てはまりがよいといえます.以下,空間的自己相関モデルの当てはまり具合を比較する際には,AICを参考にしてください.
> summary(COL.errW.eig)
Call:
errorsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = nb2listw
(COL.nb, style = "W"), method = "eigen", quiet = FALSE)
Residuals:
      Min        1Q    Median        3Q       Max 
-34.81174  -6.44031  -0.72142   7.61476  23.33626 
Type: error 
Coefficients: (asymptotic standard errors) 
             Estimate Std. Error z value  Pr(>|z|)
(Intercept) 59.893219   5.366163 11.1613 < 2.2e-16
INC         -0.941312   0.330569 -2.8476 0.0044057
HOVAL       -0.302250   0.090476 -3.3407 0.0008358
Lambda: 0.56179 LR test value: 7.9935 p-value: 0.0046945 
Asymptotic standard error: 0.13387 z-value: 4.1966 p-value: 2.7098e-05 
Wald statistic: 17.611 p-value: 2.7098e-05 
Log likelihood: -183.3805 for error model
ML residual variance (sigma squared): 95.575, (sigma: 9.7762)
Number of observations: 49 
Number of parameters estimated: 5 
AIC: 376.76, (AIC for lm: 382.75)

空間的自己回帰(SAR)モデル(lagsarlm)

  • 空間的自己回帰モデルSARには,lagsarlmコマンドを用います.
> COL.lag.eig <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, nb2listw
(COL.nb), method="eigen", quiet=FALSE)
  • 結果は以下のように示されます.
> summary(COL.lag.eig)
Call:
lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = nb2listw
(COL.nb), 
    method = "eigen", quiet = FALSE)
Residuals:
     Min        1Q    Median        3Q       Max 
-37.68585  -5.35636   0.05421   6.02013  23.20555 
Type: lag 
Coefficients: (asymptotic standard errors) 
             Estimate Std. Error z value  Pr(>|z|)
(Intercept) 45.079249   7.177346  6.2808 3.369e-10
INC         -1.031616   0.305143 -3.3808 0.0007229
HOVAL       -0.265926   0.088499 -3.0049 0.0026570
Rho: 0.43102 LR test value: 9.9736 p-value: 0.001588 
Asymptotic standard error: 0.11768 z-value: 3.6626 p-value: 0.00024962 
Wald statistic: 13.415 p-value: 0.00024962 
Log likelihood: -182.3904 for lag model
ML residual variance (sigma squared): 95.494, (sigma: 9.7721)
Number of observations: 49 
Number of parameters estimated: 5 
AIC: 374.78, (AIC for lm: 382.75)
LM test for residual autocorrelation
test value: 0.31954 p-value: 0.57188 
  • SARモデルのAICは374.78であり,通常最小二乗法(OLS)のAIC382.75,SEMモデルのAIC376.76より小さく,これらのモデルより当てはまりがよいことがわかります.

Spatial Durbinモデル(lagsarlm)

  • Spatial Durbinモデルもlagsarlmコマンドを使います.このとき,オプションtype=”mixed”を指定します.
> COL.mixed.W <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, nb2listw(COL.nb, 
style="W"), type="mixed")
  • 結果は以下のようになります.
    • lag.INCとlag.HOVALが空間的重み付け行列を持つ説明変数に対するパラメータです.
    • この場合,Z値がいずれも5%水準で統計的に有意でないほか,上述のSARモデルと比較してAICが大きいため,モデルとしての当てはまりもよくないことがわかります.
> summary(COL.mixed.W)
Call:
lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = nb2listw(COL.nb, 
   style = "W"), type = "mixed")
Residuals:
      Min        1Q    Median        3Q       Max 
-37.47829  -6.46731  -0.33835   6.05200  22.62969 
Type: mixed 
Coefficients: (asymptotic standard errors) 
             Estimate Std. Error z value  Pr(>|z|)
(Intercept) 42.822415  12.667205  3.3806 0.0007233
INC         -0.914223   0.331094 -2.7612 0.0057586
HOVAL       -0.293738   0.089212 -3.2926 0.0009927
lag.INC     -0.520284   0.565129 -0.9206 0.3572355
lag.HOVAL    0.245640   0.178917  1.3729 0.1697756
Rho: 0.42634 LR test value: 5.3693 p-value: 0.020494 
Asymptotic standard error: 0.15623 z-value: 2.7288 p-value: 0.0063561 
Wald statistic: 7.4465 p-value: 0.0063561 
Log likelihood: -181.3935 for mixed model
ML residual variance (sigma squared): 91.791, (sigma: 9.5808)
Number of observations: 49 
Number of parameters estimated: 7 
AIC: 376.79, (AIC for lm: 380.16)
LM test for residual autocorrelation
test value: 0.28919 p-value: 0.59074 

空間的自己相関モデルを使った予測(predict.sarlm)

  • 参考文献:Haining, R. 1990 Spatial data analysis in the social and environmental sciences, Cambridge:Cambridge University Press.
  • predict.sarlmコマンドは,空間的自己相関モデルの予測値を返します.
> COL.mix.eig <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, nb2listw(COL.nb), type="mixed")
> print(p1 <- predict(COL.mix.eig))
         fit      trend    signal
1  26.044311 14.8543508 11.189960
2  44.034234 29.2632112 14.771023
    :    :    :
  • ここで,モデル構造と出力結果(trend,signal)は以下のように対応します.従って,fit = trend + signal と出力されます,
モデルtrendsignal
SEMλWy-λWXβ
SAR(lag, mixed)ρWy

地理的加重回帰モデル(spgwr)

  • 同じコロンバスのデータを使って,地理的加重回帰モデルを推定してみましょう.
  • データの読み込み
> data(columbus)
  • まず,lmコマンドを用いて通常最小二乗法による重回帰モデルを推定してみましょう.
> col.lm <- lm(crime ~ income + housing, data=columbus)
  • 結果は以下のように表示されます.説明変数に対するt値はいずれも5%水準で統計的に有意ですが,R2は0.55となっています.
> summary(col.lm)
Call:
lm(formula = crime ~ income + housing, data = columbus)
Residuals:
   Min      1Q  Median      3Q     Max 
-34.418  -6.388  -1.580   9.052  28.649 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  68.6189     4.7355  14.490  < 2e-16 ***
income       -1.5973     0.3341  -4.780 1.83e-05 ***
housing      -0.2739     0.1032  -2.654   0.0109 *  
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
Residual standard error: 11.43 on 46 degrees of freedom
Multiple R-Squared: 0.5524,     Adjusted R-squared: 0.5329 
F-statistic: 28.39 on 2 and 46 DF,  p-value: 9.34e-09
  • 次に,cross varidation scoreが最小となるようなbandwidthを計算します.
> col.bw <- gwr.sel(crime ~ income + housing, data=columbus, coords=cbind
(columbus$x, columbus$y))
Loading required package: mva 
Bandwidth: 7.311181 Score: 12.22534 
Bandwidth: 11.81793 Score: 12.23804
       :       :
Bandwidth: 3.217381 Score: 11.12129
  • となり,cross varidation scoreが最小となるbandwidthが3.217となります.
  • 重み付け行列の要素としてガウス関数を使った地理的加重回帰モデルを推定します.
> col.gauss <- gwr(crime ~ income + housing, data=columbus, coords=cbind
(columbus$x, columbus$y), bandwidth=col.bw, hatmatrix=TRUE)
  • パラメータはサンプル毎に推定されますが,その基本統計量は次のようになります.
> col.gauss
Call:
gwr(formula = crime ~ income + housing, data = columbus, coords = cbind(columbus$x, 
    columbus$y), bandwidth = col.bw, hatmatrix = TRUE)
Kernel function: gwr.gauss 
Fixed bandwidth: 3.217381 
Summary of GWR coefficient estimates:
                 Min.  1st Qu.   Median     Mean  3rd Qu.    Max.
X.Intercept. 23.23000 54.13000 63.90000 58.92000 68.76000 80.9000
income       -3.13100 -1.91300 -0.98440 -1.18600 -0.36860  1.2910
housing      -1.05300 -0.37670 -0.09739 -0.14030  0.03006  0.7946
Number of data points: 49 
Effective number of parameters: 29.61663 
Effective degrees of freedom: 19.38337 
Sigma squared (ML): 25.49075 
AIC (GWR p. 96, eq. 4.21): 403.6193 
AICc (GWR p. 96, eq. 4.22): 321.6617 
Residual sum of squares: 1249.047
  • パラメータ行列をCSV形式で書き出しましょう.
> write.table(col.gauss$SDF, file="col.gauss.csv",sep=",")
  • 出力後のCSVファイルはヘッダーが一列ずつずれているので,予めずらし,一列目のヘッダーをIDなどとしてください.ArcGISでcolombus.shpを表示し,テーブルの結合を使ってcol.gauss.csvを結合させると,パラメータの空間分布が把握できます.
  • 出力される列のうち,incomeとhousingが回帰係数です(エクセルで平均値などの統計量を計算したところ,有効数字の関係上(?)上の出力結果とは小数点以下3桁以下が異なっていました.)

地理的加重回帰(GWR)モデルと通常最小二乗(OLS)モデルとの当てはまり具合の比較

  • 最後に,OLSとGWRとの比較を行い,二つのモデルの推定結果に統計的有意差があるかどうかを確かめましょう.ここでは,以下の参考文献1)の方法を紹介します.
  • 参考文献(いずれも,慶応大学の電子ジャーナルを使って閲覧可能です):

1) Leung, Y., Mei, C.-G., and Zhang, W.-X. 2000, Statistical tests for spatial non-stationarity based on the geographically weighted regression model, Environment and Planning A, 32, 9-32

2) Brunsdon, C., Fotheringham, A.S., and Charlton, M.E., 1999, Some notes on parametric signficance tests for geographically weighted regression, Journal of Regional Science, 39, 497-524

  • 文献1)には,GWRとOLSを比較する指標として,以下の3つの指標が示されています.詳細は,授業時に配布予定の参考文献を参照してください.

 F1値

  •  left
  • F1値が小さければOLSよりGWRの当てはまり具合がよいことを意味します.出力結果に含まれるp値(p-value)を用いて統計的な検定を行います.
  • F1値の算出方法
    > LMZ.F1GWR.test(col.gauss)
            Leung et al. (2000) F(1) test
    data:  col.gauss 
    F = 0.4928, df1 = 26.867, df2 = 46.000, p-value = 0.02634
    alternative hypothesis: less 
    sample estimates:
    SS OLS residuals SS GWR residuals 
            6014.856         1249.047

◆F2値

  •  left
  • F2値が小さければモデル間に統計的に有意な差がないことを意味します.出力結果に含まれるp値(p-value)を用いて統計的な検定を行います.
  • F2値の算出方法
    > LMZ.F2GWR.test(col.gauss)
            Leung et al. (2000) F(2) test
    data:  col.gauss 
    F = 1.3694, df1 = 33.39, df2 = 46.00, p-value = 0.1599
    alternative hypothesis: greater 
    sample estimates:
      SS OLS residuals SS GWR improvement 
              6014.856           4765.809 

◆F3値

  •  left
  • F3値は,モデルiのk番目の説明変数に対するパラメータ について,各モデルのパラメータ間に統計的有意差がないかどうかを検定するために用いられる指標です.
  • F3値が大きければパラメータ間に統計的な有意差があることを意味します.
> LMZ.F3GWR.test(col.gauss)
Leung et al. (2000) F(3) test
            F statistic Numerator d.f. Denominator d.f.  Pr(>)
(Intercept)     1.63378       16.81802            34.97 0.1083
income          1.35494       15.29487            34.97 0.2228
housing         0.95822        6.31335            34.97 0.4703
(参考)
  • 文献2)の方法には,以下のコマンドを使います.F値が大きければ,モデル間に統計的有意差があることを意味します.
> BFC99.gwr.test(col.gauss)
        Brunsdon, Fotheringham & Charlton (1999) ANOVA
data:  col.gauss 
F = 2.7787, df1 = 33.390, df2 = 26.867, p-value = 0.004027
alternative hypothesis: greater 
sample estimates:
SS GWR improvement   SS GWR residuals 
          4765.809           1249.047 

演習問題のためのTips:

Error in solve.default(inf, tol = tol.solve) : 
        system is computationally singular: reciprocal condition number = 2.31075e-020

などとエラーが出たら,tol.solveを2.31075e-020以下に設定してやるとよい. 以上で今回の演習は終了です.


トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2009-07-05 (日) 20:25:18 (3784d)