2016年4月4日 月曜日 曇り
データ解析環境「Ｒ」
参考文献 船尾・高浪 データ解析環境「Ｒ」 工学社 2005年
（以下の表示では行列の要素間のスペースが上手く並んでいません。ワードプレスの編集画面では綺麗に並んでいるのですが、ブラウザ画面では上手くいきません。見苦しい状況、どうかご容赦ねがいます。)
＊＊＊＊＊
ためしに今年の三笠市の降雪や積雪量のデータを見ながら整理してみた。2016年2月の1日から10日までの日付をDATE=day、それぞれの一日の降雪量（cm）をDAILYSNOW=daylysnow、各日の積雪量をSNOWDEPTH=snowdepth（cm）、シーズン中の累積降雪量をYEARLYSNOW=yearsnow（cm）、とおいてみた。
day < - 1:10
dailysnow <- c(18,7,18,2,3,55,15,0,6,2)
snowdepth <- c(128, 120, 126,115,110,149,137,125,120,113)
yearsnow <- c(634,641,659,661,664,719,734,734,740,742)
#2月になってからの累積降雪量FEB_SNOW=feb_snowは、cumsum(dailysnow) である。（参考書、p82の関数表より）
feb_snow <- cumsum(dailysnow)
feb_snow
#feb_snow = yearsnow - 616 + cumsum(dailysnow) としても同じことである。
measurement <- data.frame(DATE=day, DAILYSNOW=dailysnow, SNOWDEPTH=snowdepth, YEARLYSNOW=yearsnow, FEB_SNOW=feb_snow)
measurement
＊＊＊＊＊
以上を「Ｒ」で計算させてみると以下のようになる。
> day < - 1:10
> dailysnow < - c(18,7,18,2,3,55,15,0,6,2)
> snowdepth < - c(128, 120, 126,115,110,149,137,125,120,113)
> yearsnow < - c(634,641,659,661,664,719,734,734,740,742)
>
> #2月になってからの累積降雪量FEB_SNOW=feb_snowは、cumsum(dailysnow) である。（参考書、p82の関数表より）
>
> feb_snow < - cumsum(dailysnow)
> feb_snow
18 25 43 45 48 103 118 118 124 126
>
> #feb_snow = yearsnow – 616 + cumsum(dailysnow) としても同じことである。
>
> measurement < - data.frame(DATE=day, DAILYSNOW=dailysnow, SNOWDEPTH=snowdepth, YEARLYSNOW=yearsnow, FEB_SNOW=feb_snow)
#データ・フレームmeasurementの中身を確認すると・・（参考書、p99）
> measurement
DATE DAILYSNOW SNOWDEPTH YEARLYSNOW FEB_SNOW
1 1 18 128 634 18
2 2 7 120 641 25
3 3 18 126 659 43
4 4 2 115 661 45
5 5 3 110 664 48
6 6 55 149 719 103
7 7 15 137 734 118
8 8 0 125 734 118
9 9 6 120 740 124
10 10 2 113 742 126
>
＊＊＊＊＊
データフレームmeasurement から、2月初旬の一日の降雪量と積雪量、それぞれの平均を「Ｒ」で計算すると
> mean(dailysnow)
12.6
> mean(snowdepth)
124.3
> median(dailysnow)
6.5
> median(snowdepth)
122.5
標準偏差はそれぞれ
> sd(dailysnow)
16.37206
> sd(snowdepth)
11.75727
>
補注
一日の降雪量の標準偏差が16cmと、平均値12cmに比べて非常に大きい。・・ということは、どか雪が降ることがあるということだ。一方で積雪量の標準偏差は12cm。120cmの積雪量に対して大したことはなさそう・・・10日間で均してしまうと、この時期の積雪量は一向に減らないというような解釈でよいのだろう・・・と思うのは、やっぱり大間違い。実際には2月6日の大雪のために積雪量は一晩で40cmも増えており、これだと、パイプハウスのフレームはおろか木造建築の屋根さえ壊れてしまう危険が迫っているのだ・・・。パイプハウスや木造家屋の屋根の雪に関しては、10日間で平均したり標準偏差をとったり、という作業は迂遠で、役立たずである。局所（その一日だけ、あるいはある数時間だけ）のデータがパイプハウスの倒壊・非倒壊を決める頼りのデータということになる。雪国の厳しさを思い知らされる・・・しかしまあ、・・・「Ｒ」の使い方とは離れてしまった。
＊＊＊＊＊
補注 さて、
#各一日で積雪量の変化に関して表を作り、それが降雪量と有意に相関するかどうか調べてみよう。
#各一日で積雪量の変化をDIFF_DEPTH=d_depthとおいてみる。
> d_depth < - diff(snowdepth)
＃diff()は前進差分のベクトルを生成する関数である。（参考書、p82の表を参照）
> d_depth
-8 6 -11 -5 39 -12 -12 -5 -7
と無事に数列はでてくる。が、差分なので個数は9個になる。
DATE=day, DAILYSNOW=dailysnow, SNOWDEPTH=snowdepth, それぞれの個数は10なので、工夫が必要である。とりあえずは、DIFF_DEPTH=d_depth の最後に10個目のデータとしてNA(not available、参考書、p88参照)を入れられたらよい。どうするか？ 参考書p92で行列の要素（成分）を取り出す方法が一覧に示されているが、要素を加える方法が書かれていない。後の課題ということにしたい。
今回は、姑息的ではあるが、とりあえず最初の日に（公開されている表を読み込んで）2を入れてみる。マニュアルでコンマを入れてベクトルデータとした。余り賢くない。たとえば、データの個数が膨大になれば不可能、つまり応用の利かないやり方である。
> d_depth < - c(2, -8, 6, -11, -5, 39, -12, -12, -5, -7)
> measurement < - data.frame(DATE=day, DAILYSNOW=dailysnow, SNOWDEPTH=snowdepth, DIFF_DEPTH=d_depth)
> measurement
DATE DAILYSNOW SNOWDEPTH DIFF_DEPTH
1 1 18 128 2
2 2 7 120 -8
3 3 18 126 6
4 4 2 115 -11
5 5 3 110 -5
6 6 55 149 39
7 7 15 137 -12
8 8 0 125 -12
9 9 6 120 -5
10 10 2 113 -7
>
とりあえずこれでデータ・フレームは作れた。次は、降雪量と積雪増加分との相関を出したいのだが、上記の表を一見しただけでも、大ざっぱなことを言えるかどうか、という程度の相関のように見える。ともかく、せめてグラフに散布図でプロットする。
> plot(dailysnow, d_depth) #参考書、p157参照。
相関係数を求めてみる。（#参考書、p186）
cor(dailysnow, d_depth, method=”pearson”)
> cor(dailysnow, d_depth, method=”pearson”)
0.945296
#余計かもしれないがスペアマン法でもデータをだしてみると・・
> cor(dailysnow, d_depth, method=”spearman”)
0.6932515
検定を行ってみる。（#参考書、p186）
> cor.test(dailysnow, d_depth, method=”pearson”)
Pearson’s product-moment correlation
data: dailysnow and d_depth
t = 8.1962, df = 8, p-value = 3.667e-05 （t値、自由度、p値）
alternative hypothesis: true correlation is not equal to 0 （対立仮説）
95 percent confidence interval:
0.7797842 0.9872987
sample estimates: （データから求めた相関係数）
cor
0.945296
結果は、p-value = 3.667e-05 となっているので、「相関係数が0である」という帰無仮説は棄却され、「相関係数が0でない」ことが分かった。
回帰直線 Ｙ＝Ａ・Ｘ＋Ｂ を求める
d_depthを縦軸Y、dailysnowを横軸Xと考えるのが普通だろう。
「Ｒ」では関数lm(Y~X, data=データ名)を用いることで、回帰直線を求めることができる。「Ｒ」で計算すると、
> result < - lm(d_depth ~ dailysnow)
> result
Call:
lm(formula = d_depth ~ dailysnow) （実行した式）
Coefficients:
(Intercept) dailysnow
-12.4553 0.8853
(intercept: 回帰直線Ｂの値 dailysnow：回帰直線Ａの値）
Ｙ（y=d_depth）＝ 0.885・X（x=dailysnow) – 12.5
この直線で、y=0 となるのは x=12.5/0.885=14.1 すなわち、一日に14cm以上の降雪があれば積雪量は増える。それ以下だと積雪量は減る、という計算結果。
考察： 今年の2月の雪は重たい雪であった。最高気温もほとんどすべて氷点下以下。雪は溶けたのではなく、コンパクトに圧縮されたのだ。一見厚みが増さなくても、どんどん重くコンパクトになって、屋根が潰される危険が高まる。
1メートル以上もの雪が屋根に積もっている場合、一日に雪が14cmも降っても、積雪量はそれほど変化せず、まだ大丈夫かと・・・思うのは、危険。積雪の厚みが増さなくても実はどんどん重くなっている・・・怖ろしい重さで屋根に乗っかっていると思わなくてはならない。実に、この2月、そして3月で、私の家の屋根や前室その他が大きく壊れてしまったのであった。
＊＊＊＊＊
関数plot(mydata)で散布図を描いた後、関数lm()の実行結果resultを用いて、abline(result)とすると、散布図に回帰直線を追記することができる。赤字で回帰直線を追記すると・・・（参考書、p188）
plot(dailysnow, d_depth)
abline(result, col=”red”)
＊＊＊＊＊
＊＊＊＊＊＊＊＊＊＊
