R スクリプト （日本語版）

英語版はここに。
Konishi, Tomokazu (2023):
Drastic changes before the 2011 Tohoku earthquake, revealed by exploratory data analysis.
figshare. Dataset.
https://doi.org/10.6084/m9.figshare.22010279.v1

#######
# R script
#######

# Rはこうしたスクリプトを書いて使うものである、エクセルのように触っていれば使えるようになるわかりやすさはないが、
# いちどスクリプトを書けばそれを使いまわしていくことができる。

# まことに残念ながら、数週間もするとスクリプトは「他人が書いたもの」と同程度の読みやすさになりがちである。
# だから（将来の使いまわしのことも考えて）読みやすく書くべきだ。
#
# スクリプトを編集するためには好みのテキストエディタ―を使えばよい。
# ワープロソフトはしばしば１バイト文字と２バイト文字を入れ替えてしまう、また大文字と小文字を勝手にいれかえがちである。
# Rはそれらを別のものと認識する。だからワープロで書くことはできない。
#
# Rは # のサインがあるとそこから右の文章を無視する（色指定で例外はあるものの）。
# だから # を注釈をいれるときに使える、ここでもそれを使う。
#
# まとまったタブと空白は無視される。そこでスクリプトを分かち書きするときにこれを使える。
# 注釈と分かち書きはスクリプトの読みやすさを決定するので重要である。
#
#######
# Download and Install R
#######

# RはCRANというサイトからダウンロードできる。
# https://cran.r-project.org/
#
# Linux, MacOS, そしてwindowsに対応している。
#
# ダウンロードしたプログラムを解凍してダブルクリックすればインストールされる。
# インストール先はサジェストされるが、ただし、もしコンピュータに「２バイト文字の名前」をつけてあるばあい、
# これがスムースに行われないことがある。そのときは、コンピュータのCドライブの直下に、たとえばRというフォルダをつくり、
# そこにインストールすると良い。
#
#
#######
# R 基礎知識
#######

# 使用するフォルダ（Linuxの文化圏なのでこれをdirectoryと呼ぶ）を決めておいて、そこを指定して計算する。
# 計算結果と過程はそのフォルダに半ば自動的に保存される。
# 特に保存したい結果は指定することで保存できる（あとで例を示す）。
# また読み込みたいデータは、テキストの形で読み込むことができる。
#
# もっともかんたんなディレクトリの指定方法は、空の .rdata というファイル（記録）をそのフォルダに用意しておいて、
# ダブルクリックしてRを立ち上げることだ。ここではデータとともにそれを用意して圧縮してある。
# Rをインストールしたら、解凍したフォルダのなかの.rdataをダブルクリックすればいい。
#
# Rで使われるデータのひとまとめのことをオブジェクトと呼ぶ。これは様々なものがありえる、数値、文字列、日付、ベクトル、行列、リストなど。
# 代入は次のいずれかで行う。Rコンソールに以下を書き込むか、テキストエディタ―からコピーアンドペーストすればよい。
#
aaa <- c(2, 3)
3 -> bbb
ccc = 3

オブジェクトに使う名前はアルファベットと数字であるが、数字を接頭語にはできない。
# aa1はOKだが1aaは不可である。大文字と小文字は別のものとして認識される。
# たまに pi のように使用されているものがある、またTRUEとFALSEの略称であるTとFは使うべきではない。
# いずれにしても、オブジェクトにわかりやすい名前をつけておくのはスクリプトの既読性に役立つ。
#
# R ができる機能はfunctionというオブジェクトで提供されている。
# f()
# という形で、fはfunctionのなまえ、カッコの中はそれに供する数値なりオブジェクトなりを指定する。
# たとえば上の c() というfunctionはCombine Values into a Vector or Listという機能で、ここではベクトルをつくっている。
#
# ここで
ls()
というfunctionを使うと、いま登録されているオブジェクトをリストにして返してくれる。
#
# またいま自分がどこにいるのか迷ったときは
getwd()
# で、現在のWorking Directoryの場所を教えてくれる。
#
# これらfunctionの情報はたとえば
?getwd
# と入力するとヘルプページが開かれる。またほとんどのケースで、google searchをすれば誰かがチュートリアルを書いてくれている。
# Rはあらゆるプログラミング言語のなかでもきわめて親切なコミュニティをもっている。
# そもそも、基本的に高価であるはずの統計パッケージが無料で配布されているところが、たいへんありがたいことだ。
#
#######
# データの読み込み
#######

# ここで用いるデータは2022年2月のものだ。
# 地震のID、Time、そしてmagnitudeがタブ区切りテキストで提供されている。
# これそのものはエクセルでつくり、テキストエディタ―から出力されている。
#
# これをRに読み込むときはread.table()関数を用いる。
#
data_example <- read.table(file="data2202.txt", sep="\t", header=TRUE)
#
# read.tableのカッコの中身は、ファイルの名前、区切り方、ヘッダーのあるなしの指定である。
# カッコのなかでいろいろなものを指定できる。
#
# ちなみに、地震のID、Time、area, 北緯、東経、深さ (km)、そしてmagnitudeがカンマ区切りテキストで提供した例をのせる。
# areaが2バイト文字で書かれている、そしてmagnitudeにひとつ欠損がある。これは2022年3月のデータだ。
# 2バイト文字などが含まれる場合すこし工夫をしないとエラーがでやすい。
data_example3 <- read.csv(file="data2203.txt", sep=",", header=TRUE, fileEncoding="sjis")
# オブジェクトも２バイト文字でつくれないことはないが、やめておくほうが無難だとおもう。

# ちなみに2022年3月は、2011東北震災の余震とおもわれるM7程度の地震があり、またその余震が多くあった。
# こうしたときにどうなるのかを見るためにも確認をおすすめする。
# ただし余震が減衰したこと、およびパラメーターがあまり変化しなかったことから、危険性は低かったと私は考えている。

# いまつくったdata_exampleはリストという、もっとも柔軟性のたかいオブジェクトになっている。
# is.list(data_example) # その確認
#
# Rとしては現時点でこれを行列とは認識していないが行列とおなじようにデータを指定できる。
# たとえばdata_example[3,2] で、data_exampleの2列、3行目のデータを呼び出せる。

#######
# 時間間隔
#######

# データの数はいくつだろう？
n_data <- length(data_example[,1])

# このn_dataをつかってデータの行を指定し、ひとつ隣との差をもとめる。
time_interval <- ( as.numeric(as.POSIXlt(data_example[2:n_data,2]))-as.numeric(as.POSIXlt(data_example[2:n_data-1,2])) )/60/60

# 3:5 はRのプログラミング上で便利な記法で、これで3,4,5という数列になる。
# 3:5-1で2,3,4になるので、ひとつ前のデータを指定できる。

# as.POSIXltは与えられた日付と時間をRに認識させる関数である。それをas.numericにかけると秒をあらわす整数値になる。
# その差はだから秒になる。これで時間間隔(h)を求めることができる。ちなみにas.POSIXltのままで計算させると、
# 計算結果の単位（attributionとして与えられる）がminになったりhになったりするので間違いのもとである。

####
# ヒストグラム
hist(time_interval)
# 左にかたよったグラフであることがわかる。

####
# QQplot
# 理論値は
ideal_exp <- qexp(ppoints(n_data-1), rate=1)

# これとソートしたデータを散布図にする
plot(ideal_exp, sort(time_interval))

# 一時近似した線をいれる
regression<-line(sort(time_interval)[50:140] ~ ideal_exp[50:140])
abline(coef(regression))
# ablineは直線を描くfunction, coefはlistであるregressionから数値を抜き出すfunctionである。

# 切片と傾きを表示する　　すこしだけ切片がずれていることがわかる
text(4, 5, regression[[2]][1])
text(4, 6, regression[[2]][2])

# 以上をpngフォーマットで打ち出す
png(width=1000, height=1000, pointsize=30, file="202202_interval_QQ.png" )
par(lwd=2, mex=0.7, cex=1, mai=c(1.8,1.8,1.2,0.5))
plot(ideal_exp, sort(time_interval), xlab="ideal", ylab="sorted data" , main="202203")
abline(coef(regression))
text(4, 5, regression[[2]][1])
text(4, 4, regression[[2]][2])
dev.off()

# png function にいれる注釈で画面の大きさ（ポイントサイズ）と字の大きさを指定している。
# parというfunctionで線の太さ、マージンの大きさを指定している。
# dev.off() は描画を終了するようにという指示　

####
# moving averages
# その時間から５日前の移動平均を求める

# 入れものを用意する
movingM<-movingnum<-movingInt<- 1:n_data+NA

# いつから計算できるか？
now<-1
fromwhich <- length( which(as.numeric(as.POSIXlt(data_example[,2] ))< as.numeric((as.POSIXlt(data_example[1,2])) +60*60*24*5)) ) #5日あと

for(i in fromwhich:n_data) {
now<-i
fromwhich <- length( which(as.numeric(as.POSIXlt(data_example[,2]) )< (as.numeric(as.POSIXlt(data_example[i,2])) -60*60*24*5)) )#5日まえ
avg<- mean (data_example[now:fromwhich,3], na.rm=T)
movingM[i] <- avg
avg<- mean (time_interval[now:fromwhich], na.rm=T)
movingInt[i]<-avg
movingnum[i] <--fromwhich+now+1 }

# ここで使っているのがfor ループ。iをfromwhichからn_dataまでひとつずつ増やしていって計算させている。{ }で囲まれているあいだが繰り返される。

# 結果をテキストで書きだす
write.table(cbind(movingInt, movingnum, movingM), "202202movingaverages.txt", sep="\t")

# 結果をRの記法で書きだす
save(movingInt, movingnum, movingM, file="202202moving.rdata")


####
# time course

plot(as.POSIXlt(data_example[2:n_data,2]), time_interval )

# 20 Febあたりに群発地震がある。これをさけてパラメータをみたい。そのためにこの位置をグラフをクリックして確認する。

locator(1) # 図の場所をクリックするとその位置を答える
# $x
# [1] 1645231637　　ｘ軸は秒でかかれているのでとても大きな数字になっている。

length(which(as.numeric(as.POSIXlt(data_example[,2])) < 1645231637))
# [1] 90
# 1番から90番までが影響を受けていなそうに見えた

# この間で５日間の平均回数はいくつくらいだろう
mean(movingnum[1:90], na.rm=T)
# [1] 27.8209

# このあいだの平均はいくつだろう　これの逆数がλだ
mean(time_interval[1:90], na.rm=T)
# [1] 4.862593 これはまたSDでもある

# ここで２σひくいのはどのくらい？
4.862593- 2* 4.862593/sqrt(27.8209)
# [1] 3.018799

# Rではこうした電卓的な計算はそのまま入力すると答えがかえってくる

# 以上をグラフにする　まずグラフの範囲をきめる
xlim=range(as.numeric(as.POSIXlt(data_example[,2])), na.rm=T)
ylim=range(time_interval, na.rm=T)

# 　na.rm=T は、NAがあっても無視するようにという指示

png(width=1000, height=1000, pointsize=30, file="202202_interval_time.png" )
par(lwd=2, mex=0.7, cex=1, mai=c(1.8,1.8,1.2,0.5))
plot(as.POSIXlt(data_example[2:n_data,2]), time_interval, xlim=xlim, ylim=ylim, xlab="date" , main="202203")
par(new=T) # 重ね書きを指示
plot(as.POSIXlt(data_example[,2]), movingInt, xlim=xlim, ylim=ylim, type="l", col="blue", xlab="", ylab="", axes=F)
abline(h=3.018799, col="green4")
dev.off()

#######
# マグニチュード
#######
##
# ヒストグラム
hist(data_example[,3])

####
# QQplot
# 理論値は
ideal_norm <- qnorm(ppoints(length(sort(data_example[,3]))), mean=0, sd=1)

# これとソートしたデータを散布図にする
plot(ideal_norm, sort(data_example[,3]))

# 一時近似した線をいれる 小数の、しかしSDが高いデータが混じっているのが読み取れる（たぶん沖縄）。


regression<-line(sort(data_example[,3]) ~ ideal_norm)
abline(coef(regression))
# 以上をグラフにする

png(width=1000, height=1000, pointsize=30, file="202202_M_QQ.png" )
par(lwd=2, mex=0.7, cex=1, mai=c(1.8,1.8,1.2,0.5))
plot(ideal_norm, sort(data_example[,3]), xlab="ideal", ylab="sorted data" , main="202203")
abline(coef(regression))
text(1,3, regression[[2]][1])
text(1,2.5, regression[[2]][2])
dev.off()


####
# time course

plot(as.POSIXlt(data_example[,2]), data_example[,3] )

# 群発地震のまえの平均とSDをロバストに求める
mean(data_example[1:90, 3], na.rm=T, trim=0.1)
# [1] 3.542466
# これはトリム平均という計算方法　データ数が少ないときにmedianより有利
mad(data_example[1:90, 3], na.rm=T)
# [1] 0.59304

# ここで３σ高いのはどのくらいだろう？
3.542466+ 3*0.59304/sqrt(27.8209)
# [1] 3.879769


# 以上をグラフにする　まずグラフの範囲をきめる
xlim=range(as.numeric(as.POSIXlt(data_example[,2])), na.rm=T)
ylim=range(data_example[,3], na.rm=T)

png(width=1000, height=1000, pointsize=30, file="202202_M_time.png" )
par(lwd=2, mex=0.7, cex=1, mai=c(1.8,1.8,1.2,0.5))
plot(as.POSIXlt(data_example[,2]), data_example[,3], xlim=xlim, ylim=ylim, xlab="date" , main="202203", ylab="magnitude")
par(new=T) # 重ね書きを指示
plot(as.POSIXlt(data_example[,2]), movingM, xlim=xlim, ylim=ylim, type="l", col="blue", xlab="", ylab="", axes=F)
abline(h=3.879769, col="#33aa33aa")
dev.off()

# この１か月のあいだに２度も３σに迫っていた。これらはどうも沖縄で活発になっているのを反映したらしい。
# 最後のablineのカラー指定に # を使う、半透明な色を使用してみている。このときは#は無視されない。
# この指定方法をつかうと色盲に対応するような色を選べる、たとえば"red"のかわりに金赤 "#e8380d"など。

# ここではR本体だけで計算できるものを紹介した。しかしRには多くのパッケージプログラムが用意されている。
# それらは、なにか特殊な用途で便利なように作られて維持されている。たとえば緯度と経度を扱いたいときは
# parzer と geosphere
# というパッケージが便利である。グーグルサーチすればこれらの入手方法や使い方のインストラクションが手に入るはずだ。

##############
# making Table 1
##############

# P-values in R’s notation is 1-pnorm(q= M, sd= scale, mean= location), for reference.

1-pnorm(q= 5, sd= 1.4, mean= 3.8)

# [1] 0.195683

# to make a table,
write.table(cbind(6:18/2, 1-pnorm(6:18/2, mean=3.3008, sd=0.5883), 1-pnorm(6:18/2, mean=3.8, sd=1.4)), file="normal.txt", sep="\t")