2011年12月23日金曜日

survreg と weibreg

Rです。


survreg
> library(survival)
> library(MASS)
#library(MASS)はデータ"gehan"を使うため

> survreg(Surv(time,cens)~treat,dist="weibull", data=gehan) #デフォルトが"weibull"なので dist="weibull" は省略可
Call:
survreg(formula = Surv(time, cens) ~ treat, data = gehan, dist = "weibull")

Coefficients:
 (Intercept) treatcontrol
    3.515687    -1.267335

Scale= 0.7321944

Loglik(model)= -106.6   Loglik(intercept only)= -116.4
 Chisq= 19.65 on 1 degrees of freedom, p= 9.3e-06
n= 42

グラフはこう。
折れ線は survfit による生存関数のカプランマイヤー推定量(ノンパラメトリック)。なめらかな曲線がsurvreg による生存関数のパラメトリックモデル(ワイブル分布を仮定)。

グラフのコード:
> sr <-survreg(Surv(time,cens)~treat,dist="weibull", data=gehan) 
> sf <-survfit(Surv(time,cens)~treat, data=gehan)
> plot(sf)
> curve(1-pweibull(x,shape=1/sr$scale, scale=exp(coef(sr)[1])),add=TRUE)
> curve(1-pweibull(x,shape=1/sr$scale, scale=exp(coef(sr)[1]+coef(sr)[2])),add=TRUE)
survreg の scale = 1/(pweibull の shape)
survreg の Coefficients = log(pweibull の scale)
である。survreg はなぜか、ふつう shape と呼ぶパラメータが Scale という名前になっている。

weibreg
ehaパッケージをダウンロードしてインストールすると、weibreg という関数が使える。
> library(eha)
> weibreg(Surv(time,cens)~treat, data=gehan)
Call:
weibreg(formula = Surv(time, cens) ~ treat, data = gehan)

Covariate           Mean       Coef Exp(Coef)  se(Coef)    Wald p
treat
            6-MP    0.664     0         1           (reference)
         control    0.336     1.731     5.646     0.413     0.000

log(scale)                    3.516    33.639     0.252     0.000
log(shape)                    0.312     1.366     0.147     0.034

Events                    30
Total time at risk           541
Max. log. likelihood      -106.58
LR test statistic         19.6
Degrees of freedom        1
Overall p-value           9.29141e-06

グラフ。
weibreg による上と同じもの

グラフのコード:
> wr <- weibreg(Surv(time,cens)~treat, data=gehan)
> plot(sf)
> curve(1-pweibull(x,shape=exp(coef(wr)[3]), scale=exp(coef(wr)[2])),add=TRUE)
> curve(1-pweibull(x,shape=exp(coef(wr)[3]), scale=exp(coef(wr)[2]+-coef(wr)[1]/exp(coef(wr)[3]))),add=TRUE)
両者は、パラメータの置き方がちがう。

survreg の係数 = −(weibreg の係数)×(weibreg の shape)
survreg の Scale = 1/(weibreg の shape)

である。ぼくはかなり混乱した。

参考リンク
survreg vs. aftreg (eha)

R: weibreg
R: survreg

eha(イベントヒストリー分析)パッケージ中のオブジェクト一覧
survival(生存分析、ペナルティ付き尤度を含む)パッケージ中のオブジェクト一覧

2011年12月18日日曜日

Twitterウィルス(スパム)について

追記 2012/05/30

本家ヘルプセンターにちゃんとした記事が書かれていた↓



ぼくも前に一度ひっかかったのでメモっときます。


注意
ツイッターでこんな感じのダイレクトメッセージが来て、

誰かがブログでお前のこと悪く書いてるぞ、みたいな感じ

URLをクリックするとこんな感じにログイン画面に行った


としたら、それはスパムです。よくみるとURLが違う。twitter.comじゃなくて、twvitterとかになってる。

ここにIDとパスワードを入れると、今度は自分がスパムダイレクトメッセージをばらまく側に回ってしまいます。


対処
もしまちがってIDとパスワードを入れてしまったら、ツイッターの「設定」から、「アプリ連携」を選んで、許可したおぼえのないヤツとか、よく分からないヤツを全部消してけば(たぶん)大丈夫です。


おれのツイッターのアカウント名が写ってるが…

ついでにパスワードも変えとくといいかも。

スパムメッセージの内容は他にも「あのときの写真ができたよー、みてみて」みたいなのとか色々あります。

R: zero-inflated(ゼロ過剰)モデル

ほとんどが0でときどきピョコっと数字が入っているようなカウントデータ(なにかの数を数えたデータ)に対して回帰分析をするときは、zero-inflated(ゼロ過剰)モデルというやつを使うといいらしい。
zero-inflated モデルは、R で pscl パッケージをインストールしてzeloinflという関数を使うとできる。

例)
zeroinfl(y ~ x, dist = "poisson")


分布は "poisson"(ポアソン分布), "negbin"(負の二項分布), "geometric"(幾何分布)が使えるようである。

R: glmのオプション

R の glm 関数で


イテレーション(繰り返しの計算)回数を増やすには
glm() の引数 maxit に適当な数字をいれる。「アルゴリズムは収束しませんでした」みたいなメッセージが出ても、maxit を増やせばとりあえず収束はする(場合もある)。


イテレーションの収束判定(許容差)を変えるには
glm() の引数 epsilon に適当な数字をいれる。これを変えたいと思ったことはあんまりない。


例)
>  x <- c(18,17,15,20,10,20,25,13,12)
>  y <- c(0, 0, 0, 0, 1, 1, 0, 1, 0)
> glm(y ~ x, family=binomial(logit), maxit=1000, epsilon=0.0001) #この計算においてイテレーション回数とか変える必要は別にないが、いい例が思いつかなかった。
 
Call:  glm(formula = y ~ x, family = binomial(logit), maxit = 1000, 
    epsilon = 1e-04)

Coefficients:
(Intercept)            x  
     2.6828      -0.2111  

Degrees of Freedom: 8 Total (i.e. Null);  7 Residual
Null Deviance:     11.46 
Residual Deviance: 10.09  AIC: 14.09 


また trace = TRUE にすると各イテレーションごとに出力してくる。
例)

> glm(y ~ x, family=binomial(logit), maxit=1000, epsilon=0.0001, trace=TRUE)
Deviance = 10.09889 Iterations - 1
Deviance = 10.09096 Iterations - 2
Deviance = 10.09096 Iterations - 3

Call:  glm(formula = y ~ x, family = binomial(logit), maxit = 1000,
    epsilon = 1e-04, trace = TRUE)

Coefficients:
(Intercept)            x
     2.6828      -0.2111

Degrees of Freedom: 8 Total (i.e. Null);  7 Residual
Null Deviance:    11.46
Residual Deviance: 10.09 AIC: 14.09


参考URL: http://rss.acs.unt.edu/Rdoc/library/stats/html/glm.control.html

2011年12月3日土曜日

R:& と &&の違い

最近ようやくちょっと R に慣れてきた感がある

& と &&の違い
どちらも論理積なんだけど、

& だとベクトルの要素それぞれに対して
&& だとベクトルの最初の要素に対して

の演算になるみたいです。

例)
> c(TRUE,FALSE,FALSE) & TRUE
[1]  TRUE FALSE FALSE
> c(TRUE,FALSE,FALSE) && TRUE
[1] TRUE
> c(FALSE,TRUE,FALSE) & TRUE
[1] FALSE  TRUE FALSE
> c(FALSE,TRUE,FALSE) && TRUE
[1] FALSE

&& を使わなきゃいけない場面って少ない気がするので、とりあえず & をつかっておけばオッケーだと思います。| と || についても同様。

その他一般原則
  • R はとにかくベクトル‐行列が基本
  • なので for 文はなるべく使わないほうが計算が速くなる
  • このことはちょうど、数学でベクトル‐行列の計算をするとき、成分をいちいち書き下すよりもベクトル‐行列表記のままの方がかんたんだったりすることと対応している
  • for 文を使わないで代わりに何を使うのかというと apply とか merge とか subset とかである
  • ベクトルについての操作は大抵専用の速くてかんたんな関数があるので、そっちを使う
  • あと、一度書いたコードを後から見なおすのはけっこう大変なので最初からていねいに書くといい

2011年12月2日金曜日

養老孟司『バカの壁』への批判

動機
なんで今さらこんなエントリを書こうと思ったかというと、養老先生の「y=ax」を褒めてる人を最近になって読書メーターで見たからです。以下に養老先生の「y=ax」がなんでダメなのかを書きます。

y=ax
養老先生の「y=ax」というのはこんな感じです。

  • 人間の脳みそは式「y=ax」で説明できる
  • x は脳への入力(五感とか何らかの情報)、y は出力
  • a(の絶対値)は関心が大きい入力に対しては大きい値、関心が小さいことには小さい値をとる
まず、この「y=ax」っていう書き方が良くないです。入力に対していろんな値をとるなら a は定数じゃないことになります。 a を関数にして、y=a(x) とかって書くべきです。まあそれはいいや。この説明でなにが分かるかを見ていくことにします。

入力:関心大 → a:大
入力:関心小 → a:小

 でした。a が大きかったり小さかったりするとどうなるかというと、

入力:関心大 → a:大 → ax:大(= y:大)
入力:関心小 → a:小 → ax:小(= y:小)

そして、出力 y が大きいってことは、その入力に対しいての反応が大きいってことを表わしているんだそうです。なるほど。つまり、

入力:関心大 → a:大 → ax:大(= y:大) → 反応:大
入力:関心小 → a:小 → ax:小(= y:小) → 反応:小

人間の脳は、関心が大きいことには大きく反応するのか! でもこの「反応の大小」っていうのは何を表している量なんでしょう。どうも「反応の大きさ」はその入力に対する関心の大きさを表しているみたいです。

入力:関心大 → a:大 → ax:大(= y:大) → 反応:大 → 関心:大
入力:関心小 → a:小 → ax:小(= y:小) → 反応:小 → 関心:小

途中の y とか a とかを省略して結果をみてみます。

入力:関心大 → 関心:大
入力:関心小 → 関心:小

これですっきりしました。実に論理的。もっとも、この説明は一体なにを説明できているのかというと……たぶん何も説明できてないと思います。

(余談:数学をやってるとたまにこういうことがある。「よし、計算できそうだ。解けた! 0 = 0 ……ってこりゃ当たり前だよ。何やってんだおれ……」みたいな。)

この「y=ax」はもちろん褒めてる人だけじゃなくて貶してる人もいっぱいいたのですが、「あまりに単純化がひどい」みたいな趣旨のものが多く、私が言ってるようなことを書いてる人は少ないように見えました。「単純化がひどい」っていう批判も正しいのでしょうが、それ以前にこのモデル自身がまったく無意味だと思います。

ほかにも、
  • フロイトの無意識について述べた部分とかもおかしいです。「最近、おれの授業で寝てる学生が多い」 → 「学生は睡眠時間のことをちゃんと考えてない」 → 「都会っ子は無意識を意識してないのだ! だからフロイトは無意識を再発見しなければならかったのだ!」みたいな。フロイトの言う無意識って寝てる時間のことじゃないんですけど……。フロイトの「無意識」は意識がない状態を指す「無意識」とはまったく別もの。内田樹もなかよく対談なんかしてないで、つっこんでやればいいのに。まちがってることにまちがってると言うのも友情だと思うな、ぼくは。
  • ピカソについて述べた部分もおかしいです。「ピカソの(キュビズム時代の)絵はめちゃくちゃなようでいて、顔のパーツ部分など細部はリアルにかいてある」 → 「ピカソは見た景色を脳の中で独自に再構成できたのだ!」みたいなの。顔のパーツ部分、リアルか? まあ、ものの感じ方は人それぞれですからね……みたいな解釈もできるが、たぶん養老先生、ピカソの絵見てないだろ。それに、パーツ部分がリアルだったしても、ピカソの脳自体が特殊だと考えるのはかなり飛躍してると思います。
  • あとなんだっけ、共同体の再構築みたいな話にも不満があったような気がするが、忘れた。
  • それと、なんていうか、みんな「人それぞれいろんな意見がある」とか、「物ごとには良い面と悪い面がある」みたいな話好きだよなー。いや、それ自体はいいんだけど。それは議論の出発点であって結論じゃないよー。


おわりに
ぼくが「バカの壁」を読んだのはかなり前、中学生(いや、高校生だったっけ? 忘れた)のころだったので、このエントリは記憶をたよりに書いてます。なので、引用などは不正確です。でも大筋では間違っていないはずです。