地震発生データと都道府県地価データをRであれこれして知見を得る(1)

最近機械学習に勉強ばっかりしてたところ

ふと

やっぱたまには原点に返って何かしてみることも重要だよ

というもう一人の自分からの声が聞こえたので、おもむろにデータから作ってみることに。

テーマ作ってデータ集めて
なんかするのは大変久々な気がする。やっぱ原点に返るのは（略）

主題

「地震発生回数が増大すると、翌年の地価に影響ってあんのかな？」
　　　→特に震度のでかさよ。

データをつくる

地震発生に関するデータはここから
地震情報 - Yahoo!天気・災害

都道府県別の地価データに関してはここの【表】2から
平成26年都道府県地価調査

それぞれお借りしました。

地価データは.xlsで取ってこれたのでまあいいとして、、、
問題は地震に関するデータ処理です。

うええ、、要前処理感、、、

RでWEBスクレイピングしたかったのですが苦手なので、
とりあえず昔作ったマクロでごそっと抜いてきました。

地価データが2013年、2014年のものなので
地震データはまず2013年のものを使うことにして、

行った前処理は以下です。

①震源地を国内に限定
②震源地を「都道府県」単位に加工
　　※湾の位置に無知すぎて死ぬほど時間かかった泣きそう。
　　※具体的には「浦河沖→北海道」とかのパターンを大量に作って
　　　機械的にやった感じですin Excel。
　　※このへんの処理はRでやったりexcelでやったり日頃はいろいろですが、、、
③これらをごそっと！ピボットテーブル！！！！
　　※反省はしてます。おかげで楽でした。
④地価データは遷移を見たいので2013年と2014年の差分を取りました。
④地価データと地震データをマージ

出来たデータがこちらです。

> data=read.table("clipboard",header=TRUE)
> data
       住宅地 商業地 工業地 総回数 震度１ 震度２ 震度３ 震度４ 震度５弱 震度５強 震度6弱
北海道      0    500  -5000    157     97     39     17      3        0        1       0
青森     -600  -1300   -600     32     18      8      5      1        0        0       0
岩手     -200      0   -400    121     83     23     10      5        0        0       0
宮城      900   6600    400    301    180     90     25      4        1        1       0
秋田     -600  -1400   -600     31     21      9      1      0        0        0       0
山形     -200   -700   -300     10      9      0      1      0        0        0       0
福島     1100    300    100    286    165     81     30      8        0        2       0
茨城     1400    500  -1700    317    182     97     24     10        4        0       0
栃木     -600  -2700   -400    131     79     39      9      3        0        1       0
群馬     -500  -1000  -1200     11      6      3      1      1        0        0       0
埼玉     1400   3300    200     10      7      2      1      0        0        0       0
千葉      800   4100   1300    138     90     33      9      6        0        0       0
東京     7300  67700   8500     37     25     10      2      0        0        0       0
神奈川   1900  21200   1300      9      7      2      0      0        0        0       0
新潟     -400  -1600   -300     36     19     15      1      1        0        0       0
富山      100    600   -200      2      1      1      0      0        0        0       0
石川     -200    700   -500      8      4      3      0      1        0        0       0
福井     -800  -1400   -100     10      8      2      0      0        0        0       0
山梨     -600  -1600  -1400      8      6      1      1      0        0        0       0
長野     -500  -1400   -600     57     34     17      6      0        0        0       0
岐阜     -500    400   -100     14      8      5      1      0        0        0       0
静岡     -500   -200   -800     20     12      6      0      2        0        0       0
愛知      400  31200  -1000      9      5      4      0      0        0        0       0
三重     -500   -300   -500      7      6      1      0      0        0        0       0
滋賀      100   1600      0      4      2      1      1      0        0        0       0
京都      300   8700   -200     13     10      3      0      0        0        0       0
大阪     -300  81400   -600     12      8      1      3      0        0        0       0
兵庫      700   4300   -900     24     15      6      2      0        0        0       1
奈良        0    300    300      7      4      2      1      0        0        0       0
和歌山   -800  -1700  -1400     54     32     14      4      3        0        0       0
鳥取     -700  -2100   -600      3      3      0      0      0        0        0       0
島根     -600  -1200   -400     10      9      1      0      0        0        0       0
岡山     -200    300   -100      2      2      0      0      0        0        0       0
広島     3900  18900  -2300     14      8      6      0      0        0        0       0
山口     -500  -1900  -1400      4      3      1      0      0        0        0       0
徳島     -700  -2400   -400     12      8      4      0      0        0        0       0
香川     -800  -1800   -800     14      8      6      0      0        0        0       0
愛媛     -900  -2500   -800     11      6      4      1      0        0        0       0
高知     -700  -2300   -600      7      4      3      0      0        0        0       0
福岡      200   3900    100      8      8      0      0      0        0        0       0
佐賀     -400  -1000   -200      5      3      2      0      0        0        0       0
長崎     -400   -600   -500     16     14      1      0      1        0        0       0
熊本      900   -700   -200     31     25      6      0      0        0        0       0
大分     -400  -1100   -500      7      6      1      0      0        0        0       0
宮崎     -200   -700   -300     31     21      6      4      0        0        0       0
鹿児島   -600  -1500  -5600     61     41     14      5      1        0        0       0
沖縄     1100   2400  -7300     75     56     10      7      2        0        0       0

2300件くらいあったのに、、こんなに少なくなってしまって、、、

ちなみに、東京都には伊豆諸島周辺の地震情報が含まれてしまい、
地価との関連から考えると伊豆諸島と東京一緒くたにするのは
どうもなーという感じだったので、
「東京島」というくくりにしていったん外出しにしてます。

何はともあれまず可視化

最近ggplot2の勉強してるけどなんかこう、、むずかしい

library(ggplot2)

#色のセッティング
cOrange="#f39800"
cLightBlue="#0068b7"
cGreen="#009944"

#住宅地×工業地の地価を眺める
data.plot=ggplot(data=data, aes(x=ID, y=data$"住宅地")) 
data.plot=data.plot+geom_point(size=3,colour=cOrange)
data.plot=data.plot+geom_point(data=data,aes(x=ID,y=data$"工業地"),colour=cLightBlue,size=3)
data.plot+ylab("前年比地価変動額")

#住宅地×商業地の地価を眺める
data.plot2=ggplot(data=data, aes(x=ID, y=data$"住宅地")) 
data.plot2=data.plot2+geom_point(size=3,colour=cOrange)
data.plot2=data.plot2+geom_point(data=data,aes(x=ID,y=data$"商業地"),colour=cGreen,size=3)
data.plot2+ylab("前年比地価変動額")

自信無さ過ぎてレイヤー重ねるごとにステップ確認してるから糞コードすぎて見たくない、見えない。
なんでプロット重ねたりなんてしたんだろう、、、

とりあえず住宅地の地価変動を基準に、「住宅地と工業地」「住宅地と商業地」で２つ、
都道府県別に並べて作っています。

凡例つけようが無い作り方をしてしまったので

住宅地：オレンジっぽいやつ
工業地：青っぽいやつ
商業地：緑っぽいやつ

チイサクテミエネエ！！！

何はともあれ。

地震によってほんのりとでも地価に影響するとすれば。
住宅地の地価が一番影響受けやすいんじゃないのかなーというなんとなくな仮説を持っています。

住宅地（オレンジ）と比べると
工業地（青）が一部マイナス方向に触れたりしていて
商業地（緑）は住宅地と同じ動きをしつつも振れ幅がでかい感じ。特に東京近郊。

外れ値にしちゃいたくなるような状況ですが
データ数も元々少ないしあえて見て見ぬふりをします。

今考えたらggmapがあるんだから
地図上にマッピングした方が見やすいよねえ、、後でやろう、、、

マルチコ対応

さっそく重回帰、、、と思ったら、
なんにも考えてなかったマルチコ対応。

とりあえずVIF値を見てみる

>library(DAAG)
>model=lm(y1~x0+x1+x2+x3+x4+x5.1+x5.2+x6,data=data)
> vif(model)
        x0         x1         x2         x3         x4       x5.1       x5.2         x6 
320920.000 110940.000  27791.000   2680.700    296.440     20.644     12.705      2.148 

X0:地震総回数
x1～x6：各震度の発生回数

･･･うん、ばかだった。そりゃあx0とその他を同時にぶちこんじゃだめだった

ちなみにこのVIF値、
「10以上はアウトって言われてるけど10以下だったら大丈夫なんて保証は無いからな」
ってよくやり玉にあがるくらいなので、今回みたいなアホケースはさておき、
マルチコ問題は根が深いなあっていつも思う。

念のため散布図行列で俯瞰してみる

なんか震度1×震度2もふんわりと怪しくない？
余震とか絡めたら複雑なんだろうなあ、、そっちの分野の方に殴られてしかるべき感じ。

#x0を除外してやり直し
> model=lm(y1~x1+x2+x3+x4+x5.1+x5.2+x6,data=data)
> vif(model)
     x1      x2      x3      x4    x5.1    x5.2      x6 
58.8270 54.9970 30.1050  5.9521  6.5606  7.8988  1.0148 

#看過出来ないので震度3以上の発生分布のみを考える
> vif(model2)
     x3      x4    x5.1    x5.2      x6 
11.4120  4.9629  2.7962  5.0712  1.0136 

#やっぱ同じ震源で地震が頻発しまくる、みたいなのはよくありそう
#とりあえずこれで重回帰


#住宅地地価変動～震度3+震度4+震度5弱+震度5強+震度6弱

> summary(model2)

Call:
lm(formula = y1 ~ x3 + x4 + x5.1 + x5.2 + x6, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1010.2  -675.2  -398.7    59.1  7120.0 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    40.29     256.86   0.157    0.876
x3             69.86     103.55   0.675    0.504
x4           -111.41     214.60  -0.519    0.606
x5.1          184.66     594.70   0.311    0.758
x5.2         -392.03    1277.90  -0.307    0.761
x6            519.98    1468.71   0.354    0.725

Residual standard error: 1443 on 41 degrees of freedom
Multiple R-squared:  0.04124,	Adjusted R-squared:  -0.07568 
F-statistic: 0.3528 on 5 and 41 DF,  p-value: 0.8775

すがすがしいほどに棄却されなかった。笑

なんか頑張ってどうにかなる感じじゃない。
震度4以上でもだめだったのでまあ、、
そもそもデータ数が増えればどんどんマルチコしていくタイプの
データな気がしてきた、今更。

仮説設定自体に問題があったというアレ。

諦めきれず総回数で単回帰もやってみた

> model4=lm(y1~x0,data=data)
> summary(model4)

Call:
lm(formula = y1 ~ x0, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1043.2  -668.3  -430.0   153.5  7170.1 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    6.650    236.417   0.028    0.978
x0             3.332      2.650   1.257    0.215

Residual standard error: 1383 on 45 degrees of freedom
Multiple R-squared:  0.03394,	Adjusted R-squared:  0.01247 
F-statistic: 1.581 on 1 and 45 DF,  p-value: 0.2151

ふふ。だめだこりゃ。

他にも色々やってみたかったけど

単年度での地価変動に地震発生回数が影響を及ぼすかは
今回範囲では妄言に過ぎないといえそうです

久々に重回帰やったらポンコツな出来だった反省。

せっかくデータは取ったから地図上で可視化とかやろうかなあ

なーんかおもしろいひらめきないかなー

3月3日・うまくいかなかったモデルを可視化してみる

結局そもそも検討違いなモデル組をした気がする今回のテーマですが、
せっかくなので作った重回帰モデルを可視化。反省会かねて。

> par(mfrow=c(2,2)) #４つの図を並べる
> plot(model2)

うへえ、、、

左上：残差VSフィット値
　　　なんかやっぱまとまり無いなあというか、そもそも震度1，2のデータが大量、かつ
　　　震度3以上のデータが貧弱というデータ構造なのが分かりやすいというか、、、

右上：残差の正規
　　　直線状にplotが並べば正規分布に従ってる･･･て解釈をいつもしてるんですが、
　　　今回どうなんですかね。正規分布ですかね。

左下：残差変動状況
　　　フィット値plotとの感想が同じです。データ構造がなんか

右下：残差と影響力（クック距離）
　　　点線のクック距離が0.5を超えると影響力大きいよねーみたいな解釈なんですが、
　　　震度5,6らへんになってくるとそもそも発生回数自体が少ないので、
　　　それに伴って一件あたりの影響力がでかくなるのは当然よね！っていう。
　　　しかしこういうラインを描くクック距離初めて見た気がしなくもないです

以上。

そこらへんにあるデータをなんとなく引っ張ってきて始めた今回ですが
地価はともかくとして
地震って取扱い難しいね、、、どう扱っていいやらおっかなびっくり。

なんかまた閃いたらやる。

時系列で見るとか･･･
震度5以上の発生回数推移で考えるとか、、
地価の上下だけにわけてSVMとか、、、

まあ、、いろいろ。