« 2008年3月 | トップページ | 2008年5月 »

2008.04.25

[大学]地名が最後につく学校名

 大学などの名前で地名が最後に付くのは珍しい(日本語での場合である)。よく知られている例は、もちろん首都大学東京だが、他にはないかと思っていたら、日本航空大学校山梨というのを最近見かけた(ただ、上記のページから日本航空大学校山梨のウェブサイトに行くとForbiddenとなる)。日本航空大学校や日本航空大学校山梨は、航空大学校(国土交通省の所管で、いまは独立行政法人)とはまったく別物だそうだ。日本のない方の、航空大学校は、短大卒が受験資格ということでもけっこう有名である。

 日本が名前の上につくかどうかは重要なちがいということであろうか。

| | コメント (0) | トラックバック (0)

2008.04.24

[統計]一致性

 統計で一致性(consistency)といえば、推定で、サンプル数が大きくなると推定値が真の値に近づいていくこと、別の言い方をすればサンプル数無限大では推定値が真の値と一致することである。一致性は、推定量が持つべき性質と考えられていることが多い(たとえば『・・・推定量が持つべき性質として一般に要求される;』と「正規分布」(柴田義貞著、東大出版会)にある)。
 少し話を広げてサンプル数つまりデータの量(一致しない場合もあるが)が増えていくと推論の結果が正しいものに近づき無限大になれば必ず正しいものに到達できることを、一致性と考えてみる(直感的でもあり、よくある使い方で、別に珍しいものではありません)。[以下は、尤度比検定やAICを自分でデータ解析に使って、あれこれといろいろ頭をひねって考えた人はきっと通過していることだと思います]

 統計的検定には一致性があるかと考えてみる。帰無仮説にあたるモデルが正しいときに、有意水準をαとすると、100α%では対立仮説が選ばれる。100α%は正しくない方が選ばれることになるのだから、上記の意味では、統計的検定には一致性がないことになる。

 モデル選択ではどうだろうか。たとえば、モデル選択の代表格と考えられることが多いAICで、サンプル数がすごく大きくなると正しいモデルが選ばれるようになるのだろうか、と考えてみる。単純な場合を考えてみる。2つのモデルに包含関係が成り立っていて(つまり尤度比検定も可能な状況ということになる)、単純なモデルが正しく、サンプル数はとても大きいというケースを考えてみる。尤度比検定の議論から、両モデルの対数尤度の差の2倍(対数尤度比統計量)は、自由度がパラメーター数の差であるようなカイ2乗分布する。AICの差と対数尤度の差の2倍の関係はシンプルなので、整頓すると、このケースで正しいモデルが選ばれるのは、自由度=(パラメーター数の差)であるカイ2乗分布で、パラメーター数の差の2倍より小さな値の場合である。カイ2乗分布の分布関数から、このケースの正しいモデルが選ばれる確率は、たとえば、パラメーター数の差が1のとき0.84つまり84%くらいになる。サンプル数がとても大きくても正しいモデルが選ばれない確率は1割を超えており、やはり上記の意味では、一致性はないことになる。

注記1:なお、わかりきったことかもしれませんが、上に書いた統計的な推論(AICも推論に仮に含めています)が変なものだという意味で書いてはおりません(その人のその時の使用目的に合わないものであるということは当然ありえますが)。

注記2:カイ2乗分布の分布関数は、表を引いてもいいのだが、Rで以下のようにすぐ計算できます。
> pchisq(2,df=1)
[1] 0.8427008
なお、パラメーター数の差が大きくなれば、確率は大きくなりますが、ここでの議論には影響がありません。

| | コメント (0) | トラックバック (0)

2008.04.21

[本]『ゴールデン・トライアングル秘史』

  競泳の日本選手権のニュースをちらちら見ながら、仕事の合間に、タイ北部の台湾飛び地とも形容できそうな地域について書いてあるということで、『ゴールデン・トライアングル秘史』(NHK出版)を読んだーだいぶ前に読むはずだったのだが、そのままになっていた。国共内戦に敗れた国民党軍は台湾に逃れたのだが、雲南にいた国民党軍(1万2千人とWikipediaにはある)はビルマへ、そしてタイ北部に逃れたとのことだ。

| | コメント (0) | トラックバック (0)

2008.04.20

[その他]正7角形

 角の三等分は、定規とコンパスだけではできない作図問題としてあまりに有名で、三大問題の1つでもある。もっとも少し使える器具の制限を緩めれば可能で、とコンパスと定規ではなく、コンパスといわゆる物差しならできることがあちこちに説明されている。この作図法は、アルキメデスによるそうである。
 正多角形のうち、3,4,5,6角形が定規とコンパスだけで作図できることはよく知られている(はるか昔から)。8とか12とか15とかもできることがすぐわかってしまう。できない最小のものは正7角形で(たしかガウスが証明した)、これも少し制約を緩めてやると作図できることがアルキメデス以来わかっているそうである。

| | コメント (4) | トラックバック (0)

2008.04.19

[その他]ハトロン紙とパトローネ

 フィルムカメラが普通にはほとんど使われなくなって、パトローネはほぼ死語になっている(パトローネは、35mmフィルムなどが入っているケースです。フィルムのすぐ外側の、フィルムを巻き取る棒のような突起がついているものであり、その外側のプラスティックケースのことではありません)。このパトローネは、ハトロン紙と同じ語源(薬莢のようなもの)だそうだ。
 35mmのフィルム代節約のため、長巻きのフィルムを短く切って暗室でパトローネに詰めるのが、そういえば普通のことだったと思い出す。100フィートだったかが缶に入っていた。フィルムの面にさわらないよう暗室で手作業するのが、私の場合は普通だった。後で詰める機械(ローダー)があるのを知った。

| | コメント (0) | トラックバック (0)

2008.04.18

[本]『日本のソフトウェア産業がいつまでもダメな理由』

 なぜか(生態学会大会前ほどではないが)いろいろな用(範疇としてはいわゆる雑用)がやってきて、本屋にいく時間もあまりないのだが、久手堅憲之(著)『日本のソフトウェア産業がいつまでもダメな理由 』(技術評論社)を立ったまま読んだ。日本のソフトウェア産業は業界あげての英語アレルギーなのだそうだ。そのため(と思ったのは私である)、カタカナはとても好まれていて、ということで、例として、コラボする、リスケする、などがあげられている(p.110)。
 医療の現状について書かれた本を読むときの身につまされ方とは少しちがうが、ソフトウェア産業について書かれたこの本に出てくる人の模様は私が暮らしている業界に似ている(上記の英語の点ではだいぶちがいますが)。

| | コメント (0) | トラックバック (0)

2008.04.16

[統計][R]Rによる混合効果ロジスティック回帰(余談2)

 他のことに(私が)使われているあいだに、East_Scrofaさんや足軽日記さんのところで、だいぶ進んでいる。

 また、余談であるが、さきに書いた2つのデータの与え方を、glm()だけでなくglmmML()にも試してみる。結果はここまでこの話題で書いてきたことから容易に想像されるので、「glm()のAICとglmmML()のAICはこんなときだったら比べられるんだろ」とすぐ分かってしまう方は読み飛ばしてください。使ったのはこんなデータである。
まず、1行が1cluster
> mydata1g
alive dead x1 total ID
1 4 6 1 10 1
2 7 3 2 10 2
3 8 2 4 10 3
4 3 4 0 7 4

次に内容的には同じデータの1行1個体版である。
> mydata1
alivedead x1 ID
1 1 1 1
2 1 1 1
3 1 1 1
4 1 1 1
5 0 1 1
6 0 1 1
7 0 1 1
8 0 1 1
9 0 1 1
10 0 1 1
11 1 2 2
(中略)
30 0 4 3
31 1 0 4
32 1 0 4
33 1 0 4
34 0 0 4
35 0 0 4
36 0 0 4
37 0 0 4
長いので少し略した。
> res.glm.i1<-glm(alivedead~x1,family=binomial,data=mydata1)
> res.glm.g1<-glm(cbind(alive,dead)~x1,family=binomial,data=mydata1g)
とまず、glm()で、次に
> res.glmmML.g1<-glmmML(cbind(alive,dead)~x1,family=binomial,cluster=ID,data=mydata1g)
> res.glmmML.i1<-glmmML(alivedead~x1,family=binomial,cluster=ID,data=mydata1)
とglmmML()でやってみる。

 以下の2つはglm()の結果で、先に書いたとおりで当たり前だが、点推定値や検定結果などまったく同じである。
> summary(res.glm.i1)

Call:
glm(formula = alivedead ~ x1, family = binomial, data = mydata1)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.8363 -1.1767 0.6401 0.9766 1.3939

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4959 0.5638 -0.88 0.3791
x1 0.4943 0.2657 1.86 0.0629 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 49.961 on 36 degrees of freedom
Residual deviance: 46.009 on 35 degrees of freedom
AIC: 50.009

Number of Fisher Scoring iterations: 4

> summary(res.glm.g1)

Call:
glm(formula = cbind(alive, dead) ~ x1, family = binomial, data = mydata1g)

Deviance Residuals:
1 2 3 4
-0.6320 0.5252 -0.1188 0.2711

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4959 0.5638 -0.88 0.3791
x1 0.4943 0.2657 1.86 0.0629 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 4.71427 on 3 degrees of freedom
Residual deviance: 0.76285 on 2 degrees of freedom
AIC: 15.016

Number of Fisher Scoring iterations: 4

 以下の2つはglmmML()の結果で、こちらも、点推定値や検定結果などまったく同じである。
> summary(res.glmmML.i1)

Call: glmmML(formula = alivedead ~ x1, family = binomial, data = mydata1, cluster = ID)


coef se(coef) z Pr(>|z|)
(Intercept) -0.4959 0.5638 -0.8795 0.3790
x1 0.4943 0.2657 1.8601 0.0629

Standard deviation in mixing distribution: 1.041e-08
Std. Error: 0.4001

Residual deviance: 46.01 on 34 degrees of freedom AIC: 52.01
> summary(res.glmmML.g1)

Call: glmmML(formula = cbind(alive, dead) ~ x1, family = binomial, data = mydata1g, cluster = ID)


coef se(coef) z Pr(>|z|)
(Intercept) -0.4959 0.5638 -0.8795 0.3790
x1 0.4943 0.2657 1.8601 0.0629

Standard deviation in mixing distribution: 1.611e-09
Std. Error: 0.4001

Residual deviance: 0.7628 on 1 degrees of freedom AIC: 6.763

 さて、この例では、ほとんどランダム効果の項の効果はないのだが(したがってあまりわかりやすい例ではなくてすみません)、1行1個体という与え方のときのglm()とglmmML()が返すAICを見ると、これが同じ基準で計算されているらしいことがわかる。この例ではほとんどランダム効果項の効果がないので、1パラメーター多い分だけ、つまりAICなら2ちがってくるはずだが、実際にもそうなっている。
 glm()とglmmML()が返すAICは、1行1clusterというデータの与え方のときは基準値が違うのでそもそもそのままでは比較可能でなく、一方、1行1個体という与え方のときは同じ基準値で計算されているので比較しうるものになっている、ということのようだ。
 1行1clusterというデータの与え方のときは、glmmML()にclusterがないときのdevianceを計算させ(いわばglm()を代行させ)れば、同じ基準値によるAICはすぐ求められるわけである。

| | コメント (0) | トラックバック (0)

[香椎・千早]香椎浜団地入口交差点ー洋服の青山

 国道3号線と3号線バイパスを直進した道が交差する、香椎浜団地入口交差点だが、南西側角の、もとコンビニがあった場所には、洋服の青山ができると看板が出ている。建物は、コンビニくらいのサイズの平屋だけだった。これだけなのかと思っていたら、3号線沿いで別の建物らしきものを作っているーあれも青山なのだろうか。

| | コメント (0) | トラックバック (0)

[その他]Safariの縦スクロールバー

 いろいろなOSのパソコンを使っているのだが、少し前から、MacOS X (10.4)で動いているマシンの、Safari(ソフトウェアのアップデートでは一番新しいと言ってくる)で、縦スクロールバーが出ないことが多い。まだ、下にかなり表示されていないものがあっても、右側に出るはずの縦スクロールバーはなく、矢印キーも効かない。新しいウィンドウかタブを出せば、そっちではたいてい縦スクロールバーがあるので、深刻でもないが、面倒ではある。

| | コメント (0) | トラックバック (0)

2008.04.15

[その他]NBAプレーオフ進出チーム確定

 ゴールデンステイトはフェニックスに負けて(全部勝たなければ8位はなく、全部勝っても同勝率ではだめだったのだが)、あっさり、デンバーのプレーオフ進出が決まった。東西とも全チーム決まり、あとはシード順(とくに西)だけである。

| | コメント (0) | トラックバック (0)

2008.04.14

[その他]生態学会大会ー忘れていたこと・やり残し

 福岡での生態学会大会から約1月たつ。(大会実行委員会の仕事はいまも少しありますが、以下は実行委員会の仕事以外のことです)生態学会大会関係で、いくつか、すでにやった気になっていたり、いつやろうかと考えていることがある。すでにやった気になっていたことの1つが、自由集会で話した内容(PDF)のウェブサイト掲載である。いくつか指摘された内容を注記したり、重複スライドなどを整理したりして、できていたので、ウェブサイトに載せたつもりになっていた。先週載っていないのに気づきあわてて載せた。
 自由集会での話では、以前せっかく久保さんから教えてもらったのに、それよりもわかりにくいプログラムになっていた点があるのだが、ウェブサイトのPDFでは注記してある。先週、メモ類を見返してみると、昨年12月のあるところまではちゃんとしているのに、とある時点から急にわかりにくい(しかも誤りが起こりやすい)書き方になっている。理由らしきメモ書きもなく、”なぜだ”とやや自己嫌悪気味になっていた。
 すでにやった気になっていることやいつやろうかと考えていること以外に、生態学会大会でいただいた(ありがたい)提案については進めていますー遅遅としているかもしれませんが。

| | コメント (0) | トラックバック (0)

[その他]替え歌

 少し前(といっても生態学会大会よりも前なので、印象としては3年くらい前)に、ある催しのときに替え歌を披露しようかと考えていたことがあり(結局そんな時間はなかった)、いくつか考えてみた。
 新しく研究室に来た、卒論生ないし大学院生が教員に対して言うという想定の『関白宣言』(さだまさし)の替え歌は、汎用性があるのではないかと思う。以下は適当に考えたもののうちにとくにあたりさわりのない例である:

お前の指導受ける前に言っておきたいことがある
かなりきびしい話もするが俺の本音を聞いておけ
(中略)
指導は思いつきでするな
いつもよく調べておけ
できる範囲でかまわないから

忘れてくれるな
経験値でいえば
お前は俺の何万倍もあるってことを
(中略)
俺は逃亡はしない
たぶんしないと思う
しないんじゃないかな
まちょっと覚悟はしておけ

| | コメント (1) | トラックバック (0)

2008.04.13

[その他]NBAプレーオフあと1つずつ

 西はやはりダラスは逃げ切り、デンバーとゴールデンステイトの激しい争いになってきた。残り2試合で同率、勝率6割で8位を争っているわけである。フェニックス(サンズ)が1970年代に、49勝33敗でプレーオフにいけなかったことがあるそうだ。

| | コメント (0) | トラックバック (0)

[統計][R]Rによる混合効果ロジスティック回帰(余談1)

 ここ何回か、Rによる混合効果ロジスティック回帰について書いてきた(といっても出てくるのは、glm()とglmmML()だが)。この2つについても、AICや対数尤度、あるいはdevianceの計算という、いままで書いてきたことと似た点についていくつか補足しておく。
 ロジスティック回帰では、ソフトウェアへのデータの与え方にいくつかのやり方がある。ここのところ使っていた、各容器の中にいる個体たちの生存/死亡なら、容器ごとに、生存個体数と死亡個体数を与えるやり方もあれば、1個体ごとに生存/死亡(そして容器の番号なども)を与えるやり方もある。(その1)から(その4)で使っていたのは前者である。

 Rのglm()という関数はかなりよく使われているので、当然の知識かもしれないが、上記の2つのやり方で、内容的には同じデータを与えたとき、glm()の返す結果は表面上は違っているところがある。とても簡単な例でやってみる。
まず使うデータだが、以下のmydata1gである。これは”前者”にあたる。
> mydata1g
alive dead x1
1 4 6 1
2 7 3 2
3 8 2 4

もう1つはmydata1gを個体ごとに直した以下のmydata1である。glm()で二項分布という仮定で行うなら(つまりoverdispersionもanderdispersionもないなら)、内容的には同じである。
> mydata1
alivedead x1
1 1 1
2 1 1
3 1 1
4 1 1
5 0 1
6 0 1
7 0 1
8 0 1
(長いので途中は省略)
27 1 4
28 1 4
29 0 4
30 0 4
以下のようにそれぞれロジスティック回帰する
> res1g.glm<-glm(cbind(alive,dead)~x1,family=binomial,data=mydata1g)
> res1.glm<-glm(alivedead~x1,family=binomial,data=mydata1)

> summary(res1g.glm)

Call:
glm(formula = cbind(alive, dead) ~ x1, family = binomial, data = mydata1g)

Deviance Residuals:
1 2 3
-0.3984 0.6178 -0.2580

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.7318 0.8246 -0.887 0.3748
x1 0.5816 0.3518 1.653 0.0983 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 3.74390 on 2 degrees of freedom
Residual deviance: 0.60697 on 1 degrees of freedom
AIC: 12.41

Number of Fisher Scoring iterations: 4

> summary(res1.glm)

Call:
glm(formula = alivedead ~ x1, family = binomial, data = mydata1)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.886 -1.114 0.608 1.001 1.242

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.7318 0.8245 -0.888 0.3748
x1 0.5816 0.3518 1.653 0.0983 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 36.293 on 28 degrees of freedom
AIC: 40.293

Number of Fisher Scoring iterations: 3

推定値や検定結果などは当然のことだが同じなのだが、AICやdevianceなどが異なっている。devianceもNullとResidualの差はどちらでも同じである。
 まず、AICのちがいが計算されている対数尤度の違いのみによることは、
> logLik(res1g.glm)
'log Lik.' -4.205005 (df=2)
> logLik(res1.glm)
'log Lik.' -18.14627 (df=2)
としてみると、対数尤度の差がAICのちがいのちょうど半分であることから確かめられる。

 それぞれの対数尤度は以下のように計算されている。”容器”ごとの場合の、ロジスティック回帰による生存確率の予測値は
> res1g.glm$fitted
1 2 3
0.4625256 0.6062115 0.8312628
である。そこで、対数尤度のうち、生存数×log(生存確率)+死亡数×log(1-生存確率)の部分は
> sum(mydata1g$alive*log(res1g.glm$fitted)+mydata1g$dead*log(1-res1g.glm$fitted))
[1] -18.14627
となる。これに二項係数の部分、
> sum(lchoose(mydata1g$alive+mydata1g$dead,mydata1g$alive))
[1] 13.94126
を加えたのが、 -4.205である。個体ごとというときには、二項係数の部分が入っていないわけである。

 ”容器”ごとの場合の飽和モデルは、各”容器”にそれぞれ尤度がもっとも大きくなる生存確率を割り当てたもので、alive/(alive +dead)である。p0がその尤度がもっとも大きくなる生存確率で、以下のようになる。

alive dead x1 total p0
1 4 6 1 10 0.4
2 7 3 2 10 0.7
3 8 2 4 10 0.8

生存数×log(生存確率)+死亡数×log(1-生存確率)の部分は-17.84278となり、二項係数の部分(上に書いたのと同じく 13.94126)を加えると、-3.9015となる。この飽和モデルの対数尤度とロジスティック回帰の対数尤度(上記の-4.205)のちがいの2倍が、Residual devianceになっている。

 ”個体”ごとのときには、それぞれの個体の結果を完全に予測してしまうのが飽和モデルになるから、尤度は1で対数尤度は0になる。そこで、-18.14627という対数尤度の2倍がResidual devianceになっている。これは、glmmML()のときの計算の仕方そっくりである(glmmML()のデータの与え方が拡張されてきた経過を考えると、このように計算しているのがわかる気がする)。


 

| | コメント (0) | トラックバック (0)

2008.04.12

[統計][R]Rによる混合効果ロジスティック回帰-数値的な積分

 (その4)では数値的な積分をして、確認した。Rではintegrate()という数値積分する関数があり、integrate(関数名,下限,上限)で積分できる((その4)での計算ではintegrateと台形公式のプログラムを書いたものと両方でやってみた)。ただ数値積分するだけなら、Matlabやそれに近いソフトウェア(ScilabとかOctaveとか)の方が便利かもしれない。

| | コメント (0) | トラックバック (0)

2008.04.11

[その他]桜の葉も・・・

 『うろんなことをまた桜花』と連想しているうちに、昨日は大風が吹いて、桜の葉が目立つようになってきた。年度がかわって、毎年のこととはいえ、ネットワーク関係のいろいろな問題への対応に追われている。新しいマシンをつないだり、新しく変更があったりと、つながらなかったりうまく通信できなかったりする頻度が高い。
 せっかく暖かくなったのに山の中に出かける時間ができない。考えてみると、今週は、漫画としては「とめはねっ」の3巻しか読んでいないような気がする。

| | コメント (0) | トラックバック (0)

[統計][R]Rによる混合効果ロジスティック回帰(その4)

 (その3)では、「glmmML()のAICの計算に使われている対数尤度は、residual devianceを計算する基準となっているモデルの対数尤度を0として出したものらしい」という内容のことを書いた。では、この「residual devianceを計算する基準となっているモデル」とはどんなものということになるが、普通に考えると一番ありそうなのは、飽和モデル、つまり、ここでは、各”容器”に異なる生存確率というパラメーターをそれぞれ割り当てた場合である。

 (その2)の数値例では、こんな風になる

ID alive dead x1 パラメーター
1 2 8 1 0.2
2 5 5 2 0.5
3 6 4 3 0.6
4 9 1 4 0.9
5 3 7 3 0.3

このときの対数尤度は、二項係数の部分を入れないと、-28.025(ちなみにいれると-6.25)となる。

 glmmML()が返してきた推定結果、つまり回帰係数、切片、ランダム効果の項の標準偏差を使って、実際に数値的に積分してみると、対数尤度は、二項係数の部分を入れないで、-30.585であった。上の飽和モデルとの差は2.56で、2倍すると5.12と、glmmML()でのResidual devianceの値と同じである。どうやら、このような飽和モデルの対数尤度を0として、glmmML()ではAICを計算しているらしい。
 AICは相対的な大小が問題になり、1つだけであれこれ言うことはまずないから、基準をここにとるのに何の問題もないが、別のところに基準をとっている結果を返してくる関数などと比べるときには注意が必要である。基準を合わせてやらないといけないわけである。
 もっともランダム効果の項がない(ただの)ロジスティック回帰と比べるなら、glm()でランダム効果の項がない場合を計算して(これは普通はそうするかもしれないが)対数尤度計算の基準を合わせてもいいが、(その2)で述べたようにglmmML()に同じ基準を使ったときのdevianceを変えさせたほうが簡単だろう。cluster.null.devianceを使えばいいわけである。この例では、devianceで(対数尤度だとその半分になる)0.014ほどしか変わらない。パラメーターが1つ(ランダム効果の項の標準偏差)ちがうから、ランダム効果の項がない(ただの)ロジスティック回帰のモデルのAICは、glmmML()で言えば、9.135となる。

 ここまでで(その1)に書いた疑問「glmmML()のAICの値は小さすぎないか(対数尤度の値は大きすぎないか)」
は一応解け、glm()などの返す結果と比較可能にはなったと思う。[この疑問とそれ以降調べたことは、何か変だと漠然と感じていたところをEast_Scrofaさんの日記に触発されたものです]さて、まだ余談はある。

| | コメント (0) | トラックバック (0)

2008.04.09

[香椎・千早]桜

 先週末(5日)、少し足を伸ばして(といってわずかな距離なのだが)、地元香椎浜界隈の桜を見て来た。この時点ではまだ咲ききっていない樹もあり、東京よりはだいぶ遅いようだった。香椎浜2丁目には香椎浜小学校と城香中学のところから北に伸びる桜並木があり、1丁目の東西に伸びる桜並木がある(西の端が都市高速と南公園のところになる)。

| | コメント (0) | トラックバック (0)

[統計][R]Rによる混合効果ロジスティック回帰(その3)

 (その1)と(その2)では、計2つのデータセットで調べてみると、AICとresidual devianceの差が同じでしかも整数らしかった。2つではとてもではないが心許ないので、乱数を使っていろいろ作ってやってみる。

n1<-100
a1<-numeric(n1)
rd1<-numeric(n1)
x1<-seq(6)
ID1<-seq(6)
for (i1 in 1:n1){
y1<-floor(runif(6,min=1,max=20))
y2<-floor(runif(6,min=1,max=20))
tdata1<-data.frame(x1,ID1,y1,y2)
res1<-glmmML(cbind(y1,y2)~1+x1,family=binomial,cluster=ID1,data=tdata1)
a1[i1]<-res1$aic
rd1[i1]<-res1$deviance
}
で100通りのデータでやってみると、結果は、

> a1-rd1
[1] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[42] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
[83] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
となって、いつでも6だった。どうやら、residual devianceを計算する基準となっているモデルの対数尤度を0としてAICを求めているらしい。

| | コメント (0) | トラックバック (0)

2008.04.08

[統計][R]Rによる混合効果ロジスティック回帰(その2)

 先の計算を、もう少し”ばらついた”データでもやってみる。
> mydata12
ID alive dead x1
1 1 2 8 1
2 2 5 5 2
3 3 6 4 3
4 4 9 1 4
5 5 3 7 3

というデータを使ってみる。

> glmm.res12x<-glmmML(cbind(alive,dead)~x1,family=binomial,cluster=ID,data=mydata12)
> glmm.res12x

Call: glmmML(formula = cbind(alive, dead) ~ x1, family = binomial, data = mydata12, cluster = ID)


coef se(coef) z Pr(>|z|)
(Intercept) -2.281 0.9803 -2.327 0.0200
x1 0.873 0.3539 2.467 0.0136

Standard deviation in mixing distribution: 0.1856
Std. Error: 0.8312

Residual deviance: 5.121 on 2 degrees of freedom AIC: 11.12

AICとResidual devianceの差は先の例と同じである。

 ところで、このデータをglmでやってみると以下のようになる。
> glm.res12x<-glm(cbind(alive,dead)~x1,family=binomial,data=mydata12)
> glm.res12x

Call: glm(formula = cbind(alive, dead) ~ x1, family = binomial, data = mydata12)

Coefficients:
(Intercept) x1
-2.2575 0.8628

Degrees of Freedom: 4 Total (i.e. Null); 3 Residual
Null Deviance: 13.26
Residual Deviance: 5.135 AIC: 21.64

なお、glmmMLの結果からも、glmを適用したときと同じ(clusterがなかったらという場合)Residual Devianceは、以下のように引っ張りだせる。
> glmm.res12x$cluster.null.deviance
[1] 5.135012

glmmMLとglmの結果では、Residual devianceはそんなに変わらないのに、AICはずっとglmmMLが小さい。

| | コメント (0) | トラックバック (0)

[その他]NBA

 ボストンの最高勝率は決まりのようで、東のプレーオフ最後はアトランタになりそうだが、西はデンバーがいたいところで星を落としたので、最後は2チームの激しい争いになりそうだ(ダラスは逃げ切りそう)。西は上の方のシード順もまだ入れ替わる可能性が大きい。

| | コメント (0) | トラックバック (0)

[統計][R]Rによる混合効果ロジスティック回帰(その1)

 Rで混合効果ロジスティック回帰というのは、いろいろなところで取り上げられていて、今さら感もある話題なのだが、細かいところではときどき「これは何」と思うこともある。以下、シンプルなものを使って考えてみる。
 まず、サンプルデータだが、以下のようにmydata11という名前をつけたデータフレームに入れておく。容器の中に9個体か10個体をそれぞれ入れておいて、x1で表される条件に置いておいたら、何個体が生存(alive)で何個体が死亡(dead)なのかというデータである。よくあるシンプルな状況だと思うが、データ自体は計算例としての仮想的なものである(生存か死亡のどちらかだけで半死半生とかはないことにする)。x1は連続量である。

> mydata11
ID alive dead x1
1 1 2 8 1
2 2 4 6 2
3 3 6 4 3
4 4 8 2 4
5 5 4 5 3

(CRANのパッケージだけを考えても)混合効果ロジスティック回帰はたとえば、glmmMLパッケージの関数glmmMLでできる。少し前にデータの取り扱いに変更があった(久保さんのサイトなど参照)。以下のように分析してみる、
> glmm.res11<-glmmML(cbind(alive,dead)~1+x1,family=binomial,cluster=ID,data=mydata11)
> glmm.res11

Call: glmmML(formula = cbind(alive, dead) ~ 1 + x1, family = binomial, data = mydata11, cluster = ID)


coef se(coef) z Pr(>|z|)
(Intercept) -2.2550 0.9271 -2.432 0.0150
x1 0.8466 0.3296 2.569 0.0102

Standard deviation in mixing distribution: 2.947e-08
Std. Error: 0.3394

Residual deviance: 0.7813 on 2 degrees of freedom AIC: 6.781

こんな風に書いてもできる、
> glmm.res11x<-glmmML(cbind(alive,dead)~x1,family=binomial,cluster=ID,data=mydata11)
> glmm.res11x

Call: glmmML(formula = cbind(alive, dead) ~ x1, family = binomial, data = mydata11, cluster = ID)


coef se(coef) z Pr(>|z|)
(Intercept) -2.2550 0.9271 -2.432 0.0150
x1 0.8466 0.3296 2.569 0.0102

Standard deviation in mixing distribution: 2.947e-08
Std. Error: 0.3394

Residual deviance: 0.7813 on 2 degrees of freedom AIC: 6.781

もちろん、結果はまったく同じである。さて、ぱっと見て、あれれと思うのは(私が思ったということだが)AICの値である。AICは、対数尤度の-2倍にパラメーター数を加えたものだから、6.781ということは対数尤度は-0.39くらいとなる。各容器ごとの生存確率をパラメーターだと考えても、普通に考えると一番対数尤度が大きいのは、それぞれの容器の生存率が生存確率に等しいときのはずである。そのときの対数尤度として-0.39みたいなこんな大きな値になるのだろうかと感じるわけだが、計算してみるとそのときの対数尤度は-6.51くらいである(二項係数の部分も入れた値、入れなければもっと小さくなる)。-0.39はそれよりもだいぶ大きいので、そんな値をとるのは変だ、というのが、問題の1つの発端である。

| | コメント (0) | トラックバック (0)

2008.04.07

[その他]NBAプレーオフ間近

 西は、まれに見る激戦で、ごく最近までプレーオフ進出チームが決まらなかった。ようやく3チーム決まったが、あと5試合程度を残すだけになったいまでも、7,8,9位はどうなるかわからない。ダラスやデンバーはどうなるのだろうか。まして、シード順は最後まで決まりそうにない。

| | コメント (0) | トラックバック (0)

2008.04.06

[香椎・千早]香椎浜団地入口交差点

 国道3号線と3号線バイパスを直進した道(3号線バイパス自体は右にまがり水谷の方に向かう)が交差する、香椎浜団地入口交差点のまわりは、バイク屋さん(絵がうつる電光掲示板があった)が名島よりに移転し、コンビニがなくなって(洋服の青山ができると書いてある)、だいぶ空いている印象である。
 久しぶりに3号線を北九州方面から走ってきてみると、東部市場前と産大前のあいだの東側(北九州方面から来ると左側)で工事をしている。だいたいあのあたりに、3号線バイパスが出てくることのだろうと思うが、3号線バイパスとの分岐点の工事なのだろうか。

| | コメント (0) | トラックバック (0)

« 2008年3月 | トップページ | 2008年5月 »