グラフの例:トンデル論文

はじめに

トンデル(Martin Tondel)は放射線の低線量被ばくの影響を研究している人です。2011年12月28日のNHKの番組「追跡!真相ファイル」低線量被ばく 揺らぐ国際基準にも登場しました。しかし,彼の論文についてはいろいろ批判があります。

トンデルの2004年の論文

まず Increase of regional total cancer incidence in north Sweden due to the Chernobyl accident? という論文です。表1のデータをCSVの形で載せておきます。これはスウェーデン北部の地域を放射能汚染度で分類したときの人口・発がん者数・総死亡数です。1列目がセシウム137(137Cs)の密度(kBq/m2),次の3列が1988年の人口(男・女・計),次の3列が1988-1996年の発がん数(男・女・計),最後の3列が1988-1996年の総死亡数(男・女・計)です。

Cs137,Pop_M,Pop_W,Pop_T,Can_M,Can_W,Can_T,Dea_M,Dea_W,Dea_T
<3,185646,173863,359509,2810,3881,6691,5629,3029,8658
3–29,269922,257890,527812,4501,5877,10378,7294,4358,11652
30–39,48053,44270,92323,752,1075,1827,1294,811,2105
40–59,63512,61350,124862,1185,1559,2744,2048,1188,3236
60–79,11112,10513,21625,172,229,401,291,172,463
80–120,8722,8329,17051,153,215,368,288,159,447
Total,586967,556215,1143182,9573,12836,22409,16844,9717,26561

これをコピーして次のどちらかをRに打ち込めば tondel1 というオブジェクトにデータが読み込まれます。

tondel1 = read.csv("clipboard")     # Windowsの場合
tondel1 = read.csv(pipe("pbpaste")) # Macの場合

まずは発がん率を汚染度ごとに調べ,グラフにしてみます。エラーバーは95%信頼区間です。

# quartz()                      # Mac
# par(family="HiraKakuProN-W3") # Macフォント
par(mgp=c(2,0.8,0))             # 好み:マージンを狭く
attach(tondel1)
ci = matrix(NA, 6, 2)  # 6×2行列を作る。初期値はすべてNA
for (i in 1:6) ci[i,1:2] = binom.test(Can_T[i], Pop_T[i])$conf.int
plot(c(0.8,6.2), range(ci), type="n", xaxt="n",
     xlab=expression(paste({}^137, Cs, " (", kBq/m^2, ")")),
     ylab="Cancer Rate")
points(1:6, Can_T[1:6]/Pop_T[1:6], type="o", pch=16, lwd=2)
arrows(1:6, ci[,1], 1:6, ci[,2], angle=90, code=3, length=0.05)
axis(1, 1:6, Cs137[1:6])

上の CanDea に置き換えて死亡数についても調べてみましたが,微妙です。それぞれの地域性(年齢構成,都会か田舎か,煙草の消費量など)が考慮されていないからです。

発がん数と汚染度 死亡数と汚染度

論文では表5に5-59歳の発がん数の素データが1986-1987年と1988-1996年について載っています。CSV形式で表すと次のようになります:

Cs137,Can_86_87,Can_88_96
<3,712,4181
3–29,1024,6402
30–39,175,1130
40–59,254,1623
60–79,36,252
80–120,32,235
Total,2233,13823

実際の論文では,年齢調整などいろいろな統計処理をしているのですが,残念ながらそれを検証するデータを私は持ち合わせていません。しかし,この素データは,同じ地域についての比較ですので,地域性による交絡因子には影響を受けにくいと考えられます。そこで,このデータだけに基づいて,できるだけモデルに依存しないように,1986-1996年の発がん数のうちどれだけの割合が1988-1996年のものであるかを,2項分布を仮定して binom.test() で調べた95%信頼区間とともにグラフにしてみます。1986年のチェルノブイリ原発事故が起きたばかりの1986-1987年と,それによる影響が出始めた1988-1996年との発がん数の比は,その地域の汚染度によって違いがあるだろうというわけです。

発がん数,1988-1996年/1986-1996年

このグラフを描くには,上のデータをクリップボードにコピーしてから,次のように打ち込みます:

tondel2 = read.csv("clipboard")     # Windowsの場合
tondel2 = read.csv(pipe("pbpaste")) # Macの場合
attach(tondel2)
for (i in 1:6)
  ci[i,1:2] = binom.test(Can_88_96[i], Can_86_87[i]+Can_88_96[i])$conf.int
plot(c(0.8,6.2), range(ci), type="n", xaxt="n",
     xlab=expression(paste({}^137, Cs, " (", kBq/m^2, ")")),
     ylab="Cancer 88-96/86-96")
points(1:6, Can_88_96[1:6]/(Can_86_87[1:6]+Can_88_96[1:6]),
       type="o", pch=16, lwd=2)
arrows(1:6, ci[,1], 1:6, ci[,2], angle=90, code=3, length=0.05)
axis(1, 1:6, Cs137[1:6])

なんとなく右上がりの傾向が見えてきました。このあとはモデルに依存しますし,計算もおおまかですが,横軸を等間隔または区間の中央と仮定し,2項分布の一般化線形回帰をしてみます。

n = Can_86_87[1:6] + Can_88_96[1:6]
p = Can_88_96[1:6] / n
x = 1:6  # あるいは区間の中央をとって次のようにする
# x = c(3/2, (3+29)/2, (30+39)/2, (40+59)/2, (60+79)/2, (80+120)/2)
fit = glm(cbind(Can_88_96[1:6], Can_86_87[1:6]) ~ x, family=binomial)
summary(fit)

結果のサマリー summary(fit) を見ると,傾きの p 値は,等間隔で p = 0.0855,区間の中央で p = 0.0877 程度です。微妙に有意ではありませんでした。

しかし,こういうものには交絡がつきものです。トンデルの分析では(年齢調整などのほか)人口密度を交絡因子として考慮しているのが特徴です(スウェーデンでは人口密度40人/km2未満の地域ではがんが少ないようです)。その結果,ぎりぎり2σレベルで有意だというのが彼の結論です。もっとも,2σレベルではそれほど強い証拠と言えませんし,隠れた交絡があるかもしれませんので,この論文のタイトルには「?」が付けられています。

トンデルの2011年の論文

Risk of malignancies in relation to terrestrial gamma radiation in a Swedish population cohort という論文は,被ばく線量とがんの関係を調べています。

次は論文の表2から取ったものです。TGR(terrestrial gamma radiation)が線量率(nGy/h,ほぼ0.001μSv/h)で,HR(hazard ratio)がハザード比(0-60を1としたときの危険性の度合),その右の二つが95%信頼区間です。グラフ化は宿題です。

TGR,HR,CI_L,CI_H
0-60,1,1,1
61-75,1.04,1.02,1.06
76-83,1.05,1.02,1.07
84-91,1.07,1.04,1.10
92-95,1.04,1.10,1.08
96-366,0.99,0.96,1.02

参考


Last modified: 2012-01-01 07:20:23