nonukes2011

LSS(Life Span Study) 13報　被爆者データのRによる分析　

　市民にも科学を､という趣旨で分析のためのプログラムを公開します｡
といっても､今のところ玄人には簡単なことまで説明し､素人には難しくみえるかもしれないという微妙な位置づけになっています｡
　とりあえずプログラムについては､GPLライセンスで公開するので､適宜修正､公開してください｡

被曝された方から得られた大変貴重なデータであることは心して分析したいと思います｡

0.お急ぎの方

　色々説明していたら長くなったので､Rで実行したい人は下記から｡

LSS13 データのダウンロード

放影研のここから住所などを入力すればすぐにダウンロード可能｡
Rで読み込むのは下記のファイル｡エディタなどで開いてみるとよい｡

R13mort.dat

R13 data documentation.pdfが変数名ファイル｡みておく方がよい｡このファイルでの変数名とデータでの変数名大文字小文字が異なるので注意｡

Rの起動と作業ディレクトリの指定

上記の　R13mort.dat　があるディレクトリを作業ディレクトリに指定｡

私の提供するファイル

処理用プログラム

下記リンクからGoogle docで公開しています｡適宜コピー､Rにペーストして利用して下さい｡
2013/1/19追記　変数の定義ミスがあったので､読み込みスクリプト修正しました｡処理スクリプトは変更ありませんが､変数の定義が異なるので､当該変数の部分の結果も変化します｡お手数ですが､再度処理をお願いします｡こちらに修正前後の結果を比較しておきました｡線量部分は大きな変化はありません｡

@habari2011duniaさんから下記､ご指摘頂きました｡ありがとうございます｡

lssinc13forpubrd.R

R13mort.datを読み込み､必要な変数を生成してR用のデータセットを作成､保存するRスクリプト｡

lssinc13forpubprc.R

上でつくったデータセットを使って推定するためのRスクリプト｡
説明コメント､出力イメージつきなので､推定しない人もご覧になるとよいのでは｡

1.このデータについて

0)論文　
Preston, DL, Y Shimizu, DA Pierce, A Suyama, and K Mabuchi (2003), "Studies of mortality of atomic bomb survivors. Report 13. Solid cancer and noncancer disease mortality: 1950-1997," Radiation Research, 160 (4), 381-407. 例えばここからpdf　

　参考)　原爆被爆者における固形がんリスク - 放射線影響研究所

　で用いられたデータを分析するためのRスクリプト｡

ここ数報の研究について知りたい方は､ここにまとめたので参照｢放射線影響研究所のデータについて｣｡

　ここでは13報(以下LSS13論文などと称する)を扱いますが､14報も今年公開され(概要pdf)､データも既に提供中です｡線量評価13報はDS86､14報はDS02と異なっている､変数名も異なっている､疾病の分類も若干異なるといった点に注意が必要です｡

　しかし､結論=線形閾値無しモデルが最良=はかわりません｡そのうち､14報の分析スクリプトも公開するかもしれませんが､ここで公開しているものがわかれば自分でもできるでしょう｡

1)このデータの問題/限界
　以下を指摘できます｡
　a)被曝後5年から開始されたのでその間に亡くなられた方が含まれていない｡
　b)比較対象グループつまり､被曝していない両地域の人々が含まれていない(他のデータセットでは投下時には市外､その後入市された方が含まれているデータ､分析もある)｡

　c)(多分)外部被曝しか考慮していない｡
　d)年齢分布などに偏りがある｡例)被曝線量が高い人では､若い人の割合が高くなっている｡

　e)線量バッチなどつけているはずがないので､被曝量は爆心地からの距離やどこにいたかという情報に基づく推定値であり､それ自体誤差を含んでいる｡

　f)個人レベルのデータではなく年齢層､被曝量などの層別に集計されたデータである｡

　a)b)については(外部)被曝量の影響を過小評価する方向に､c)については(外部)被曝量の影響を過大評価する方向に作用する可能性が高い｡

　d)について傾向スコアで補正した推定を行ったところ､(敏感な若年層のウエイトを高線量部分について小さくするので)若干だが影響は小さくなった｡

Hamaoka, Yutaka (2011)"A Search for Better Dose-Response Function," 14th International Congress of Radiation Research, Warsaw, Poland, 27 Aug -1 Sep.

　　その後､sample selection bias を検討した論文があることに気づいたが未読｡
　e)については検証した論文があったが未読｡
　f)については､上記のPolandの学会で同研究所の児玉氏にお目にかかったときに言ったところコホートとはいえ､細かく区切られているので個人レベルでも結果はかわらないのではないかとのこと｡しかし､個人レベルのハザードモデルなど適用可能なモデルが広がるので何とかしたい｡　　
　｢放射線影響研究所のデータについて｣の下の方には､個人レベルでの推定例が紹介されている｡
追記)2013/1/18　ポアソン回帰など一般線形モデルの推定値の標準誤差は説明変数の共分散の逆数に比例｡連続変量を適当にカテゴリにまとめると当然分散は小さくなり､推定値の標準誤差は大きくなり､有意になりにくくなる｡

　これらの限界があるが､我々が簡単に入手可能(だが曲解?されている)ほぼ唯一のデータなので､分析してみる｡

　被爆者データでは線形閾値なしモデルが支持されているが､その他の研究を踏まえて国際的には100mSV以下では不確実とされているのかもしれません｡UNSCEARの報告書など読み込んでいないので､どのような根拠で､そのようになっているのか､もしくはなっていないのか現在調査中です｡

　なお､　原発従業者疫学調査について　にまとめたように年間被曝線量1mSv程度の原発従業員を対象とした調査でも線量と有意な関係が見いだされています｡

　予防原則という観点からも日本人を対象としたデータである､この結果を重視すべきであると考えます｡

2.推定されているモデルについて

1)ポアソン回帰モデル

　死亡者数がポアソン分布に従い､その平均値が各種の変数(年齢､性別､被曝量など)に依存するというポアソン回帰モデルで推定している(層別に集計したデータに対して)｡
　推定は有料のEpicure(簡単な解説)というソフトが用いられている｡しかし､オープンソースの無料･統計ソフトRを使えばポアソン回帰(やその他のモデルでも)推定可能(Mac, Win, Linux版あり)｡

　(あるコホートにおける)死亡者数yが下記のようなポアソン分布に従うと考える｡ここでλは(人年あたりの)死亡率｡

　P(Y=y)=exp(-λ)λ^y/y!

　この死亡率λが性別､都市､年齢(z)､被曝線量xなどによって影響されるとモデル化｡これを被曝線量に依存しない部分=ベースラインλ0(z)と､線量とzに依存するf(x,z)の積として表現｡λ0(z)に各種変数を導入する際は､負にならないように指数にのせる｡以下では過剰相対リスク(ERR)モデルを推定する｡

　　λ=λ0(z)*(1+f(x,z))=exp(b0+b1z1+...)*(1+f(x,z))

注　2015/7/8)λ0(z)*f(x,z)=exp(b0+b1z1+...)*f(x,z)を上記のように表記訂正(MAKIRINTAROさんのご指摘による)｡

　例えばb1>0であれば､z1が大きくなるほどλ=死亡率が高くなることになる｡

検定については下の方を参照｡

　上記論文では　f(x,z)の線量部分を以下のように定式化して推定｡線形(1次)モデルのあてはまりが最良であったと文章で記述している｡

線形(1次)

2次

βx^2

1次+2次

β1x+β2x^2
1次項､2次項の相関が高いので不適切

(手動)閾値モデル Tより低い線量では影響0というモデル｡

β(x-T) x>=T
0 otherwise

閾値を変更して推定｡それを繰り返すという原始的な方法｡閾値自体も推定量なので､そのものを推定して､統計的検定する方がよい｡

線量カテゴリ･ダミー(以下ノンパラ回帰と呼ぶ)

Σβi*xi 　　　xi　はi番目の線量カテゴリダミー｡上記のような関数形を想定せず､各線量毎に影響の大きさを推定｡カテゴリの区切り方に結果が依存しやすい｡

　･最終的には線形モデルが最良としている｡これが下記のように他の変数にも依存するというモデルを用いている｡　　x:線量

　　 f(x,z)=β*x*性別ダミー*exp{α1被曝時年齢/30+α2*log(到達年齢/70)}

#私が推定した結果　(線量部分については放影研の提供している結果出力ログとほぼ一致｡ベースラインについては多分､性別などのコーディングの違いによって推定結果が異なる｡→これはAMFITのマニュアルがないので不明｡

#[[1]]　　　この結果は上記のように変数を修正した推定結果(2013/1/19)
# beta.estimated ses t s
#(Intercept) -7.20114548 0.66308663 -10.860037 ***
#mhiro 3.78935443 0.82399285 4.598771 ***
#mnaga 3.76984999 0.82418426 4.574038 ***
#fhiro -0.03939173 0.03170574 -1.242416
#me30 -0.05684162 0.04259556 -1.334449 .
#fe30 0.21281297 0.03998091 5.322864 ***
#me30sq -0.03531584 0.02359363 -1.496838 .
#fe30sq 0.02496813 0.02235298 1.116993
#me30qsp 0.05365397 0.04115247 1.303785 .
#fe30qsp -0.08020425 0.03852640 -2.081800 *
#me50qsp -0.18927995 0.05968431 -3.171352 ***
#fe50qsp -0.06817465 0.05521142 -1.234792
#mlage70 8.94142250 1.63463851 5.469969 ***
#flage70 -2.63884525 2.26133457 -1.166942
#mlage70sq 2.12190266 1.22083589 1.738074 *
#flage70sq -4.91288034 1.81215840 -2.711066 **
#mlage40qsp -3.33159492 1.56302819 -2.131500 *
#flage40qsp 5.99872530 2.11289049 2.839108 **
#mlage70qsp -7.25067676 2.01427381 -3.599648 ***
#flage70qsp -2.05687939 1.66741196 -1.233576

#下記が線量に関連する部分

#exp(e30) -0.37602715 0.09371096 -4.012627 ***
#　　　　　　　　　　↑上式のα1　年をとって被曝するとリスクは低下
#exp(lage70) -0.69497460 0.45241326 -1.536150 .
#　　　　　　　　　　同α2
#x:dsex 0.25724104 0.09736845 2.641934 **
#　　　　　　　　　　↑性別による影響の差がある
#x 0.46774879 0.06249616 7.484440 ***
#　　　　　　　　　　↑これが線量の係数　放影研のLSS13の結果0.475とほぼ一致
#[[2]]
# LL AIC BIC
#1 -12879.08 25806.16 26263.13　→下にあるモデル1よりもAICは小さく､あてはまりは良好｡

2)この論文の問題点
　再推定したところ､線形が最良という結論は変わりませんが､以下の大小の問題があり､100mSv以下では影響がないという誤解もしくは解釈を生んでいると考えます｡

検定統計量などが明示されていない｡

線形､2次､線形+2次､手動閾値モデル､ノンパラモデル　を比較して線形が最良であった､と結論づけている｡しかし､各モデルの適合度(deviance､尤度､AIC等)が明示されていない｡(14報では線形､二次､線形+二次の比較についてはdevianceの差による検定が明示されている｡が､手動閾値モデルなどとの比較はされていない｡)
(線量以外の)個別のパラメータのt値なども明示されていない｡疫学ではt値ではなく信頼区間で表示するようですが､ベースライン部分の係数の推定値､信頼区間はまったく表示されていない｡

線量の1次項､2次項を同時に入れたモデルを推定している｡

これらの相関は0.9以上と高いので､同時に入れると多重共線性の問題が生じる｡曲線性の検定をしたい気持ちはわかるが､あてにならない可能性が高い｡

サンプルを限定した分析の限界を明示していない｡

全サンプルを用いた場合､上記のような結論｡それで終了でいいはずだが､低線量サンプルに限定した分析をしている｡サンプル数が減少すると､標準誤差も大きくなり有意でなくなりやすいが､その点についての記述が弱い｡これが100mSv以下無影響論の根拠になっている可能性がある｡

例えば下記の論文も全領域では有意｡線量域を限定した分析すると有意ではなくなることを示した｡ただし､単に統計的検定力が低下したためであると正しく解説している｡
GILBERT, E. S., KOSHURNIKOVA, N. A., SOKOLNIKOV, M. E., SHILNIKOVA, N. S., PRESTON, D. L., RON, E., OKATENKO, P. V., KHOKHRYAKOV, V. F., VASILENKO, E. K., MILLER, S., ECKERMAN, K. & ROMANOV, S. A. 2004. Lung Cancer in Mayak Workers. Radiation Research, 162, 505-516.
LSS13論文では線量を限定した推定の際も､他のパラメータは全データを用いて推定したものを用いている(そうしないと不安定になって推定できない)｡そうするならば以下に私がしているように､全サンプルを用いてある線量以下では効果が異なるというモデルを推定すべきである｡

内挿値における誤差についての記述漏れ

推定したパラメータを用いて､各線量毎の期待過剰死亡数を算出している(TABLE 10)｡これ自体も推定値なので推定の誤差があるが､それを明示していない｡
5mSV以下では過剰死亡数が0となっている｡これが､低線量では被曝による影響がないという解釈につながっている可能性が高い｡

例えばアリソン｢放射能と理性｣　のTable 4,5は　LSS13ではなく､Preston et al. (2004) Table7などを引用｡ http://www.bioone.org/doi/abs/10.1667/RR3232

以下の推定では､これについては対応していません｡

3.　このスクリプトでの推定

1)モデルの適合度

モデル全体のあてはまりの優劣

モデルの全体的あてはまりを比較する指標として､異なるタイプのモデルでも比較できるAIC(赤池の情報量基準)を用いる｡
ここでのルーチンでは　AICが出力されている｡これの値が小さい方が､データにあてはまりがよいモデルということになる｡

AIC　赤池の情報量基準　モデルのデータからのずれの大きさを表す指標なので､値が小さいモデルほど､あてはまりがよいモデルとなる｡用いたパラメータの数も考慮して補正されている｡

部分的な診断:変数が有意か否か

変数の横に　.とか*がついているか?

Rでは　.10%　*5%　**　1%　***0.1%なので注意｡　
これらの記号がついている変数は死亡率に影響を与えると考えることができる｡

厳密には係数=0:その変数と従属変数､つまり死亡率には相関がない)という帰無仮説が棄却されたことを意味する｡相関がないという帰無仮説が棄却されたということであって､相関があることを積極的に言明するものではない(背理法なので)｡厳密には因果であるともいえないことに注意｡

下記のモデル1の出力例
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) -9.1352303 0.0667514 -136.85 <2e-16 ***
　　右に. とか*がついている変数が影響を与えると考えてよい｡
#city 0.0128649 0.0229870 0.56 0.576
#sex -0.6048946 0.0208939 -28.95 <2e-16 ***
#agex 0.0012827 0.0008275 1.55 0.121
#age 0.0699884 0.0009362 74.76 <2e-16 ***
　　ここまでが線量に関係ないベースライン
#log(1 + col_d10) 0.6253448 0.0485312 12.88 <2e-16 ***
　　これが線量の係数　　*がついているので0ではないといえる｡

#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#(Dispersion parameter for poisson family taken to be 1)
# Null deviance: 27566 on 37059 degrees of freedom
#Residual deviance: 15002 on 37054 degrees of freedom
#AIC: 26533　　←複数のモデルを比較するときに使えるAIC　小さい方がよい｡

2)ここで推定したモデル
　以下､線量部分の定式化を紹介｡スクリプト内でのモデル番号､(AIC)､係数とその有意水準を示す(log(1+col_d10))などとあるのがそう｡

　モデル1-8ではベースラインには上のように､都市､被曝時年齢､到達年齢を､線量部分にはそれらを入れないシンプルなモデルを推定｡モデル9でLSS13と同じ程度の変数を投入したモデルを推定｡
　推定したモデル群のなかでAICが最小となるのはモデル9｡

より一般的なモデル

(1+x)^β
β=1のとき線形､β<1のとき下にそった形､β>1のとき上に反った形になる｡しかも1次+2次モデルのような多重共線性の問題が生じない｡

また､普通のポアソン回帰推定ルーチンでlog(1+x)を説明変数とすればよいので､指定や推定が楽｡
線量範囲を限定した推定､閾値を手動で変化させた推定もこのモデルで行っている｡
下記の線形(1次)と係数は異なるが､図示するとほぼ同じような形状になる｡　lssinc13forpubprc.R　の出力イメージ参照｡

モデル1　シンプルモデル(26533)

log(1 + col_d10) 0.6253448 ***

モデル5　同じことを自作の最尤ルーチンで(26533)　同じ結果

log(1 + col_d10) 0.625***

モデル9-0　LSS13のように変数追加｡それらは有意に説明するので､AICも(25818)と低下｡

log(1 + x) 0.474301 ***

上のモデル5にパラメータを1つ追加

(1+β1x)^β2

モデル　5.1　　( 26534.05)

β2:log(1 + col_d10) 1.0197　***
β1=b7*col_d10 0.5088　.

線量カテゴリ･ダミー(以下ノンパラ回帰と呼ぶ)

Σβi*xi 　　　xi　はi番目の線量カテゴリダミー
カテゴリの区切り方を変えて推定｡
モデル2　　LSS13と同じ23カテゴリ　(26556)
モデル2.1　0.01Svごとに細分化　　　(26799)
モデル2.2　1Sv以下では0.01区切り､それより大きいと0.1区切(26631)　
モデル2.3　23カテゴリを二つずつ併合､12カテゴリに｡(26546)

(手動)閾値モデル

β(x-T) x>=T ,0 otherwise

のように閾値以下を0にするのではなく上の一般モデルでの(手動)閾値モデルを推定｡

(1+x)^β2　 x>=T, (1+x)^β1　otherwise

ある(閾)値で傾きが異なるというモデル｡Tでの連続性無視｡

モデル3　　T=100msV　　　(26535)
log(1 + col_d10):col_d10<0.1FALSE 0.6269 ***

100mSV以上の傾き

log(1 + col_d10):col_d10<0.1TRUE　0.743(p=0.165)

100mSV未満　　以下同様

モデル3.1　T=10mSv　　　(26535)

log(1 +col_d10):col_d10<0.01FALSE 0.634　***
log(1 +col_d10):col_d10<0.01TRUE 11.750.

10mSV以下のリスクが非常に大きいという結果｡ただしあてはまりはモデル1よりも悪い｡特に理論的ななにかがなければ､モデル1を優先すべき｡

モデル3.2　T=5mSv　　　(26533)

log(1+col_d10):col_d10<0.005FALSE 0.6116***
log(1+col_d10):col_d10<0.005TRUE -47.47(p=0.112)

マイナスになる｡｡｡有意ではないが｡｡上と同様モデルあよりも劣る｡

モデル3.3　T=1mSv　　　(26535)

log(1+col_d10):col_d10<=0.001FALSE 0.627***
log(1+col_d10):col_d10<=0.001TRUE -17.35

線形(1次)

モデル6(シンプルモデル)　　(26532.98)

β　　0.5396　***

モデル9(変数追加=ほぼLSS13と同じ定式化｡結果も下記のようにLSSの0.475と一致)　　(AICも推定したモデルのなかで最小　25814.55)

β　　0.470　***

2次

βx^2

これは推定せず上記の一般的なモデルの枠組みで下記のように定式化した｡

(1+x^2)^β

モデル4.1　　　(26541)

β　→log(1 + col_d10sq) 0.4895710 ***

1次+2次

β1x+β2x^2　　　　　　　　　　推定せず
(1+x^2)^β1+(1+x^2)^β2

モデル4　#2乗項は有意ではない　(26533)

1次項､2次項の相関は0.95と高いので不適切だが参考までに｡
β1→log(1 + col_d10) 0.4348747 **
β2→log(1 + col_d10sq) 0.1651926 (p= 0.127)

(以下追加したもの)ロジスティック

1/{1+exp-(β2x-β1)}

β2の値が大きく､急に立ち上がるような形状であれば､β1/β2を閾値として解釈可能｡
モデル7　　　(26535)
logist const -3.908***
logist d 　4.584***

x/{1+exp-(β2x-β1)}

モデル7.1
分子にも線量を入れたもの｡(推定不能)

閾値直接推定

閾値を手動ではなく推定させる｡→統計的検定可能｡
β(x-T) x>=T ,0
モデル8　(26535)

slope 0.53134 ***
thresthold 0.02691 ***　←26.9mSv｡有意だがAICはモデル1よりも悪化

　ということで､AIC最小のLSS13の線形モデル9が最良｡つまり固形ガンによる死亡については､線形閾値無しモデルが､もっともあてはまりがよい｡というのが一連の分析からの結論｡

4.おまけの推定

1)線量域を限定して推定

上記の一般(シンプル)モデル1で､推定に用いるサンプルを線量によって限定｡この場合､サンプル数が異なるのでAICの比較は無意味｡線量の結果のみ示す｡

簡単のため変数の少ないモデル1の結果を紹介しますが､LSS13のモデル9でしても同様の傾向です｡ただし､低線量段階ではパラメータが推定できなくなる状況が多発します｡

10mSv以下

log(1 + col_d10) -3.727327 p=0.367
係数マイナスなのでホルミシス???ただし有意ではない

50mSv以下

log(1 + col_d10) 0.645　　　p=0.48673

100mSv以下

log(1 + col_d10) 0.621　　　p=0.25824

150mSv以下

log(1 + col_d10) 0.4599 p=0.21393

200mSv以下
log(1 + col_d10) 0.6317　　　　0.0282 *

ここでやっと有意になる｡

とまあ低線量では有意になりにくいことがわかった｡それでは高線量部分で限定するとどうか｡
2-3Svに限定

log(1 + col_d10) -0.097926 p=0.954　ここでも有意ではない｡｡｡｡

2Sv以上

log(1 + col_d10) 0.217170 p=0.878

1Sv以上だと有意にはなる

log(1 + col_d10) 0.569871 *

ということで､高線量部分でもサンプルを限定して推定すると有意ではなくなる｡

サンプル数が減少→検定力の低下
恣意的な限定
区間毎にパラメータが違うことを(暗黙のうちに)仮定しているようなもの
なので､よほどのことがなければ全サンプルを使うべき｡

参考までにsubjects数は下記の通り

sum(subjects[LSdat$col_d10<=0.05]) #47530
sum(subjects[LSdat$col_d10<=0.1]) #69108
sum(subjects[LSdat$col_d10<=0.15]) #72675
sum(subjects[LSdat$col_d10>=2 & LSdat$col_d10<=3]) #457
sum(subjects[LSdat$col_d10>=2 ]) #488
sum(subjects[LSdat$col_d10>=1 ]) #2101

2)推定時点を変更して推定

LSS14報の(非ガン)のグラフから刺激されて作成｡初期のころはU字型かつ､底の部分が1よりも小さいホルミシス的な形状だったが､時間が経過して死亡数が増加すると線形になる｡
このことをLSS13でも確認してみる｡1957年から10年ごとの死亡数を用いて､モデル1とモデル2で推定｡
下の図のように､初期の頃はモデル2(ノンパラ=○)で推定するとマイナスの係数が多発｡ホルミシス云々との根拠になったかもしれない

ただし､当時のデータでも､ノンパラモデルよりも線形モデルの方があてはまりは良好｡(例　1957年までの場合　モデル1のAIC=10874､モデル2=10899｡)

しかし､データが蓄積(死亡数が増加)されることによって､そのような減少はみられなくなる｡しかも線量の係数は大きくなる(傾きが急になる)｡
似たようなことは､台湾のコバルト混入セメントでつくられたアパートの分析でも観察されている｡
早い段階から安心するべきではないという教訓｡

注)1957(黒),67(青),77(黄),87(ピンク),97(赤)年時点までの死亡数を用いて推定した結果｡

　実線は[1+線量)^β､点は線量カテゴリダミー(23区分)を用いた｡

3)他の疾病

モデル1で白血病についても推定した例｡線量の係数は
log(1 + x) 1.833 ***
と1より大きいので上に反った形｡
固形ガンと白血病でこのように形状が異なるのはなぜか､とポーランドの学会で､放影研の児玉氏に質問した｡自分は､わからないが丹羽先生によると､固形ガンでも白血病のようになるべき､との見解とのこと｡
もう一人の放影研の研究員(米国で学位を取った台湾の方で白血病を分析)にも質問したが自分はstatisticianなのでメカニズムはわからない｡とのこと｡｡｡いずれにしてもERRは白血病の方が高い｡

5.推定の準備(上の0.と一部重複)

LSS13 データのダウンロード

放影研のここから住所などを入力すればすぐにダウンロード可能｡
Rで読み込むのは下記のファイル｡エディタなどで開いてみるとよい｡

R13mort.dat

R13 data documentation.pdfが変数名ファイル｡みておく方がよい｡

Rの準備

Rのインストール

まだの方は､例えば私のHPを参照してインストール(インストール方法､つかってみる　まで実行してみるとよいのでは)｡ライブラリなどは特に追加インストール不要｡

私の提供するファイル

処理用プログラム

下記リンクからGoogle docで公開｡GPLライセンス｡

lssinc13forpubrd.R　読み込み用Rスクリプト

lssinc13forpubprc.R　上でつくったデータセットを使って処理するための(推定用)Rスクリプト(出力画像付き)

6.Rへの読み込みとデータセットへの保存

Rを起動｡上記の　R13mort.dat　があるディレクトリをRの作業ディレクトリに指定(わかる人は当然､setwdコマンドで指定可能)｡

Mac版Rでは　その他>作業ディレクトリの変更
Win版Rでは　ファイル>作業ディレクトリの変更

Windows7ではデスクトップを指定できないかもしれない｡その場合はC:ドライブに　上記のファイルを移動して､C:を作業ディレクトリに指定してください｡

上記の　lssinc13forpubrd.R　読み込み用Rスクリプトの実行

上のリンクから画面を開いて順次､コピー､Rの画面にペースト､リターンすれば処理できます｡
(注)

#ではじまるのはコメント文なのでコピー､ペーストしてもしなくても構いません｡
なにかあるとエラーメッセージが表示されます｡このプログラムの場合､生じるとしたら､ファイルがみつからない､開けないというもの(日本語もしくは英語)｡その場合は､上記のように､R13mort.datのあるディレクトリを作業ディレクトリとしてください｡
Mac版SafariだとGoogledocからコピーするのにはコツがいるような感じ｡下の行から上にむかって行を選ぶと間違いが少ないよう｡

最後に下記のような出力がされ､フォルダに0LSdat.rda　ができていればok｡Rでのデータセット名はLSdat｡

　　被曝線量(行方向)別の人数subjects､死者数deathなど｡

subjects death cancer solid stomach lung liquid heart external
~0.005 37458 19182 4053 3833 1182 510 220 3577 1010
~0.1 31650 16132 3443 3277 1019 432 166 3020 774
~0.2 5732 3065 707 668 201 82 39 570 134
~0.5 6332 3446 812 763 256 112 49 663 156
~1.0 3299 1774 481 438 119 67 43 320 85
~2.0 1613 876 314 274 65 50 40 130 27
2.0+ 488 296 107 82 25 11 25 42 14
Total 86572 44771 9917 9335 2867 1264 582 8322 2200

7.　データの集計や推定

　上記のプログラムで作成した　0LSdat.rda　を読み込んで各種の推定を行う｡上と同様､0LSdat.rdaがある場所を作業フォルダとして指定｡

lssinc13forpubprc.R　上でつくったデータセットを使って処理するための(推定用)Rスクリプト(出力画像付き)
上記をクリックして開いた画面からコピー　Rにペースト｡

#のあとにコメントや処理結果も入れてあるので､処理できなくとも　どのようなことを(科学者は)しているのかを感じることができることを期待｡ちょっと勉強すれば､ここで推定したモデル以外も推定可能｡

8.演習方法

処理プログラムをどうしていいかわからない方は､とりあえずモデル1について

説明変数を追加する

といっても使えそうな変数はあと　distcat　被曝時の爆心地からの距離ダミー程度｡これを追加するには､下記のようにすればよい｡

カテゴリー変数なので､as.factor(distcat)とする｡

summary(res1.1<-glm(solid~ city+sex+agex+age+log(1+col_d10)+as.factor(distcat),family="poisson",offset=log(PYR),data=LSdat))

変数間の交互作用を追加する｡

例えば　被曝時年齢と性別の組み合わせで影響が違うと考えるならば､sex:agex　を加える｡もしくは読み込みプログラムで作成してある交互作用変数を利用する｡

summary(res1.2<-glm(solid~ city+sex+agex+age+log(1+col_d10)+sex:agex,family="poisson",offset=log(PYR),data=LSdat))
他の疾病について分析する｡

放影研からダウンロードした　変数名ファイル　R13 data documentation.pdfをみて変数名を指定｡　(Rは半角が基本｡また､大文字小文字を区別することに注意｡また､R13 data documentation.pdfでの変数名とデータファイル内ので変数名大文字小文字が異なっているものもある)｡
names(LSdat)

で下記のように変数名が出力されるので､それを利用｡

例えば胃がんの場合　死亡数は　stomach 胃の線量は　sto_d10

summary(res1.1<-glm(stomach ~ city+sex+agex+age+log(1+sto_d10),family="poisson",offset=log(PYR),data=LSdat))

胃の場合には胃の線量を使う方がよいと思いますが､未確認｡結腸線量のままでも結果は大きくはかわりません｡

部位別の線量の相関は高いものが多いことは下記のコマンドで確認できる｡
cor(LSdat[,70:88])

参考)変数名一覧

[1] "city" "sex" "un4gy" "distcat" "agxcat" "agecat" "dcat" "time"
[9] "PYR" "subjects" "gdist" "agex" "age" "year" "death" "cancer"
[17] "solid" "liquid" "oralca" "digestca" "esoph" "stomach" "colon" "rectum"
[25] "liver" "gallbldr" "pancr" "respca" "lung" "bone" "connect" "melanoma"
[33] "nmskin" "breast" "uterus" "cervix" "corpus" "utrNOS" "ovary" "othfem"
[41] "prost" "othmale" "bladder" "kidney" "othurin" "brainm" "brainb" "othcnsm"
[49] "othcnsb" "thyroid" "benign" "blood" "noncadis" "infect" "endo" "nervous"
[57] "heart" "resp" "digest" "urinary" "tb" "stroke" "pneu" "livcir"
[65] "external" "suicide" "kerma" "kerma_g" "kerma_n" "col_d10" "col_gam" "col_neu"
[73] "mar_d10" "mar_gam" "mar_neu" "bra_d10" "bre_d10" "liv_d10" "lun_d10" "ova_d10"
[81] "pan_d10" "ske_d10" "ski_d10" "sto_d10" "tes_d10" "thy_d10" "uri_d10" "ute_d10"
[89] "index" "lage70" "lage70sq" "lage70qsp" "lage40qsp" "e30" "e30sq" "e30sp"
[97] "e30qsp" "e50sp" "e50qsp" "msex" "female" "dsex" "fcity" "fsex"
[105] "mhiro" "mnaga" "fhiro" "fnaga" "mlage70" "mlage70sq" "mlage70qsp" "mlage40qsp"
[113] "me30" "me30sq" "me30sp" "me30qsp" "me50sp" "me50qsp" "flage70" "flage70sq"
[121] "flage70qsp" "flage40qsp" "fe30" "fe30sq" "fe30sp" "fe30qsp" "fe50sp" "fe50qsp"
[129] "col_d10sq" "d1" "d2" "d3" "d4" "d5" "d6" "d7"
[137] "d8" "d9" "d10" "d11" "d12" "d13" "d14" "d15"
[145] "d16" "d17" "d18" "d19" "d20" "d21" "d22" "dcat3"
[153] "agxcat3" "dcat4"
注)黒字はもともとの変数｡青は読み込みプログラムで作成したもの｡

9.期待すること

線形モデルのあてはまりが最良であることを実感する｡
部分的にサンプルを取り出して推定することの不適切さを実感する
論文や彼らの提供しているログの出力､ここでの出力をみて､それを上回る結論を提示できるようになる｡
ベースラインの結果が放影研のログと異なる原因の解明｡

性別2水準､都市2水準　これらをくみあわせて4つの変数を作成｡通常はどれかを0に固定しないと推定できないはずなので私は長崎･女性を0に固定｡しかし､放影研のAMFITからの出力では推定できている｡総合平均などからの乖離のようなものでしょうが､よく理解していない部分｡その他も若干異なるなので､Rでもそれができるようにしてもらえるとありがたい｡AMFITのマニュアルがあればわかるかもしれない｡

他の定式化､モデルで推定し､よりあてはまりのよいモデルを探索する｡

論文をみるとわかるように､分析はオーソドックス｡ちょっといじればLSS13モデルよりも優れたモデルがつくれるかもしれない｡

もともとはそれをめざして分析開始｡変数が少なければ､他のモデルの方があてはまることもありますが､LSS13のように変数を入れ込むとなかなか難しい｡

これよりも使いやすいスクリプトなどを誰かが開発してくれる｡

ここで紹介したoptimのプログラムは､変数等を変更すると､書き直さなければならず､面倒｡説明変数群行列xや線量反応部分を式として与えれば､推定してくれるように発展させるような奇特な方が出現する｡

より難易度の低い推定ツールの開発｡

エクセルの最適化ルーチンでも推定可能なはず｡市民が推定するのは困難なので､とっつきやすいものを｡

sourceforge.jpにプロジェクトを開設したが､CVS,SVNなど使ったことがないので､このように公開｡よりよい方法があれば｡

10.発展系

　下記のようなデータが存在するので､再分析などしてみる｡コホートデータならば多分､適用可能｡

放射線影響研究所　これを含む被爆者コホートデータ　住所など入力すればその場でDL可能

http://www.rerf.or.jp/library/dl/index.html

放射線医学総合研究所放射線安全研究情報DB

http://www.nirs.go.jp/db/anzendb/AnzenkenkyuDB.php
環境中の核種濃度､動物実験など｡　　オンラインで申請したらすぐにIDがメール返送された｡

STOREプロジェクト

http://www.storedb.org/
欧州における過去の放射線実験のデータのアーカイブ　誰でもダウンロード可能

･ドイツ　ウラン鉱山Wismut労働者のコホートデータ

http://www.bfs.de/en/bfs/forschung/Wismut/wismut.html
研究計画が認められればデータ提供も可能とある

米　Comprehensive Epidemiologic Data Resource

https://www.orau.gov/cedr/

ブラウザはIE/FFのみ対応?Safariではうまく作動せず｡
Hanfordなど米国核施設従業者の個人レベルデータ提供｡要書面でのユーザー登録｡公表前には了解を得る必要がある｡

11.参考

　いずれも私のこのブログの他エントリー｡

　疫学入門リンク集

　放射線影響研究所のデータについて
　低線量WG報告書の問題点

　原発従業者疫学調査について

12.謝辞
出版物ではないが､データの利用条件となっているので付加｡

用いたデータは広島および長崎の放射線影響研究所(放影研)から入手したものである。放影研は、日本の厚生労働省(厚労省)ならびに米国のエネルギー省(DOE) により資金提供を(後者については、その一部を米国学士院に対する DOE 研究助成金 DE-HS0000031 を通じて)受けている公益法人である。この報告書に示した結論は著者のものであり、必ずしも放影研またはその資金提供機関の判断を反映するものではない。

2012年5月6日日曜日