So-net無料ブログ作成

エクセルでやったらどえらい簡単にできたが間違ってないか? [統計学もどき]

EZRの記事がコメントが増えたのと、今月も1つくらいは書いておかねばと思いましたので、新たに設けたカテゴリーで行ってみたいと思います。

アワプラさんが貴重なまとめをしていただいています。

甲状腺がん悪性・悪性疑い152人〜福島県民健康調査
http://www.ourplanet-tv.org/?q=node/2004

検討委員会での記者会見にて、基本調査の①-8ページの福島県住民の外部被ばく線量分布と、甲状腺検査(本格検査)の②-6ページの甲状腺がん症例の外部被ばく線量分布が大きく異なり、甲状腺症例の外部被ばく線量が高い方に分布しているという指摘があった。このオッズ比に関して、性別を調整して表2として示した。この調整オッズ比とその95%信頼区間の推定結果は、EpiInfo 7TMのStatCalcの最尤推定値(Mid-P)を用いた。95%信頼区間の下の値が1より大きいと、いわゆる「統計的有意差がある」ということになる。
  その結果、明瞭な量反応関係が観察された。つまり、外部被ばく線量に関連する何らかの要因(恐らく内部被ばく)が甲状腺がんの発症に関与していることが、より明瞭に推察できることになった。1mSv未満もまた、被ばくしていると考えられるので、示されたオッズ比は過小評価されている。先行検査では、平成23年度、24年度、25年度の別に、甲状腺がん症例の線量分布が発表されておらず、平均有病期間が平成26年度としてほぼそろっている本格検査とは異なり、明瞭ではなかった。基本調査の住民の性・年齢・市町村別線量分布、先行検査および本格検査における甲状腺がん症例の性・年齢・市町村別線量分布の発表が待たれるところである。
2015113014.jpg 
これだけでも十分なのですが、ちと統計学の相関関係が気になってあれこれ探して、Rで試すのはなんだかまだ私には無理そうだったので、エクセルのはないかと見たところ

散布図と相関係数(相関分析)
http://hitorimarketing.net/tools/correlation-analysis.html

これに倣ってやってみたところ(発症者数は100万人当たりに換算してます。)

相関係数.png

案外 楽にピアソンの積率相関係数(単相関係数)とやらが表示されました。簡単すぎたのでこれでいいのかと思ってしまいましたが、点が三つしかない散布図でも、正の相関になっていますので、多分これでいいのでしょう。Rに慣れていけば、P値なんかもわかるので、相関係数の正確さもわかってくるでしょう。

毎月3日のアベ政治を忘れておった。流行語大賞ノミネートまでいったんでした。
sho_fj.jpg
nice!(9)  コメント(25)  トラックバック(0) 
共通テーマ:blog

nice! 9

コメント 25

ryu-ron

やっぱRでやるのは難しいですね。
エクセルのCSVデータで読み込めたけど、空白のデータをnulかなにかに書き換えないといけないのでしょうが、なぜか編集でエラーがでます。
まあゆっくりやります。
by ryu-ron (2015-12-04 09:06) 

ryu-ron

なんとかなったような。

> x <- read.csv("soukan3.csv", header=TRUE)

> x
mSv 人数
1 0 NA
2 1 28
3 2 75
4 3 NA
5 4 NA
6 5 146

> cor(x)
mSv 人数
mSv 1 NA
人数 NA 1

> cor(x, use="complete.obs")
mSv 人数
mSv 1.0000000 0.9865587
人数 0.9865587 1.0000000

> cor.test(x $ mSv, x $ 人数 )

Pearson's product-moment correlation

data: x$mSv and x$人数
t = 6.0374, df = 1, p-value = 0.1045
alternative hypothesis: true correlation is not equal to 0
sample estimates:
cor
0.9865587

これが積率相関係数でエクセルのとほぼ同じでした。残念ながらP値は5%を大幅に上回りました。

> cor.test(x $ mSv, x $ 人数, method="s")

Spearman's rank correlation rho

data: x$mSv and x$人数
S = 0, p-value = 0.3333
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
1

こちらは順位相関係数とやらであんまり意味がわかってませんが、係数は1ですか。うーん わからん。
1mSvごとにデータがでないと、正確なのは難しそうです。
by ryu-ron (2015-12-04 11:49) 

ryu-ron

今日「確率のはなし」が来るかと思ったら、県立の方がきのうまで休刊日で、1週間伸びてしまいました。
うー 来週以降は、朝の目覚めも悪そうで本を読む元気が残っておるかなあ。

そういえば
t = 6.0374, df = 1, p-value = 0.1045

これはt値で自由度1ですねえ。データそのものは3あるので、3-1で2じゃないのか?なんだかいろいろ難しい。

alternative hypothesis: true correlation is not equal to 0
sample estimates:

これを訳してみないといけないなあ。
by ryu-ron (2015-12-04 12:30) 

ryu-ron

> x <- read.csv("soukan3.csv", header=TRUE)

2mSvから5mSVまでを1mごとにして、数値を割当たら傾きが変わるのでこれだと相関関係が弱いですねえ。

> x
mSv 人数
1 0 7
2 1 21
3 2 75
4 3 45
5 4 45
6 5 56

> cor.test(x $ mSv, x $ 人数 )

Pearson's product-moment correlation

data: x$mSv and x$人数
t = 1.6212, df = 4, p-value = 0.1803
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3719271 0.9538215
sample estimates:
cor
0.6297055

ただ今度は95%信頼区間が出てもっともらしくはなります。
by ryu-ron (2015-12-04 13:08) 

ryu-ron

そうか点が少なくても、山本医師の
http://ryukei-rondan.blog.so-net.ne.jp/2015-11-12

オッズ比をとった単純回帰分析
http://ryukei-rondan.blog.so-net.ne.jp/upload/detail/m_image003.jpg.html

これを今回の数値でやれば精度を上げられそうですね。
by ryu-ron (2015-12-04 14:25) 

ryu-ron

弁当屋さんサイトを参考に回帰分析に挑戦したのだけど
http://rmecab.jp/ranko/chap3.html

> x <- read.csv("soukan3.csv", header=TRUE)

> x
mSv 人数
1 0 NA
2 1 28
3 2 75
4 3 NA
5 4 NA
6 5 146

> soukan3.lm <- lm (mSv ~ 人数, data = soukan3)

> x.lm <- lm (mSv ~ 人数, data = x)

> summary(x.lm)

Call:
lm(formula = mSv ~ 人数, data = x)

Residuals:
2 3 6
0.2347 -0.3901 0.1554

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.202701 0.550467 -0.368 0.775
人数 0.034571 0.005726 6.037 0.104

Residual standard error: 0.4811 on 1 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.9733, Adjusted R-squared: 0.9466
F-statistic: 36.45 on 1 and 1 DF, p-value: 0.1045

error: 0.4811 on 1 degrees of freedom
間違ってたのか?
by ryu-ron (2015-12-04 15:47) 

ryu-ron

これがロジスティック回帰分析?

> x <- read.csv("soukan3.csv", header=TRUE)

> x
mSv 人数
1 0 NA
2 1 28
3 2 75
4 3 NA
5 4 NA
6 5 146

> x.glm <- glm (mSv ~ 人数, data = x)

> summary(x.glm)

Call:
glm(formula = mSv ~ 人数, data = x)

Deviance Residuals:
1 2 3 4 5 6
0.2347 -0.3901 0.1554

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.202701 0.550467 -0.368 0.775
人数 0.034571 0.005726 6.037 0.104

(Dispersion parameter for gaussian family taken to be 0.2314159)

Null deviance: 8.66667 on 2 degrees of freedom
Residual deviance: 0.23142 on 1 degrees of freedom
(3 observations deleted due to missingness)
AIC: 6.8272

Number of Fisher Scoring iterations: 2

lmをglmに変えただけなんだけど。
by ryu-ron (2015-12-04 16:31) 

ryu-ron

合計だけ合わせて、空白のないデータを作ってロジスティック回帰分析をするといろいろと数字がでてきますね。

> x <- read.csv("soukan4.csv", header=TRUE)

> x
mSv 人数
1 0 8
2 1 16
3 2 30
4 3 45
5 4 60
6 5 90

> x.glm <- glm (mSv ~ 人数, data = x)

> summary(x.glm)

Call:
glm(formula = mSv ~ 人数, data = x)

Deviance Residuals:
1 2 3 4 5 6
-0.47685 0.04001 0.19451 0.28863 0.38274 -0.42904

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.006289 0.303177 -0.021 0.984445
人数 0.060392 0.006075 9.942 0.000575 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 0.1701724)

Null deviance: 17.50000 on 5 degrees of freedom
Residual deviance: 0.68069 on 4 degrees of freedom
AIC: 9.9688

Number of Fisher Scoring iterations: 2
by ryu-ron (2015-12-05 09:06) 

ryu-ron

なかなか難しい。空白の行をRの方でカットするやりかただとエラーが出てうまくいかず、予めエクセルの方でカットしたものを読みこんでロジ回(長いので次からこれでいきます。)すると
> x
mSv 人数
1 1 28
2 2 75
3 5 146

> x.glm <- glm (mSv ~ 人数, data = x)

> summary(x.glm)

Call:
glm(formula = mSv ~ 人数, data = x)

Deviance Residuals:
1 2 3
0.2347 -0.3901 0.1554

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.202701 0.550467 -0.368 0.775
人数 0.034571 0.005726 6.037 0.104

(Dispersion parameter for gaussian family taken to be 0.2314159)

Null deviance: 8.66667 on 2 degrees of freedom
Residual deviance: 0.23142 on 1 degrees of freedom
AIC: 6.8272

Number of Fisher Scoring iterations: 2

いろいろ勉強せんとわからんことが多いですね。
by ryu-ron (2015-12-05 09:33) 

ryu-ron

統計的にしっかり理論づけて書いておられる方がいます。
http://fiddledad.blog.fc2.com/blog-entry-67.html

わたしもこんな風に書けるようになりたいです。
by ryu-ron (2015-12-05 10:19) 

ryu-ron

ロジ回はこのアスタリスクに意味があったはずで
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ある程度のデータ数があったときに、影響のある要素のアスタリスクが多くなるはず。
このケースは外部被ばく線量のみなので、ほんとはいろいろ考えられる原因について、統計値が必要になる。

しばし中断でまた仕事モードに戻ります。
by ryu-ron (2015-12-05 10:52) 

ryu-ron

仕事モードのはずが、あとひとつ。
http://diamond.jp/articles/-/82512

広瀬さんとマコさんの対談もすごく大事です。
茨城北部は要注意です。
by ryu-ron (2015-12-05 11:31) 

ryu-ron

すぐ昼休みになってるが、気合が足りません。

「基本調査の住民の性・年齢・市町村別線量分布、先行検査および本格検査における甲状腺がん症例の性・年齢・市町村別線量分布の発表が待たれるところである。」

記事の引用文にありますが、結局 細かいデータを残してない。まあ被曝に原因を特化されるのを恐れてあいまいにしたいのでしょう。
by ryu-ron (2015-12-05 12:09) 

ryu-ron

ロジ回で、二項分布のときは、
x.glm <- glm (mSv ~ 人数,family=binomial,data = x)

と指定しないとけないとなってるのですが、

[5] エラー: y values must be 0 <= y <= 1

なにやら、データにyになるものを指定しないときけないみたいですね。
やっぱ難しい。
by ryu-ron (2015-12-06 11:25) 

ryu-ron

奥村教授の世話になりたっくないのだが
http://oku.edu.mie-u.ac.jp/~okumura/stat/logistic.html
ここに書いてることを理解するのが近道かもしれませんね。

それと大分が誇るいや誇った過去形かな、KOKIKOKIYAさんが、内部被ばく心配なしといいだした。

私はホットパーティクルをあきらめきれませんね。

名前はセシウムボールでもいいのですが、肺腺がんに特化して統計学的に追及していきたいと考えています。

by ryu-ron (2015-12-06 12:21) 

ryu-ron

よくわかんねえな。

奥村教授のやり方を真似てみたら

> x <- read.csv("soukan6.csv", header=TRUE)

> x
y mSv 人数
1 0 0 NA
2 1 1 28
3 1 2 75
4 1 3 NA
5 1 4 NA
6 1 5 146

> x.glm <- glm (y~mSv+人数, data = x,family=binomial(link="logit"))

> summary(x.glm)

Call:
glm(formula = y ~ mSv + 人数, family = binomial(link = "logit"),
data = x)

Deviance Residuals:
[1] 0 0 0

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.357e+01 9.690e+04 0 1
mSv -1.625e-08 1.652e+05 0 1
人数 3.412e-10 5.788e+03 0 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 0.0000e+00 on 2 degrees of freedom
Residual deviance: 3.4957e-10 on 0 degrees of freedom
(3 observations deleted due to missingness)
AIC: 6

Number of Fisher Scoring iterations: 22

> fitted(x.glm)
1 2 3 4 5 6
NA 1 1 NA NA 1

まあことさらエラーもなかったのですが、意味がわかっておりません。
by ryu-ron (2015-12-06 12:39) 

ryu-ron

詳細データがわからないままでロジ回をしても意味がないとは思ったのですが、奥村教授の
http://oku.edu.mie-u.ac.jp/~okumura/stat/logistic.html
カブトムシの問題の表は使えそうなので、
甲状腺がんのデータを次のようにして

data1<- read.csv("soukan7.csv", header=TRUE)
data1
result3 = glm(cbind(y,n-y) ~ x, data=data1, family=binomial(link="logit"))
summary(result3)


> data1<- read.csv("soukan7.csv", header=TRUE)

> data1
x n y
1 1 283286 8
2 2 145455 11
3 5 27327 4

xがmSv、nが住民線量分布、yが癌症例の線量分布です。
mSvは、上限の数値を取るので、不正確ではあるのですが、まあやってみようということで

> result3 = glm(cbind(y,n-y) ~ x, data=data1, family=binomial(link="logit"))

> summary(result3)

Call:
glm(formula = cbind(y, n - y) ~ x, family = binomial(link = "logit"),
data = data1)

Deviance Residuals:
1 2 3
-0.7617 1.0542 -0.3690

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -10.6087 0.3649 -29.074 < 2e-16 ***
x 0.3917 0.1377 2.844 0.00446 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 8.1610 on 2 degrees of freedom
Residual deviance: 1.8278 on 1 degrees of freedom
AIC: 17.282

Number of Fisher Scoring iterations: 4

よくはわからんのですが

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

が出たので、一応エラーなしということで、この意味をしらべて、ロジ回実験は取りあえず年内はこれが最後ということにしておきましょう。

ある程度のデータ数になると

Deviance Residuals:
Min 1Q Median 3Q Max
-1.5941 -0.3944 0.8329 1.2592 1.5940

というように最小、最大、中央値なんぞが出るみたいです。
by ryu-ron (2015-12-07 12:12) 

ryu-ron

カルトフェンさんの時にいろいろ教わったKokikokiyaさんですが、

http://pfx225.blog46.fc2.com/blog-entry-2734.html#comment3939

思いがけずに内部被曝は少なかったのに驚きと共に汚染地の子供はチェルノブイリの二の舞にはなりそうもないのが、
一安心です。
「見捨てられた初期被曝」は相変わらず甲状腺ガンが出ていますので初期被曝された近場の方は要注意ですが。
避難された方は内部被曝の打撃は今後ますますリスクが低下するので人生を前向きに生きるべきです。
ただ汚染地には相も変わらずガンマ線が飛び交っていますので、気を付けなければなりません。
とりあえずは内部被曝は一安心の5年目の客観的事実です。

http://pfx225.blog46.fc2.com/blog-entry-2731.html#comment3935

私が言っているのは、木下黄太よりもより恐ろしい厳しさなのです。

まあ、言わんとすることはわからんでもないし、原発容認派に転向したわけでないから、批判をしても何かが変わるものではないのですが、なんとなく読み方を誤れば推進派からみたら、内部分裂と悪宣伝に使われかねない内容ではあります。

ホットパーティクル、セシウムボールがある以上まだ安心できんのが神奈川以北の地域です。
by ryu-ron (2015-12-07 15:35) 

ryu-ron

R ポアソン分布 で検索してなかったから気がつきませんでしたがこんな記事がありました。
http://www.geek.sc/archives/679#i-2

Rでポアソン分布で、100万回に1回起こる可能性が8万回に1回起こる累積確率を求めても7%と、信頼度を95%とするならば、有意とは結論づけられない結果が得られました。

また、この甲状腺癌の問題に関しては、全国調査では潜在的患者の母集団が隠れているので、実際には全国でも甲状腺癌の潜在患者はかなりいて、発症していないか何らかの強制的な検査によって発見されていないだけで、福島の子供の確率とは比較できないとする、標本に関する問題提起もあります(http://togetter.com/li/241058)。

2013.03.12の記事ですか、ちと古いな。

by ryu-ron (2015-12-07 17:01) 

ryu-ron

あれ、WHOといえば一度コメントで紹介してましたね。
http://tkj.jp/takarajima/contents/blog/p/1016/

もともと幼少期の甲状腺ガン発症率は非常に低い。従って、幼少期に被曝した場合のリスクを、原発事故発生からの15年間に絞って計算すると「小児甲状腺ガンと被曝との関係性がより明白になる」と、WHO報告書は言う。
 ひょっとするとこの部分が、原発事故による健康被害は「ない」とする評価を続ける環境省や専門家会議の癇に障ったのかもしれない。
 多発が予測されたのはそれだけではない。
 小児甲状腺ガンほどではないにせよ、小児白血病も多発するという。被災時に1歳だった男児の場合、浪江町では事故発生からの15年間で発症率は1.8倍(同0.03%→同0.055%)に増え、飯舘村では15年間で1.5倍(同0.03%→同0.044%)に増える。1歳女児の場合、浪江町では事故発生からの15年間で発症率は1.6倍(同0.03%→同0.047%)に増え、飯舘村では15年間で1.3倍(同0.03%→同0.04%)に増える(同報告書62ページ。【図2】)。
 そして、乳ガンも増える。被災時に10歳だった女児の場合、浪江町では事故発生からの15年間で発症率は1.5倍(同0.01%→同0.015%)に増え、飯舘村では15年間で1.3倍(同0.01%→同0.013%)に増える(同報告書63ページ。【図3】)。

うーん 読み返す時間がないなあ。
by ryu-ron (2015-12-08 09:49) 

ryu-ron

あれ宝島 2014年12月14日 10:00 去年の記事か。
林先生のツイッターからの情報でした。
https://twitter.com/SciCom_hayashi/status/673993014797197312
by ryu-ron (2015-12-08 09:56) 

ryu-ron

日本語の白書が
http://www.unscear.org/docs/reports/2015/Fukushima_WP2015_web_jp.pdf
宝島が見たのは
http://apps.who.int/iris/bitstream/10665/78218/1/9789241505130_eng.pdf?ua=1
去年の分でしょうが、実質的には今年の白書の中身はこれがベースでしょう。
都合の悪い部分はカットして、日本語の白書にしたのかいな。

ムラにとっては、1.8倍とか1.5倍ならsクリーニング効果か、過剰診断で逃げられる数字だと思っているのでしょう。
by ryu-ron (2015-12-08 11:41) 

ryu-ron

もうすぐ確率のはなしがきますので、年内はポアソン分布に集中して、
http://cse.fra.affrc.go.jp/kiyo/home/etude/scRibbles/entori/2011/12/22_poason_fen_buwo_linishita_zui_you_fano_shuo_ming.html
このあたりを読みこんで、津田先生が分析に使った手法をRで再現しようなどと大それたことを考えていますが、実現できるかどうかは未知数です。たぶん挫折するでしょう。
うまくいかないと気は来年の宿題にします。
by ryu-ron (2015-12-08 15:07) 

ryu-ron

さきほどのページは色が薄く目に悪いので、テキストをワードにコピペして、いくらかEZRで試しましたが、イミフが多すぎて、越年決定。
まあゆっくりいきましょう。

われわれは、“Epi Info 7”ソフトのStatCalcを用いて、Pを中間値とする最大尤度オッズ比によって、PORの95% CI、ガイギー科学表に記載されているポアソン分布にもとづいてIRRを計算した23。

同じことはRでできそうな気はしてるのですが。
by ryu-ron (2015-12-08 16:25) 

ryu-ron

岩上さんも元気になったみたいですが、無理しないで下さいね。
津田先生のインタビューはどうなったかわかりませんが、甲状腺がんの話題も下火になりつつあります。
また来年年が明けて火をつけるといたしましょう。

久々にラジウムホットスポットの新ネタを探してみたけど、出てこない。
何か作戦を考えましょう。

あ 本の話題はPKOの方に書いておくか。
by ryu-ron (2015-12-09 08:51) 

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

※ブログオーナーが承認したコメントのみ表示されます。

トラックバック 0