QGIS, Rを用いた公示地価データの空間統計分析

ここでは、QGIS, Rを用いて空間統計学の基礎的な演習を説明しています。

    参考文献:
  1. 瀬谷創・堤盛人、空間統計学-自然科学から人文・社会科学まで-、朝倉書店、2014
  2. 古谷知之、Rによる空間データの統計分析、朝倉書店、2011
    演習内容:
  1. GRASSのインストール、標高データを用いた地滑り危険度マップの作成
  2. 反射率、輝度温度、標高データを用いた土地被覆分類
  3. QGIS, Rを用いた公示地価データの空間統計分析
    本演習内容:
  1. Rのインストール
  2. 使用するデータ
  3. データの前処理
    1. (市区町村ポリゴン)関西二府四県の抽出、面積計算 (QGIS)
    2. (公示地価データ)緯度、経度座標の出力 (QGIS)
    3. (公示地価データ)特定の地点までの距離の計算 (QGIS)
    4. (人口データ)昼夜間の人口密度の計算 (Excel)
    5. (公示地価データ)人口密度データの付与 (Excel)
  4. 空間的自己相関指標の計算
  5. 重回帰モデルのパラメータ推定
    1. (通常の)回帰モデル
    2. 同時自己回帰モデル (simultaneous autoregressive model: SAR)
    3. 条件付き自己回帰モデル (conditional autoregressive model: CAR)
    4. 空間ラグモデル (spatial lag model: SLM)
    5. 誤差項の空間的自己回帰モデル (spatial error model: SEM)
  6. 演習課題

Rのインストール

以下のサイトからインストーラーをダウンロードし、インストールする。

使用するデータ

データの前処理

空間統計解析の前に、上記のデータを処理する必要がある。

(市町村ポリゴン)関西二府四県の抽出、面積計算 (QGIS)

QGISを起動し、「japan_ver80.shp」を表示させる。その上で、下記の処理で関西二府四県の市町村を選択する。

「レイヤパネル(QGISの左下の窓)」中の「japan_ver80.shp」で右クリック
→ 「属性テーブルを開く」
→(属性テーブルウィンドウの左下のボタンから)「アドバンストフィルター」

「アドバンストフィルター」内で記述する式

"KEN" = '京都府' OR "KEN" = '大阪府' OR "KEN" = '兵庫県' OR "KEN" = '滋賀県' OR "KEN" = '奈良県' OR "KEN" = '和歌山県'

その後、選択した市町村を別のshapeファイルとして保存する。

「レイヤパネル(QGISの左下の窓)」中の「japan_ver80.shp」で右クリックし「名前を付けて保存」
「アドバンストフィルター」で選択された関西二府四県の市区町村 別名ファイルでの保存後の関西二府四県の市区町村
選択された関西二府四県の市区町村 抽出後の関西二府四県の市区町村

この後、各市町村の面積を計算する。但し、「japan_ver80_Kansai.shp」の「空間参照システム」が「JGD2000」であるため、そのまま計算すると面積が度の単位で返されてしまう。km2の単位での面積を求めるため、下記の処理を行う。

  「プロジェクト」
→ 「プロジェクトプロパティ」
→ 「'オンザフライ'CRS変換を有効にする」にチェックを入れる。
→ 「フィルター」に「JGD」というキーワードを入れて、
「JGD2000 / Japan Plane Rectangular CS VI」を選択する。

すると、平面直角座標系第6系で表示されるようになる。その後、 属性テーブルを開き、下記の処理を行う。

(属性テーブルウィンドウの左上の)「編集モード切替」ボタン(鉛筆のアイコン)をクリック
→ 「フィールド計算機を開く」(メニューバーの右寄りのアイコン)
空間参照システムの変更 フィールド計算機での計算式
空間参照システムの変更 フィールド計算機での計算式
km2単位の面積
km2単位の面積

(公示地価データ)緯度、経度座標の出力 (QGIS)

公示地価データ「L01-10_26-g_LandPrice.shp」の「空間参照システム」は「JGD2000」である。上記の処理と同様に、一度別ファイルとして、平面直角座標系第6系が指定されたshapeファイルを保存する。(例えば、「L01-10_26-g_LandPrice_VI.shp」)

その後、属性テーブルを開き、下記の処理をフィールド名「X」と「Y」(共にm単位)について行う。

(属性テーブルウィンドウの左上の)「編集モード切替」ボタン(鉛筆のアイコン)をクリック
→ 「フィールド計算機を開く」(メニューバーの右寄りのアイコン)

(公示地価データ)特定の地点までの距離の計算 (QGIS)

後の重回帰分析において特定地点との直線距離(本来なら移動に要する時間距離が望ましい)を用いる為、ここで計算する。以下、四条駅、大阪駅の平面直角座標系XI系での座標値(単位はm)を記す。これを利用して、上記で作成した「L01-10_26-g_LandPrice_VI.shp」の属性テーブル内に新たなフィールドを設け、距離を格納する。

(属性テーブルウィンドウの左上の)「編集モード切替」ボタン(鉛筆のアイコン)をクリック
→ 「フィールド計算機を開く」(メニューバーの右寄りのアイコン)

(人口データ)昼夜間の人口密度の計算 (Excel)

まず、「昼夜間人口データ.xls」をExcelで開き、「平成22年」のワークシートを開く。各班で決定する課題においては種々の人口データを活用してほしいが、ここでは夜間人口と昼間人口に着目する。

一方、「japan_ver80_Kansai.shp」の属性テーブル「japan_ver80_Kansai.dbf」もExcelで開く。表1のように市区町村名で昇順に並び替えておく。この表を「昼夜間人口データ.xls」の新たなワークシート(例えば「面積」)に貼り付けておく。

表1:市区町村名で昇順に並び替えた「japan_ver80_Kansai.dbf」のデータの一部

SIKUCHOSON CITY_ENG	P_NUM(人口)	H_NUM(世帯数)	Area
綾部市	Ayabe-shi	35419	15832	348.6 
伊根町	Ine-cho	2339	941	61.2 
井手町	Ide-cho	7899	3375	18.2 
右京区	Kyoto-shi, Ukyo-ku	195921	92735	293.4 
宇治市	Uji-shi	190856	81882	67.9 
宇治田原町	Ujitawara-cho	9735	3572	58.5 
下京区	Kyoto-shi, Shimogyo-ku	76814	42940	6.8 
笠置町	Kasagi-cho	1529	675	23.6 
亀岡市	Kameoka-shi	91548	37993	225.8 
久御山町	Kumiyama-cho	16390	6809	13.9 
宮津市	Miyazu-shi	19452	8694	174.2 
京丹後市	Kyotango-shi	58514	22678	504.6 
京丹波町	Kyotamba-cho	15555	6405	304.5 
京田辺市	Kyotanabe-shi	66697	27001	43.1 
向日市	Muko-shi	54340	23226	7.8 
左京区	Kyoto-shi, Sakyo-ku	156069	77067	248.0 
山科区	Kyoto-shi, Yamashina-ku	131699	63369	28.8 
上京区	Kyoto-shi, Kamigyo-ku	76947	40677	7.1 
城陽市	Joyo-shi	78560	34328	32.8 
精華町	Seika-cho	37443	14269	25.9 
西京区	Kyoto-shi, Nishikyo-ku	151827	66038	59.6 
大山崎町	Oyamazaki-cho	15468	6427	6.1 
中京区	Kyoto-shi, Nakagyo-ku	104345	55138	7.5 
長岡京市	Nagaokakyo-shi	80222	34842	19.3 
東山区	Kyoto-shi, Higashiyama-ku 37183	20688	7.5 
南区	Kyoto-shi, Minami-ku	98433	48872	15.9 
南山城村	Minamiyamashiro-mura	2961 1249 64.4 
南丹市	Nantan-shi	33418	13958	619.2 
八幡市	Yawata-shi	73038	31944	24.3 
舞鶴市	Maizuru-shi	86996	40796	343.6 
伏見区	Kyoto-shi, Fushimi-ku	278230	132889	61.9 
福知山市	Fukuchiyama-shi	80682	35689	555.1 
北区	Kyoto-shi, Kita-ku	112006	53017	95.2 
木津川市	Kizugawa-shi	73095	27560	85.4 
与謝野町	Yosano-cho	23318	9158	108.8 
和束町	Wazuka-cho	4357	1751	65.3 

ここで、「昼夜間人口データ.xls」の新たな列に、各市区町村の面積を参照するものとする。Excelでよく用いられるvlookupという関数を利用する。vlookupやhlookupは便利である反面、指定通りにデータを準備しないとエラーが返されるだけなので注意を要する。

(京都市)「北区」の面積を参照する関数

=VLOOKUP(B1224,面積!$A$1:$E$36,5)

vlookup時の注意点

夜間人口、昼間人口を上記の面積データで割って、人口密度データを計算する。下記の表2に一例を示す(関係のない列は割愛している)。

表2:市区町村名で昇順に並び替えた「昼夜間人口データ.xls」のデータの一部
(面積:km2、夜間人口密度・昼間人口密度:人/km2)

地域 夜間人口 昼間人口 面積 夜間人口密度 昼間人口密度
北区	122,037	129,464	95.2 	1281.4 	1359.3 
上京区	83,264	97,833	7.1 	11781.7 13843.2 
左京区	168,802	176,443	248.0 	680.7 	711.5 
中京区	105,306	155,123	7.5 	14123.2 20804.4 
東山区	40,528	56,253	7.5 	5389.5 	7480.6 
下京区	79,287	135,656	6.8 	11653.5 19938.6 
南区	98,744	136,137	15.9 	6208.6 	8559.8 
右京区	202,943	189,617	293.4 	691.7 	646.3 
伏見区	284,085	276,715	61.9 	4585.8 	4466.8 
山科区	136,045	121,223	28.8 	4726.3 	4211.4 

(公示地価データ)人口密度データの付与 (Excel)

公示地価データの属性テーブル「L01-10_26-g_LandPrice.dbf」をExcelで開く。このデータに上記で生成した、昼夜間人口密度データを付与していきたい。

各地点が属する市区町村名は「L01_018」のフィールドに記載されている。見ればわかるように、「京丹後」「福知山」のように、「〜市」「〜町」「〜区」といった表示が割愛されている。また「京都市北区」「京都市南区」は「京都北」「京都南」と記述されている。今回も「vlookup」を用いて該当するデータを参照していく処理を考える。ただ参照元と参照先のデータ名を一致させておく必要がある。今回は参照先である「昼夜間人口データ.xls」の表記を変更するとする。

京都府のデータに限ると高々35個のデータなので、手動でもそれほど時間を要しない。ただ筆者が実施した処理を下記に紹介するので、参考にしてほしい。

参考:秀丸での検索、置換(演習問題)

  1. 「昼夜間人口データ.xls」の京都府内の地域名(例えば、「北区」「上京区」の列)をコピーして、テキストエディタ「秀丸」の新規ファイルに貼り付ける。
  2. まず、「市」を除去する。滋賀県の「八日市市」のように市名に「市」が含まれる場合もある。末尾の「市」だけ除去するために、「秀丸」の「置換」の機能を使う。
    • 「正規表現」をチェックする
    • 検索:市\n (「\n」は改行マーク)
    • 置換:\n
  3. 同様に、「区」「町」「村」についても上記の処理を繰り返す。(「市」「区」「町」「村」の順番でなくても不具合はないはず)
  4. 「北」「南」については、手動で「京都北」「京都南」に置換する。
  5. その他に、公示地価データの属性テーブル「L01-10_26-g_LandPrice.dbf」と一致していない表記がないか確認する。
  6. 確認後に、コピーして、「昼夜間人口データ.xls」の地域名の列に貼り付ける。

この後、公示地価データの属性テーブル「L01-10_26-g_LandPrice.dbf」に「vlookup」を使って、昼夜間の人口密度データを付与する。vlookupの使い方は上記の説明と同じなので割愛する。

最後に、Excel上で今回の演習に関係のない列を削除して、表3に示すように列の項目名を変更する。その上で、CSVファイル「Kyoto10_Price.csv」として保存する。

表3:今回の演習に関係のない列を削除した最終的な「Kyoto10_Price.csv」(一部)

Price Name Station Dist2Station LU_Type X Y Dist2Shijo Dist2Osaka Pop_Day Pop_Night
20500 京丹後 網野 2800 都市 -88727.760400 -34220.500800 101389.614 117519.498 117.0063648 113.855172
35500 京丹後 網野 700 都市 -88036.598900 -35927.693800 99651.729 115677.869 117.0063648 113.855172
27700 京丹後 網野 2000 都市 -87853.618600 -34644.832400 100495.227 116808.790 117.0063648 113.855172
31500 京丹後 峰山 1400 都市 -85322.116500 -41517.504500 93683.096 109485.344 117.0063648 113.855172
42900 京丹後 峰山 500 都市 -84414.309900 -40967.711300 93479.503 109678.295 117.0063648 113.855172
25000 京丹後 丹後大宮 850 都計外 -81716.056600 -46111.195400 87861.835 103918.760 117.0063648 113.855172

空間的自己相関指標の計算

Rを起動し、「Kyoto10_Price.csv」が存在するディレクトリへ移動する。

メニューバーの「ファイル」→「ディレクトリの変更」で指定

以下、Rでの実行画面を示す。赤字が入力文字を表す。

注意: spdepというライブラリを読み込んでいるが、一連の処理の最初に一度だけ実行すればいい。もしエラーが返されれば、「メニュー」→「パッケージ」からダウンロードして入手する。

> library(spdep)
> price <- read.table("Kyoto10_Price.csv", sep=",", header=TRUE)
> summary(price)
     Price              Name         Station     Dist2Station  
 Min.   :   1150   伏見   : 60   北大路  : 24   Min.   :    0  
 1st Qu.:  87375   左京   : 51   福知山  : 21   1st Qu.:  500  
 Median : 147000   宇治   : 50   桂      : 17   Median :  900  
 Mean   : 180155   右京   : 43   大久保  : 16   Mean   : 1321  
 3rd Qu.: 218250   山科   : 36   長岡天神: 15   3rd Qu.: 1600  
 Max.   :3260000   京都北 : 33   亀岡    : 14   Max.   :13000  
                   (Other):407   (Other) :573                  
       LU_Type          X                Y             Dist2Shijo      
 1低専     :152   Min.   :-88728   Min.   :-141494   Min.   :   256.3  
 1住居     : 87   1st Qu.:-27913   1st Qu.:-119458   1st Qu.:  4643.1  
 1住居,準防: 66   Median :-22850   Median :-112067   Median :  8369.8  
 近商,準防 : 48   Mean   :-29012   Mean   :-109713   Mean   : 16753.7  
 調区      : 39   3rd Qu.:-19985   3rd Qu.:-106960   3rd Qu.: 18214.8  
 1中専     : 38   Max.   :-11800   Max.   : -34221   Max.   :101389.6  
 (Other)   :250                                                        
   Dist2Osaka        Pop_Day           Pop_Night       
 Min.   : 23617   Min.   :   51.67   Min.   :   47.58  
 1st Qu.: 33620   1st Qu.:  680.71   1st Qu.:  647.55  
 Median : 38864   Median : 2567.84   Median : 2091.09  
 Mean   : 42718   Mean   : 3402.11   Mean   : 4043.03  
 3rd Qu.: 43643   3rd Qu.: 4585.75   3rd Qu.: 4466.78  
 Max.   :117520   Max.   :14123.18   Max.   :20804.42  

> coords <- matrix(0, nrow=nrow(price), ncol=2)
> coords[,1] <- price$X
> coords[,2] <- price$Y
> price.tri.nb <- tri2nb(coords)

# Moran's I test
> moran.test(price$Price, nb2listw(price.tri.nb, style="W"))

        Moran I test under randomisation

data:  price$Price  
weights: nb2listw(price.tri.nb, style = "W")  

Moran I statistic standard deviate = 22.865, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.4737820920     -0.0014727541      0.0004320206 

# Geary's C test
> geary.test(price$Price, nb2listw(price.tri.nb, style="W"))

        Geary C test under randomisation

data:  price$Price 
weights: nb2listw(price.tri.nb, style = "W") 

Geary C statistic standard deviate = 9.1825, p-value < 2.2e-16
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic       Expectation          Variance 
      0.492007962       1.000000000       0.003060513 

重回帰モデルのパラメータ推定

(通常の)回帰モデル

以下の手順に従って入力していく。

注意: 推定した結果はまとめて、「lm_result」というオブジェクトに保存されているが、個別の要素は「names()」を使って確認できる。Rでも図示できるようだが、筆者が確認した限りでは十分に鮮明な結果が得られなかった。そのため、本演習では推定結果の一部をファイルに出力し、QGIS上で図示するものとする。

> library(spdep)
> price <- read.table("Kyoto10_Price.csv", sep=",", header=TRUE)

> lm_result <- lm(Price~Dist2Shijo+X+Y+Pop_Night, data=price)

> summary(lm_result)

Call:
lm(formula = Price ~ Dist2Shijo + X + Y + Pop_Night, data = price)

Residuals:
    Min      1Q  Median      3Q     Max 
-300994  -55912  -11388   32065 2784660 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.962e+05  9.116e+04   5.443 7.34e-08 ***
Dist2Shijo  -3.166e+00  7.200e-01  -4.397 1.28e-05 ***
X            9.232e-01  1.085e+00   0.851    0.395    
Y            2.708e+00  6.334e-01   4.276 2.18e-05 ***
Pop_Night    1.507e+01  1.305e+00  11.548  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 167600 on 675 degrees of freedom
Multiple R-squared:  0.332,     Adjusted R-squared:  0.328 
F-statistic: 83.86 on 4 and 675 DF,  p-value: < 2.2e-16

> names(lm_result)
 [1] "coefficients" "residuals" "effects"  "rank"
 [5] "fitted.values" "assign" "qr" "df.residual"
 [9] "xlevels" "call" "terms" "model"

> out.file <- file("Kyoto10_Price_rst.csv", open="w")

> x_str <- sprintf("%f,%f,%f,%f", price$X, 
 price$Y, price$Price, lm_result$residuals)(本当は1行で記述)

> writeLines(x_str, out.file)
> close(out.file)

推定結果をQGIS上で表示する。まず、上記の出力結果をGIS上に取り込む。

  メニューバーの「レイヤ」の「レイヤの追加」
→「デリミティッドテキストレイヤの追加」

「デリミティッドテキストレイヤの追加」 「デリミティッドテキストレイヤの追加」の様子

次に、他のデータと重ね合わせて表示させるためにも、上記のファイルに座標系を与えて保存する。

「レイヤパネル」中の「Kyoto10_Price_rst」で右クリックし、「名前を付けて保存」

この時点で表示されている「Kyoto10_Price_rst」は表示画面から削除する。最終結果を分かりやすく表示させるために、下記のデータも事前にダウンロードし表示させておく。(注意:プロジェクトのCRSは「JGD2000」とし、「'オンザフライ'CRS変換を有効にする」にチェックを入れる

最後に、上記で作成した「Kyoto10_Price_rst.shp」を表示させ、残差に応じて表示させる色を変えていく。

  「レイヤパネル」中の「Kyoto10_Price_rst」で右クリック
→「プロパティ」
→「スタイル」

線形回帰モデルの推定結果 線形回帰モデルの推定結果。正負共に残差が大きい地点が多数散見され、推定に使用した線形モデルが不適切であると予想できる。

同時自己回帰モデル (simultaneous autoregressive model: SAR)

SARモデルパラメータ推定後の残差は、「SAR_result$fit$residuals」に格納されている点に注意する。

> sar_result <- spautolm(Price~Dist2Shijo+X+Y+Pop_Night, 
  data=price, nb2listw(price.tri.nb, style="W"), 
  family="SAR")
> summary(sar_result)

Call: spautolm(formula = Price ~ Dist2Shijo + X + Y + Pop_Night, 
  data = price, listw = nb2listw(price.tri.nb, style = "W"), 
  family = "SAR")

Residuals:
      Min        1Q    Median        3Q       Max 
-368324.8  -39917.1   -7905.6   19640.0 2727348.8 

Coefficients: 
               Estimate  Std. Error z value  Pr(>|z|)
(Intercept)  4.8938e+05  1.4055e+05  3.4820 0.0004977
Dist2Shijo  -3.3432e+00  1.1068e+00 -3.0206 0.0025225
X            7.0138e-01  1.6624e+00  0.4219 0.6731003
Y            2.6339e+00  9.8222e-01  2.6816 0.0073277
Pop_Night    1.4370e+01  1.8322e+00  7.8427 4.441e-15

Lambda: 0.42335 LR test value: 72.818 p-value: < 2.22e-16 
Numerical Hessian standard error of lambda: 0.046141 

Log likelihood: -9105.795 
ML residual variance (sigma squared): 2.4202e+10, (sigma: 155570)
Number of observations: 680 
Number of parameters estimated: 7 
AIC: 18226

> out.file <- file("Kyoto10_Price_SAR_rst.csv", open="w")

> x_str <- sprintf("%f,%f,%f,%f", price$X, 
 price$Y, price$Price, sar_result$fit$residuals)

> writeLines(x_str, out.file)
> close(out.file)

通常の回帰モデル推定後の処理と同様に、QGIS上で下記の処理を行うが、誤差分布の表示において注意を要する。

  1. デリミティッドテキストレイヤの追加
  2. 「JGD2000 / Japan Plane Rectangular CS VI」のCRSを持つshapeファイルとして保存
  3. 「スタイル」を変更し、残差に応じて色を変更
    注意:自動的に各階級に対応する上限値、下限値が決められる。他の結果と比較するために、通常の回帰モデルで与えられた上限値、下限値を手動で入力する。

条件付き自己回帰モデル (conditional autoregressive model: CAR)

> car_result <- spautolm(Price~Dist2Shijo+X+Y+Pop_Night, data=
  price, nb2listw(price.tri.nb, style="W"), family="CAR")
> summary(car_result)

Call: spautolm(formula = Price ~ Dist2Shijo + X + Y + Pop_Night, 
  data = price, listw = nb2listw(price.tri.nb, style = "W"), 
  family = "CAR")

Residuals:
      Min        1Q    Median        3Q       Max 
-747665.1  -25729.3   -3382.2   15782.7 2663876.8 

Coefficients: 
               Estimate  Std. Error z value  Pr(>|z|)
(Intercept) 487846.0530 172832.4135  2.8227  0.004763
Dist2Shijo      -3.4955      1.3613 -2.5677  0.010237
X                0.5323      2.0235  0.2631  0.792508
Y                2.5498      1.2155  2.0976  0.035937
Pop_Night       14.4365      2.0666  6.9856 2.836e-12

Lambda: 0.81798 LR test value: 88.678 p-value: < 2.22e-16 
Numerical Hessian standard error of lambda: 0.056908 

Log likelihood: -9097.865 
ML residual variance (sigma squared): 2.2514e+10, (sigma: 150050)
Number of observations: 680 
Number of parameters estimated: 7 
AIC: 18210

> out.file <- file("Kyoto10_Price_CAR_rst.csv", open="w")

> x_str <- sprintf("%f,%f,%f,%f", price$X, 
 price$Y, price$Price, car_result$fit$residuals)

> writeLines(x_str, out.file)
> close(out.file)

上記のSARモデル推定後と同様に処理する。

空間ラグモデル (spatial lag model: SLM)

> slm_result <- lagsarlm(Price~Dist2Shijo+X+Y+Pop_Night, 
  data=price, nb2listw(price.tri.nb, style="W"))
> summary(lag_result)

Call:lagsarlm(formula = Price ~ Dist2Shijo + X + Y + Pop_Night, 
  data = price, listw = nb2listw(price.tri.nb, style = "W"))

Residuals:
      Min        1Q    Median        3Q       Max 
-385059.3  -37767.1   -7370.8   20333.0 2736621.4 

Type: lag 
Coefficients: (numerical Hessian approximate standard errors) 
               Estimate  Std. Error z value  Pr(>|z|)
(Intercept)  2.7271e+05  8.8270e+04  3.0895  0.002005
Dist2Shijo  -1.7003e+00  6.9302e-01 -2.4534  0.014150
X            5.2089e-01  1.0165e+00  0.5125  0.608330
Y            1.4850e+00  6.0401e-01  2.4587  0.013946
Pop_Night    9.6840e+00  1.3386e+00  7.2344 4.676e-13

Rho: 0.42173, LR test value: 76.224, p-value: < 2.22e-16
Approximate (numerical Hessian) standard error: 0.044985
    z-value: 9.3749, p-value: < 2.22e-16
Wald statistic: 87.889, p-value: < 2.22e-16

Log likelihood: -9104.092 for lag model
ML residual variance (sigma squared): 2.4088e+10, (sigma: 155200)
Number of observations: 680 
Number of parameters estimated: 7 
AIC: 18222, (AIC for lm: 18296)

> out.file <- file("Kyoto10_Price_SLM_rst.csv", open="w")

> x_str <- sprintf("%f,%f,%f,%f", price$X, 
 price$Y, price$Price, slm_result$residuals)

> writeLines(x_str, out.file)
> close(out.file)

誤差項の空間的自己回帰モデル (spatial error model: SEM)

> sem_result <- errorsarm(Price~Dist2Shijo+X+Y+Pop_Night, 
  data=price, nb2listw(price.tri.nb, style="W"))
> summary(err_result)

Call:errorsarlm(formula = Price ~ Dist2Shijo + X + Y + Pop_Night, 
  data = price, listw = nb2listw(price.tri.nb, style = "W"))

Residuals:
      Min        1Q    Median        3Q       Max 
-368324.8  -39917.1   -7905.6   19640.0 2727348.8 

Type: error 
Coefficients: (asymptotic standard errors) 
               Estimate  Std. Error z value  Pr(>|z|)
(Intercept)  4.8938e+05  1.4055e+05  3.4820 0.0004977
Dist2Shijo  -3.3432e+00  1.1068e+00 -3.0206 0.0025225
X            7.0138e-01  1.6624e+00  0.4219 0.6731003
Y            2.6339e+00  9.8222e-01  2.6816 0.0073277
Pop_Night    1.4370e+01  1.8322e+00  7.8427 4.441e-15

Lambda: 0.42335, LR test value: 72.818, p-value: < 2.22e-16
Approximate (numerical Hessian) standard error: 0.046149
    z-value: 9.1735, p-value: < 2.22e-16
Wald statistic: 84.153, p-value: < 2.22e-16

Log likelihood: -9105.795 for error model
ML residual variance (sigma squared): 2.4202e+10, (sigma: 155570)
Number of observations: 680 
Number of parameters estimated: 7 
AIC: 18226, (AIC for lm: 18296)

> out.file <- file("Kyoto10_Price_SEM_rst.csv", open="w")

> x_str <- sprintf("%f,%f,%f,%f", price$X, 
 price$Y, price$Price, sem_result$residuals)

> writeLines(x_str, out.file)
> close(out.file)
線形回帰モデルの推定結果 SARモデルの推定結果
線形回帰モデルの推定結果 同時自己回帰モデル(SAR)の推定結果
CARモデルの推定結果 SLMモデルの推定結果
条件付き自己回帰モデル
(CAR)の推定結果
空間的自己回帰モデル(SLM)の推定結果
SEMモデルの推定結果
誤差項の空間的自己回帰モデル
(SEM)の推定結果

演習課題

指定日に各班の代表が発表します。それまでに各班で共同で作業をするように。

  1. 京都の中でも京都市のデータに限定して、上記の空間統計学の処理(空間的自己相関指標の計算、重回帰モデルのパラメータ推定)を実施し、結果を示せ。
  2. 京都以外の複数の府県のデータを統合して処理した場合(例えば、京都府、大阪府、滋賀県のデータを統合)の結果を示せ。
  3. 他の国土数値情報を説明変数に取り込んで、回帰モデルのパラメータを推定せよ。
  4. 地価を説明する回帰モデルに関して、その他の班独自の解析を試みよ。

[GRASSのインストール、標高データを用いた地滑り危険度マップの作成]
[反射率、輝度温度、標高データを用いた土地被覆分類]
[QGIS, Rを用いた公示地価データの空間統計分析]
[須崎純一のページへ]
須崎純一 京都大学大学院 工学研究科社会基盤工学専攻 空間情報学講座