「国民の信頼醸成に向けた取組について（見解案）」に対する御意見を募集しています。 （１２月１２日（水）～１２月１８日（火））
　への意見　　12/16-17にかけて送信｡(800文字以内｡トピック分けの指定があるため)｡

1)総括の必要性

　原子力大綱以降に限定して､思いつくだけでも国民の信頼を裏切る行為として以下が挙げられる｡

･東京電力

　2002年に発覚したトラブル隠蔽

･その他電力会社など

　活断層の恣意的な解釈

　(国会事故調でも指摘された)土木学会津波審議への委員の多量選任に代表される規制への干渉

･電事連

　国会事故調に記載されているような､審議会などへの働きかけ

･各自治体､電力会社

　公聴会などにおける､いわゆるやらせ質問

･原子力委員会､福島健康調査など

　会合に先立つ秘密会

　まずはこれらについて記述｡その原因､反省は最低限必要｡これらを無視した記述､体質こそが信頼を損なわせた最大の原因である｡

2)対話､コミュニケーションの大前提としての意見反映プロセスの明示

　対話､コミュニケーションは､相互に譲歩する可能性があることが前提である｡

p.3｢国民参加の機会をこれまで以上に 整備するとともに、その結果、如何なる反映がなされたかについて丁寧に説明すべ きである。｣

　これまでは国や電力会社が立地などの結論を決めており､そのことを理解させることにのみ注力してきた｡信頼､コミュニケーションといっても､結局パブリックアクセプタンス(公衆に受容させる)ためのイベントであり､規定の路線を変えることはなく､文書の文言修正に取り入れる程度であった｡

　国民の意見をきいて根本的に国の方針を変更させることもあるような仕組みにすべきである｡そのためには､国民投票､住民投票などによる立地の意思決定手続きを明確化すべきである｡

3)政府のエネルギー・環境会議の結論の正しい理解を

p.1　｢政府のエネルギー・環境会議による国民的議論において、国民の多くが脱原子力 発電依存を望んでいると結論された背景には、こうした不信があるものと考えられ る。｣

→コスト(廃炉や補償は含まない)を比較しても原子力に優位性がないこと､(大きな議論にはならなかったが)核燃料サイクルが破綻していることなどにいての理解が深まったことが大きい｡信頼以前に原子力はだめであるという理解が深まったと考えるべきである｡

同　｢原子力発電に依存しない社会をできるだけ早期に実現することを目指す｣

　このことを最優先にすべきであり､｢安 全性が確認された原子力発電所を重要電源として活用すること｣　は副次的な手段である｡

　原発なしでもやっていけることはここ1年半で明確になった｡火力発電の効率改善､再生エネルギーの活用など､原発抜きの撤退戦略を原子力委員会では描くべきである｡

4)原子力委員会が現在行うべき(だった)こと
　この1ヶ月で人材､研究計画､信頼確保ついての見解をまとめている｡12月で任期が切れ､またそのあり方自体が議論されている原子力委員(会)が､将来を規定するようなことを提言すること自体が無意味である｡
　まずは､これまでの原子力政策､研究の課題をまとめるべきである｡具体的には発刊されていない原子力白書2011の中でこれを行い､最後の役目を果たすべきである｡

原子力委員会　「原子力人材の確保・育成に関する取組の推進について（見解案）」に対する御意見の募集　　　2012/11/15 送信


･海外について
　20年後には原発ゼロが原則である我が国が、海外において原子力を推進した
り、そのための人材を育成、支援する必要はない。

・国内について
　東大原子力を例にとると80年代中盤には既に、進学振り分けで底抜け状態=希
望すれば誰でも入れる低レベル学科となっていた。むろん、なかには優秀なもの
もいたが。
　人の数だけでなく質についても、過去30年近く原子力工学の質は低下し続け、
量子システムなどに名称変更されてきた。このことは当時、東大原子力の教授で
あった近藤委員長がよくご存じであろう。
　実際、特許の数や人員数をまとめてみても、原子力のピークは80年代初頭で
あった。
http://nonuke2011.blogspot.jp/2011/09/80.html
　技術的にも魅力のない分野なのである。除染技術の高度化などは、せいぜい表
面を洗い流す、削るといったローテクで行えばいい。
　若い者をこのような魅力のない分野に引きずり込むような方策はとるべきで
はないし、ここ10年程度、すでに原子力は、これらの策をとり、大して効果も上げ
ていないことを自覚すべきである。

・必要な人材のタイプについて
　これを読むと大学、研究レベルを念頭においているようだが、現場の作業員が
もっとも必要となるであろう。それらの人材について言及すべきである。
・倫理の重要性
　倫理についての言及は一カ所しかない。国会事故調が指摘したように、東電や
電事連によるロビー活動、その結果としての甘い基準、運用がなされているのが
現状である。
　つまらない技術的教育よりも倫理性を強化した教育をすべきである。

原子力災害対策特別措置法施行令の一部改正（案）に対する意見募集について へのパブリックコメント　8/12送信

１）対象について

　商業炉のみに限定する必要性がない。再処理設備なども含めるべきである。

２）距離、範囲について

　避難、屋内退避、安定ヨウ素剤の予防服用等の緊急防護措置を準備する区域（以下「緊急防護措置準備区域」という。）

であれば、今回の福島原発の例からみても飯舘村は５０ｋｍ圏内、福島市は８０ｋｍ圏内に入る。

よって、８０ｋｍ圏内とすべきである。

３）二つの条件ともにではなく、いずれかにすべきである。

　様々な要因によって、地域防災計画を策定できない市町村も存在する。それらが自動的に含まれるようにすべきである。

４）７条について

・住民の意思の尊重

　都道府県知事への協議のみとしているが、被害を受ける可能性があるのは住民である。住民の合意を得るべきである。

これを踏まえて、７条については以下のように書き換えるべきである。

２　原子力事業者は、前項の規定により原子力事業者防災業務計画を作成し、又は修正しようとするときは、政令で定めるところにより、あらかじめ、原子力設備（商業用発電所、再処理工場などを含む）から８０ｋｍ以内もしくは原子力防災計画を立案した市区町村内の住民の合意を得なければならない。

エネルギー･環境会議へのパブリックコメント　内容編　　　2012/7/15送信

1)不適切な前提

･電力需要量の過大な見積もり

　2030年の人口は10%程度減少が見込まれている｡それにも係わらず､｢2030 年までにGDP が２割以上増える｣といった実現性のない目標を前提としている｡

･一方で既存技術の改良､向上などの軽視

　化石燃料系の発電効率の改善､コジェネの活用などの要素が充分に取り入れられていない｡

･熱､電気の有効活用

　原発の発電効率は30%程度｡熱の7割は廃棄されている｡一方､ガスコンバインドでは70%程度の効率｡かつ都市部にも立地できるため､廃熱もさらに利用が可能である｡このようなトータルでみたエネルギー利用を考慮していない｡

　上述のようなGDPに挑戦的な目標を設定するならば､これら技術効率の改善にむけた目標の設定､インセンティブをもたせる制度設計をすべきである｡

2)3つのシナリオの非現実性

･25%という方向性を誤った選択肢

　0%　15%　25%を挙げているが事故前の原発依存度は26%である｡脱原発依存の方向性をめざすのであるから､これよりも低い目標となるのが当然であり､25%のシナリオはあり得ない｡

･原発設備の現状からみた非現実性

　今後､新規増設せず､40年で廃炉を前提とすると､2030年時点で操業40年未満となるもの､つまり､1990年以降に営業運転が開始されたのは､19基､合計2111万kWである｡

　ただし､これらには東海地震の影響が懸念される浜岡4-5号､中越地震で大きな被害を受けた柏崎刈羽2~7号機が含まれている｡これらは即時廃炉にすべきであるから､差し引くと1160万kWとなる｡さらに､これはあくまで設備容量であり､点検や不正によって停止している期間を除外した稼働率は70%程度とすると812万kwとなる｡これが2030年時点で可能な最大限の原発による発電量である｡

　設備容量1160万kWは1978年ごろ､参考までに2111万kWは1985年頃の水準である｡ただし､現在も大飯1基のみ､その一方で8基の火力発電所を停止した｡現在ただちに0基にしても発電能力としては問題がない｡

･再処理と組みあわせることの不適切性

　再処理､高速増殖炉など､40年取り組んできたが､前者はガラス溶融段階での煉瓦の混入､後者はナトリウム漏れ､燃料棒の落下など極めて初歩的な段階にとどまっている｡米国､フランスなどは90年代にはこのような方策はあきらめた｡再処理を前提とした核燃料サイクルは既に実現不可能であることを前提にすべきである｡

3)偏った情報提供

･核廃棄物の発生の無視とCO2についての偏重した記述

　原子力発電所の運転によって､高レベル､低レベル核廃棄物が発生する｡使用済み燃料プールに13,530

トンU が存在する(電事連　「原子力・エネルギー図面集　7−28 各原子力発電所の使用済み燃料の貯蔵量」　http://www.fepc.or.jp/library/publication/pamphlet/nuclear/zumenshu/digital/index.html ) 

運転すると､処理技術も確立していない廃棄物が増加するわけである｡このことを明示すべきである｡

･間接的オフセット取引の無視

　CO2の削減については､途上国への省エネ､環境技術の供与を提供する間接的なオフセットも考慮される｡このことがまったく呈示されていない｡自国で上述のように効率を改善するだけでなく､他国への支援を進めることによって充分に国際的に貢献できる｡

4)選択肢の呈示方法

　上述のように%表示では､実際に必要な電力量といった重大な情報が隠されてしまう｡%ではなく電力もしくは原発の基数(120万kW換算)すべきである｡

エネルギー･環境会議へのパブリックコメント　手順編　　2012/7/15送信

1)不適切な選択肢の設定

　3つの選択肢を示しているが､脱原発の方向といいながら現状の26%と同様の25%という選択肢を設定している｡また､これとあわせて､再処理の方法も組み合わせている｡このような不適切な選択肢を呈示するべきではない｡

　そもそも%で表示するのではなく､原発による発電設備量(例　0基､　5基:700万kW程度:　10基　)という､より具体的でわかりやすい方法で明示すべきである｡

　なお､さいたま説明会での高原長官の説明によると2010年策定エネルギー長期計画での目標値45%と比較しているようだが､2010年実績の26%と比較すべきである｡

2)一方的な情報の提供

　CO2の問題を強調しているが､一方で原発の運転に伴う核廃棄物の発生｡再処理や高速増殖炉などの技術的行き詰まりなどの情報が提供されていない｡3つの選択肢について､再処理の方法も抱き合わせている以上､これらの情報も提供すべきである｡

3)期間の短さ

　6月末の方針決定から7月末までのパブリックコメント(8月上旬まで延長)､8月上旬までの意見聴取会､討論型世論調査といった手順を想定しているようである｡このような重大な問題を､2ヶ月未満という短期間で決定すべきではない｡

　選択肢の設定､それらへの評価､判断などじっくりと時間をかけるべきである｡

4)参加可能な人数の少なさ

　上述の手順では､パブリックコメントは自由に誰でも寄せられるものの､意見聴取会200名×11会場､討論型世論調査では3000名にアンケート､うち300名程度しか参加できない｡これでは参加可能な者の数が圧倒的に少なく､国民的議論とはいえない｡

　これらの情報はいずれも､ホームページ上で公開されているだけであり､広く告知されていない｡がれき受け入れ広告と同等程度の広報を行うべきである｡

　実際にさいたまでの意見聴取会に出席したが､政府担当者や意見表明者への質問すらさせていない｡結局政府担当者の説明をきき､意見表明者の意見をきいているだけで､自宅でネット中継をみているのと同じである｡

5)民意の反映の方法の不明確さ

　上述のように､いくつかの方法で意見を聴取するようだが､それらをどのように判断するのかの基準､方向性すら示されていない｡パブリックコメントは募集しても､｢取り入れることはできない｣｢今後の課題｣といった形で無視されてきたものがおおい｡期間の短さもあり､さらにその危惧は大きい｡

　国民的な論議を取り入れるならば､国民投票もしくは少なくとも立地･消費地での住民投票を行うべきである｡

原子力委員会への意見　　2012/5/26　送信


1)事務局の構成の明示､
「原子力委員会及び原子力安全委員会設置法」 （PDF:144KB）
「原子力委員会及び原子力安全委員会設置法施行令」
http://www.aec.go.jp/jicst/NC/about/index.htm

をみると､安全委員会については事務局についての項目があるが､原子力委員会については存在しない｡事務局の構成などについて明示すべきである｡

2)原子力委員会議事運営規則
http://www.aec.go.jp/jicst/NC/about/hourei/4.pdf
によれば､
｢第六条 委員会の議事録は、速記録として作成し、発言者の確認を経て、原則次の回の定例会議又は臨時会議において配布するものとする。｣　とあるが､5/8以降の議事録が未だに公開されていない｡公開という意味でも原則に立ち返り､毎回､即座に公開すべきである｡

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次)
• βx
• 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次)
• βx
• モデル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 を通じて)受けている公益法人である。この報告書に示した結論 は著者のものであり、必ずしも放影研またはその資金提供機関の判断を反映するも のではない。