R言語 (2)

過去にこの血圧データを分析するのを色々な言語でやったぞ。scheme,python,haskell,go-lang,awk等だ。そして今回はR言語って訳。自分の分身データだと気合が入るってもんです。

前回のスクリプトを少し整理した。

bld <- read.csv("bld.csv", header = F)
colnames(bld) <- c("yy","mm","dd","hh","hi","lo","pl")
am  <- subset(bld, bld$hh < 12)

m.am <- aggregate(cbind(hi,lo,pl) ~ mm + yy, data=am, mean)

まず、元データのbld.csv内の年の扱いが西暦下２桁だったのを４桁に改めた。sedを少し改変して年の頭に"20"を追加するだけ。ヘッダーは無し。だから、Rスクリプト中で読み込む時、header = F を指定。 FALSEが正式名称だけど、省略してます。

ヘッダーに変わる列名を colnamedsで追加してます。unix文化では、データに余計な物を付けないんで、受け取った方で何とかするのが正しい方法だと思うよ。

月毎に平均を取る方法、前回から少し変えて、最高血圧、最低血圧、脈拍も同時に計算させてます。更に、グルーピングの順番を月、年の順にした。こうしておけば、後で並べ替えは必要なくなります。

> head(m.am, n=3)
  mm   yy       hi       lo       pl
1  1 2012 125.0645 79.32258 59.29032
2  2 2012 123.2759 79.48276 59.13793
3  3 2012 120.9667 78.43333 61.83333

副作用が一つあって、列順が、月年になってしまっています。まあ欧米風と笑って許しましょう。それから、R語での変数名は、キャメルケースでも無く、C語みたいにアンダーバーで結合するでもなく、ドットで結合するって方式です。

なんか今までの風潮だと、ドットはメソッドチェーンだとか、pythonみたいにパッケージ名と関数名の分離に使うと思ってた。慣れるとR語風が快適と思えるようになってきましたよ。

ここまで、たった４行書くだけで、上のようなデータが手に入り、R語の破壊力に感心しましたよ。これってベクトルを中心にするR語の威力なんですなぁ。

それから、emacsでESSを使ってる人用のTipです。実行したい行にカーソルを持って行って　C-c C-cします。すると、空行に出会うまで、連続して実行してくれるんです。これで、途中々、実行結果を確認しながら、コードが書けます。何気ない気遣いに感謝んです。

ここまでが、前回からの復習。いよいよ時系列って事で、 Rと時系列（１） を、教科書にします。オブザーバーとして、実データを時系列解析して結果を考察

まずは、時系列形式を作るんだな。

hi.ts <- ts(m.am$hi, start=c(2012,1), frequency=12)
pl.ts <- ts(m.am$pl, start=c(2012,1), frequency=12)

データ列は、年月順に並べてあるもの。そしてデータの開始年、月を指定。月は12で繰り上がるよってヒントを与えてあげるんだな。

> str(hi.ts)
 Time-Series [1:104] from 2012 to 2021: 125 123 121 118 119 ...
> end(hi.ts)
[1] 2020    8

そのように解釈してくれたか確認。strって関数で要約を確認。普通strったら文字列関連かと思うけど構造ってのの略なんだな。で、年は2021年まで続いているって。構造の詳細を確認。そう、そこまでのデータが有るんよ。

> hi.ts
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2012 125.0645 123.2759 120.9667 118.1667 118.9032 115.8276 113.2903 117.5484 ..
2013 124.0323 119.7857 116.8065 117.0333 118.6154 120.6333 118.9355 119.4194
 :

一応、検算してみるかな。

 > colMeans( subset(m.am, yy == 2012 & mm == 8) )
         mm         yy         hi         lo         pl
    8.00000 2012.00000  117.54839   76.29032   58.90323
!> colMeans( subset(m.am, yy == 2013 & mm == 2) )
        mm        yy        hi        lo        pl
    2.0000 2013.0000  119.7857   75.0000   57.0000

どうやら問題なさそう。でも問題は、 こんな表を見せられたって、眼がチカチカするだけ。図にしましょ。

ts.plot(hi.ts, pl.ts)
plot(decompose(hi.ts))

最初の図は、取り合えずの出来具合を、同一の画内に表示。ふと思ったんだけど、最高血圧と脈拍って、単位が全く違う。そんな数値を同一に表示していいのか。

２番目の指示で、前回の教科書と同じ形式の図が出てきた。元データ、トレンド、周期性データ、誤差分。こうしてトレンドが年齢と共に上昇してるのが残酷なまでに表示されちゃったぞ。周期性では、１月が一番高く春が一番低いって結果。寒い時は気を付けましょうって事だな。

acf(hi.ts)
acf(hi.ts, lag.max=96,ylim=c(-.2,1.0))

自己相関。下のやつは、lagを広範囲に振った時。いずれも前月との強い相関が見られますた。

ccf(hi.ts, pl.ts)

自己相関だけじゃつまんないので、最高血圧と脈拍の相関も取ってみた。X軸はlagになってて、０のポイントで、強い負の相関って結果になった。まあ、上で並べて表示した時、予想はされていたけどね。

血圧が低いと脈拍を多くして、体内に送り出す血液量を増す必要がある。血管内にゴミが溜まって細くなると、圧力を高める必要が有る。血圧が上がると脈拍を低下させてエンジンである心臓をいたわろうと言う挙動になるのかな？　素人の推測だけど。

流量なんてのが出てきたんで、R語に備わっている流量データに注目してみる。有名なのはUKgasだ。これ、あちこちで使われているんで、少々秋田。で、他のを探したら、ナイル川の年間流量データってのが有った。

そんなに大きな川の流量なんて、どうやって測るねん？　年間って事だと、例えば１時間あたりのの流量を求めて、それを加算してくのかな。気が遠くなる。たとえ１時間の流量を求めるとしたら、どうやる？　根本原理を考察（素直に信じるのは良くないと思う）。

素人考えだけど、川の断面積を求める。川の流れるスピードを求める。それを掛け算すれば流量になるな。でも、流れは水面近くと底の方で違うだろうに。なんかそういうのの公式みたいなの有るのかな？　軽く検索したら、河川の土木事務所とかで使われているEXCELのアプリが有るみたい。公式なり、表を引いて計算してるんだろうね。

で、このナイル川の流量データには少々問題が有るようである。それを統計でカバーしましょって例みたい。

> example(Nile)

Nile> require(stats); require(graphics)

Nile> par(mfrow = c(2, 2))

Nile> plot(Nile)
Hit <Return> to see next plot:

Nile> acf(Nile)
 :

こんな風に、取り合えず実行してみてから、コードを眺めてみるのがいいのかな。これぞ現場実習です。で、現場ではついて行けず。なんせナイル川は遠すぎます。日本の河はどうよ。

ナイル川には負けるけど、日本で一番長い川、信濃川の流量データでもないかしら？

データで見る信濃川 って事で、信濃川スゲーってのが強調されてるけど、年間推移は見つからず。 ああ、ここに有った。日流量年表検索 。

なお世界一の河はアマゾン。日本もその流域の一部になってて、既に水没してますよ。中小企業は、軒並み被害にあってる。政府よ、アマゾン被害から企業を守る堤防を早く築けよ。

で、信濃川のデータは細かすぎる。年を指定してDLして、、、なんてやってられないね。まとめたやつは、日本河川協会からDVDになって発売されてる。昔から水の利権は争いの種だったけど、データにも利権が発生してるのね。

しょうがない河岸を変えよう。

同じ国土交通省の管轄だけど、気象庁はどうよ？　気象庁で分かり易いのは、気温だな。

よく天気予報のおねーさんが言ってる、先月は平年に比べて２度高かったとかの平年。 過去30年のデータを集めてきて、並べ替え。それを３等分して真ん中の組のデータだけを使って計算するとか。極端に寒い（暑い）データを除外した上で求めているのね。そして、このデータは、５年だか１０年おきき更新するらしい。

こういう知識は気象予報士の試験に出るのかな。たまたま読んでいた小説に、リタイアした人が趣味で予報士試験に挑戦する話が出て来てた。難しくて３度の挑戦でやっと合格したとか。合格率5%とからしい。1000から1500時間の勉強が必要とか。暇を持て余した人はボケ防止に挑戦してみるのも面白い鴨。

ここにDL用のページが用意されてた。 過去の気象データ・ダウンロード

東京と富士山と言う、日本を代表する地域を指定。Fuji-yamaは、静岡県所属？　それとも山梨県？　よく出るクイズだな。観測所が有る所は山梨県でした。他に南極なんてのも有ったけど、無視。データの種類は、取り合えず温度データだけ。期間は1990年から先月までを指定。

結果を表示させて、問題なければCSVでダウンロード。優しいインターフェースで、小学生にも使えるぞ。夏休みの宿題のネタにどうぞって訳だな。

ちゃんとDL出来たか lireofficeで確認。文字コードをSHIFT-JISに指定しないと、文字化けするぞ。どこまで行ってもWindows仕様な訳ね。年・月を分離する記号は "/" になってるし、月も１桁だったり２桁だったり、EXCEL以外では不便な風。

ダウンロードした時刻：2020/09/10 05:30:27

,東京,東京,東京,富士山,富士山,富士山
年月,平均気温(℃),平均気温(℃),平均気温(℃),平均気温(℃),平均気温(℃),平均気温(℃)
,,品質情報,均質番号,,品質情報,均質番号
1990/1,5.0,8,1,-18.7,8,1
1990/2,7.8,8,1,-12.1,8,1
 :
2020/7,24.3,8,2,4.8,8,1
2020/8,29.1,8,2,7.6,8,1

生データを確認。品質情報とか付いてるけど、そんなのをオミットして取り込みたいぞ。 元データを改変する事なくね。R言語にはpipeが有る事を知った。

con <- pipe("sed '1,5d' data.csv | sed 's#/#,#' | cut -d',' -f1-3,6")
tmp <- read.csv(con, header=F)
colnames(tmp) <- c("yy", "mm", "tokyo", "fuji")

R言語でパイプから流れて来るデータを読めるのね。これは素敵な事。して、パイプでunix側に依頼したのは、まずsedで、余計な冒頭の５行を削除、データだけをむき出しにする。

次の段のsedで年・月を分ける"/"をカンマに変更。そうしておいて、cutコマンドで区切り記号をカンマに設定して、必要な項だけを抽出。

なお、２段目のsedで正規表現の区切りに"#"を使ったのは、見易くする為です。

!> summary(tmp)
        yy             mm             tokyo            fuji
  Min.   :1990   Min.   : 1.000   Min.   : 4.70   Min.   :-22.50
  1st Qu.:1997   1st Qu.: 3.000   1st Qu.: 9.20   1st Qu.:-14.55
  Median :2005   Median : 6.000   Median :16.90   Median : -5.70
  Mean   :2005   Mean   : 6.457   Mean   :16.53   Mean   : -5.94
  3rd Qu.:2012   3rd Qu.: 9.000   3rd Qu.:23.10   3rd Qu.:  2.45
  Max.   :2020   Max.   :12.000   Max.   :29.60   Max.   :  8.30
                                                  NA's   :1

東京暑いと言っても平均気温だと30度を割ってるのか。そんな事より気になるのは富士山のデータにNA（欠損値）が有る事。一体何時なんだ。

!> subset(tmp, is.na(tmp$fuji))
       yy mm tokyo fuji
 353 2019  5    20   NA

去年の５月か。全面的な改修工事でもあって、閉鎖されてたのかな？

sakae@pen:~/ANA$ grep '2019/5' data.csv
2019/5,20.0,8,2,,1,1

元データに当たると、確かにデータが欠損してた。品質情報（右から２番目）も、確かに、統計値がない（欠即）って説明されてたよ。

 > max(tmp$fuji)
 [1] NA
!> max(tmp$fuji, na.rm=T)
 [1] 8.3

欠損値が影響してくるような関数では、安全なために結果がNAになる。そのような関数には、NA除外オプションが用意されていたりする。NA値をどうするかってのは、悩ましい問題。色々な処理方法が提案されてる。例えば、Rによる欠損値処理 。やり出したらキリがないので先に進む。

tokyoのトレンドを表示させてみた。内心、じわじわと温度上昇してると思ったけど、結構長期 周期でふらふらしてた。

このトレンドって、観測データを平均化したもの（ですよね）。もう一回平均化したらどうなる？

 > x <- decompose(t.ts)$trend
!> plot(decompose(x))

t.tsはtokyoの時系列形式。これを一回処理してトレンドデータだけをｘに取り出す。それをもう一回decompose。結果は、更に滑らかになった。って、当たり前の事だな。 　 今回はこれぐらいにするか。やや、新しいページを見つけたぞ。

Rではじめての時系列分析（AirPassengers）

このページを辿って行くと、pythonでは、どうやるなんてのに案内された。pythonめんどくせーな。