今回のデータセットには日本人のお名前が冠せられてます。調べてみたら、「廣瀬先生」を知らなかったら「信頼性、統計業界」じゃモグリと言われそうな大先生でした。知らなかった私は確実にモグリです。廣瀬先生のせいじゃありませんが、データセットの解説ページの一文字が誤っているみたいです。そのフェイントにモグリの老人は右往左往。
※「データのお砂場」投稿順Indexはこちら
hiroseサンプルデータセット
Rのパッケージ「Boot」に含まれるサンプルデータセットをabc順に経めぐってます。今回はhiroseとな。以下がサンプルデータセットの解説ページへのリンクです。
ガス絶縁の変圧器の絶縁材料として使われるPETフィルム（ポリエチレンテレフタレートということでよいのだよね、書いてないけど）の「寿命」加速試験の結果のデータセットみたいです。voltの欄に5とか15とかさらっと書いてますが、単位kVです。10V電圧にたじろぐ弱電な老人は恐れ入るばかりであります。元データは、H. Hirose先生の1993年論文で、IEEE Transactions on Reliabilityに掲載されたものみたいです。ただ、上記のページに説明されている以下引用の「１」という値に騙されましたぞ。詳しいことは後で。
The censoring indicator;
1means right-censored data.
廣瀬英雄 先生
調べてみたら、先生のホームページが直ぐにヒットしてきました。
https://hirosehideo.com/welcome.html
九工大に長らくおられた（現在は既に退官）先生で九工大の名誉教授です。なお上記のホームページからサンプルデータセットの元になった論文の日本語版を拝見することが可能です。URLは以下です。
https://hirosehideo.com/TakaokaReview/1995_143_1.pdf
なお、上記の日本語論文は「高岳レビュー」なるコーナーに配置されてますが、これは先生が、株式会社高岳製作所の出身であるためみたいです。高岳製作所といえば、断路器、変圧器、など電力伝送系分野のオオモノを製造している会社であったと（2014年に、東光電気株式会社と完全統合され、東光高岳となっているみたい。）
元論文、さらに日本語版が拝見できるけれど、まずはいつもの通り、サンプルデータセットに飛び込んであがいてみることにいたしましたぞ。
まずは生データ
いつものように生データをロードして概要をみてみます。
データはシンプルなデータフレームで、電圧を記したvolt、時間（hour）を記したtime、そして右打ち切りデータであるか否かを記した censの３カラムからなります。
この時点でじゃ、解説ページの１が右打ち切りのお印であるとの情報にとらわれておる老人です。
とりあえずプロットしてみました（この時点では廣瀬先生の意図などまったく理解していなかったので、X、Yヒックリ返ってます。）良くないプロット操作が以下に。
plot(time ~ volt, data = hirose, main="hirose, Failure Time of PET Film", ylab="time[hour]", xlab="volt[kV]")
プロット結果が以下に。
あちゃ～、5kVと7kV以上のところでデータが懸絶してます。それに上記では右打ち切りとそうでないデータが混ざっているし。。。
そこでaggregateしてみることに。操作はこんな感じ。
aggregate(hirose$time, by = list(volt=hirose$volt, cens=hirose$cens), mean)
結果が以下に
むむっ、右打ち切を示すcens=1ばかりじゃない、0なのは5kVのときだけ？（これが大間違いです。）
上記のvoltとcensの組み合わせで「層別」して、各分類ごとにボックスプロットしてみることに。dplyr使用。
v5_0 <- hirose %>% filter(volt==5, cens==0) %>% select(time) v5_1 <- hirose %>% filter(volt==5, cens==1) %>% select(time) v7_1 <- hirose %>% filter(volt==7, cens==1) %>% select(time) v10_1 <- hirose %>% filter(volt==10, cens==1) %>% select(time) v15_1 <- hirose %>% filter(volt==15, cens==1) %>% select(time)
まず値が他と懸絶している5kVのときのboxplot。操作は以下に
boxplot(c(v5_0, v5_1), names=c("5kV_0", "5kV_1"), ylab="time[hour]", main="hirose, Failure Time of PET Film")
プロット結果が以下に。
なんじゃ、右打ち切りのハズの1の方がばらついていて、cens=0の方が、スッパリ一つの時間でキッチリ打ち切られている感じ！
念のため7kV以上についてもboxplot。
boxplot(c(v7_1, v10_1, v15_1), names=c("7kV_1", "10kV_1", "15kV_1"), ylab="time[hour]", main="hirose, Failure Time of PET Film")
結果が以下に。
ここに至りて、モグリの老人もようやく気が付きました。
右打ち切りデータはCEN=0のときじゃん
bootパッケージの中には、「右打ち切りデータを取り扱うためのブートストラップ法関数」censbootというものがあるので、多分それを適用せよとの、サンプルデータセット編集の中の人の思し召しなのでしょうが、CEN=1はフェイント。
廣瀬先生の論文自体は、勝手な理解でまとめてしまうと「いったい何Vだったら壊れれずに使えるのか」を計算するためのもの。グラフで描けば
plot(volt ~ time, log="xy", data = hirose, main="hirose, Failure Time of PET Film", ylab="volt[kV]", xlab="time[hour]", pch=16, col="red")
以下のように、縦軸の電圧を下げていくと、どんどん寿命が延びていく（だから5kVのときに右打ち切りデータが現れる）けれども、多分、5kVよりも下に、ほぼほぼ寿命が実質無限に見えるような電圧がありそう、ということかな。勝手なこと書いているけれど。
まあ、詳しいことは先生の元の論文をご参照あれ。
またまたブートストラップ法のブの字も使っておらんな。モグリだし、しかたないか。