不運なRejectを避けるために

主著の論文が10本を越えてから、やっとIF > 4以上の中堅雑誌(Ecologyの分野では)に載せることができるようになってきた。それまでもレビューに回ることはあったが、一度とてアクセプトにこぎつけることはできなかった。さて、その勝敗の差を分けたのはなんなのであろうか、と自分なりに整理してみる。ただし、以下は十分良い結果が得られている前提で、いかにその結果を美味しく料理するか、という段階であるのであしからず。

一般性の担保

これは言わずもがなかもしれないが、IF > 4になると、特定の生態系、あるいは特定の分類群にしか当てはまらないことが明らかな現象は、Editor Rejectを食らう運命にある。ただ、普通は一種、よくても数種しか扱えないと思うので、いかに現象論として一般化できるのかをイントロとディスカッションで大きく、かつ説得力のある形で打ち出せるかが勝負どころ。大きく出るだけであれば簡単だが、やはり説得力を持たせるためにはかなり綿密な既存文献のレビューが必要になる。

新奇性の見せ方

新奇性があるにも関わらず、ぐだぐだと長い文章を書き連ね、新奇性がイントロの後のほうに出てくるのはよくない気がする。重要な知見の欠落を早い段階で指摘し(できれば2段落目、遅くとも3段落目)、そのギャップをいかにインパクトのある形で見せれるかが違いを生み出す気がする。

手法の説明

ここは見落としていたポイントだった。適切な手法を使っていようと、それが読者に伝わらなければ、それはトンチンカンな理由からリジェクトを食らう原因となる。私はベイズモデルをよく利用するのだが、その利点をあまり説明しなかったがために、「古典的なモデルを使ったほうが良い」という理由からリジェクトされた経験が何度かある。最近は、GLMMなどと比較して、こうした大きな利点がある、というのを明記するようにしており、この理由によるリジェクトはグッと減った。何度も「おれの方が正しいのに!」と思ったが、思い返してみれば「自分の手法の正しさをわかりやすく記述していない」ことが原因にあったのだ。自業自得であった。

結果の書き方

これも盲点だったのだが、単に「~の有意な効果が得られた」ではいけない。「コントロールと比べ、3~4倍もの~が生じることがわかった」のような、定量的なインパクトを示す記述が必要。いわゆる普通のジャーナルとトップジャーナルのResultの書き方を見比べると、この点は雲泥の差があることがわかる。印象が全然違うのだ、きっと。

弱点のカバーの仕方

どんな研究にも弱点はある。だが、弱点をただ認めるだけでいいのは下のレベルのジャーナルだけのような気がする。いかに弱点を「弱点らしからぬフリ」で守り抜くかが実は大事なテクニックであるように感じる。いかにさりげなく、だけどしっかりと研究の穴をガードするのか、これはなかなか難しく、数をこなさないとできないのだろう(まだ検討中)。

Supporting Informationの有効活用

想定されるネガティブなコメントで対処できるものであれば、面倒くさがらずにSIのなかで「こういう場合も想定されるが、それでも同様の結果が得られる」といったような記述を加えておく。上のJournalになるほど、ちょっとしたネガティブなコメントが入るだけで、一発リジェクトの危険性が高い。そうなる前に、予防線を張っておくのが無難である。

英語力

とにかく文献を読み、片っ端からいい表現を盗み取るに限る(文自体コピペはいかんけど、あくまで英語表現)。つたない英語表現より、読みやすい表現のほうが絶対アクセプトされやすい。また、経験上、つたない文章は、英文校閲による改善も限界がある。一定ライン(ゼロベンでTOEIC700~800点くらい?)を越えた英文は、校閲者がきちんと意図をくみ取ってくれるので、よりよい英語表現を的確にSuggestしてくれる(これは自分の英語力の改善と共にひしひしと感じている)。結局自分自身の英語能力を上げるのは必要不可欠で、常日頃から意識して鍛えるべき能力という他ない。

SPONSOR LINK

行列操作

Rの行列操作になれる

Rの行列指定になれると、いろいろと作業が簡略化できます。備忘録もかねて、整理したいと思います。

数字による行列指定

Rでは「行列[行番号,列番号]」の形で行列が指定されます。例として、3×3の行列を作り、特定のデータを行列の番号で指定してみます(Xは行列ではなく、データフレームでも可)。

> X <- matrix(rnorm(9,0,1),3,3)
> X
          [,1]       [,2]       [,3]
[1,] 2.3587491 -0.8524698  1.8697275
[2,] 0.8162371  0.4589377 -0.3638491
[3,] 1.8183986  1.2378938 -0.7109136

> X[2,3]
[1] -0.3638491

複数のデータを同時に指定する場合は、”1:3″(1から3まで)や”c(1,3)”(1と3を指定)などを使います。また、空欄”[1,]”とすると、その空欄の部分(行もしくは列)のすべての要素が指定されます。

> X[c(1,2),1]
[1] 2.3587491 0.8162371

> X[1:3,1]
[1] 2.3587491 0.8162371 1.8183986

> X[1,c(1,2)]
[1]  2.3587491 -0.8524698

> X[1,]
[1]  2.3587491 -0.8524698  1.8697275

特定の行や列を取り除きたいときは、数字を負にします。

> X[,-2]
          [,1]       [,2]
[1,] 2.3587491  1.8697275
[2,] 0.8162371 -0.3638491
[3,] 1.8183986 -0.7109136

> X[,c(-2,-3)]
[1] 2.3587491 0.8162371 1.8183986

要素名による行列指定

要素名を指定することも可能です。

> colnames(X)<-c("a","b","c")
> X
             a          b          c
[1,] 2.3587491 -0.8524698  1.8697275
[2,] 0.8162371  0.4589377 -0.3638491
[3,] 1.8183986  1.2378938 -0.7109136

> X[,"a"]
[1] 2.3587491 0.8162371 1.8183986

複数ある行のうち、特定の列名(行名)に一致する列番号を取り出したい場合は、colnames(rownames)を使います。複数の列名(行名)に一致する番号を知りたい場合は、”==”ではなく”%in%”を使います。

> colnames(X)
[1] "a" "b" "c"

> which(colnames(X)=="a")
[1] 1

> which(colnames(X) %in% c("a","b"))
[1] 1 2

> X[,which(colnames(X) %in% c("a","b"))]
             a          b
[1,] 2.3587491 -0.8524698
[2,] 0.8162371  0.4589377
[3,] 1.8183986  1.2378938

> X[,-which(colnames(X) %in% c("a","b"))]
[1]  1.8697275 -0.3638491 -0.7109136

複数ある行のうち、特定の列名(行名)に一致する列番号だけをのぞきたい場合は、”!=”も使えます。ただし、これは指定する行もしくは列が一つだけの場合に限られるので、「-which(colnames(X) %in% c(“a”,”b”))」などとするほうが汎用性はあるように思います。

SPONSOR LINK

function関数の使い方

function関数

Rで複数の関数を組み合わせて値を算出しなければならないとき、function関数で定義してしまうと楽です。しかし、いったい何をやってるのかわかりにくいところもあります。ここでは、function関数を使った関数の定義について簡単に紹介します。

変動係数

変動係数(CV)とは、標準偏差を平均で除したもので、変動性の比較などでよく使われる指標です。RのデフォルトでCVは定義されていないのですが、比較的よく使うものなので、例として挙げてみます。まず、function関数の仕組みですが、引数の指定と、その引数を使った計算式の指定、という形になっています。

CV <- function(x){ sd(x)/mean(x)}

#example
> x <- rnorm(100,10,10)
> CV(x)
[1] 1.042615

上記のfunction()のかっこ内に引数(今回はx)として何を使うのか指定し、{}内に引数を使った計算式を定義します。

ベクター中の0要素をカウントする関数

少し込み入ったものとして、ベクター中に0がいくつかあるのかカウントし、そして0が何番目の要素かを返す関数を定義してみます。

zerocount <- function(x){
 n_zero <- length(which(x==0))
 id_zero <- which(x==0)
 return(list(n_zero=n_zero, id_zero=id_zero))
}

#example
>x <- rpois(100, 2)
> zerocount(x)
$n_zero
[1] 17

$id_zero
 [1] 12 16 21 22 25 31 35 40 43 67 74 77 78 79 86 87 95
 

こちらでは、返す要素が二つあるので、返り値をリスト化して返すよう指定しています(4行目)。

応用

発展的には、乱数発生などのプロセスを組み合わせることで、簡単なシミュレーションモデルを作ることができます。そのほか、自作の作図関数なども作ることができます。興味のある人はぜひトライしてみてください。

最頻値の推定関数

最頻値

平均、中央値、分散など、確率変数を特徴づける指標はいくつもあるけれども、多くのものについてはRで標準実装されている。しかし、地味に知りたい「最頻値」なるものは実装されていない。

最頻値の推定の仕方

離散値の場合は以下のコードで簡単に推定できる。

x <- rpois(100,10)
names(which.max(table(x)))

しかし、サンプルサイズが比較的小さい場合や、確率変数が連続値の場合にはとてもいい方法とは思えない。そこで、table関数の代わりに、density関数をつかった関数を定義する。

mfv <- function(x,from=min(x,na.rm=T),to=max(x,na.rm = T)){
  if(is.numeric(x)==F)stop("use numeric")
  tmp <- density(x,from=from,to=to)
  Mode <- tmp$x[which.max(tmp$y)]
  return(Mode)
}

# example
> x <- rnorm(1000,10,20)
> mfv(x)
[1] 12.16744

density関数は、確率変数の確率密度(のようなもの;実際はカーネル)を変数の傾度にそって推定してくれるもの。上で定義したものは、このカーネルの値が最も高くなる確率変数xの値を返す、という関数になる。細かいけど、いちいち計算するのは面倒なので、関数化してしまうのが楽。。。ちなみに、fromとtoに値を指定すると、その範囲の中だけでカーネルを推定する。例えば、カウントのような正の値しかとらない変数の最頻値を知りたいときは、from=0とすることで、正の実数の範囲の中だけで最頻値を計算してくれる。また、確率など0-1しかとらない、とわかっているのであれば、from=0、to=1とする。

だが、そもそもdensity関数がなにをどうやって推定しているのかようわかっていないので、正確ではない可能性がある。確認次第追記します。

「役立てるだけの能力がなかった」

ベイズの定理

統計の基本定理の一つである「ベイズの定理」が発表されたのは1763年。いまではGoogle検索のアルゴリズムにも採用されるほど重宝されているが(検索履歴を事前分布として使っているものと思われる;詳細は不明)、実際に実用化が進んだのは早くとも2000年以降だと思われる。ベイズの定理が発表された当時は、膨大な計算量が必要とされることから、「役に立たない」ものとみなされていたらしい。ベイズの定理の後に、Fisherにより最尤法なる画期的な推定手法が発表されると、その簡便性と正確性の高さから広い支持を得た。長い長い間、ベイズの定理が表立って利用されることはほとんどなかった。

マルコフ連鎖モンテカルロ

しかし、ごく最近になって、マルコフ連鎖モンテカルロなる逐次的な計算手法が開発された。これは、ベイズの定理と至極相性がよい。というのも、ベイズの定理における計算量の問題が、この計算手法によってかなり緩和されるからだ(※)。もちろん、手計算で行える範疇ではないが、近年のPCの計算能力をもってすれば、十分計算できるほどのものになる。「役に立たない」と思われていたベイズの定理に、やっと人類が追いついた。計算量の問題さえ切抜けてしまえば、ベイズのその他もろもろの利点が浮き彫りになり、各種研究分野で実用化が進んでいる。Google検索におけるベイズ推定の応用もそのひとつになるのだろう。

「今」役に立つ

「役に立つ」という言葉は、今目に見えている問題を相手にとった言葉なのだろう。役に立つ研究というものが、今後起こりうる未知の問題に対し、どれだけの示唆を与えてくれるのだろうか、と思う(「役に立つ研究」は、今の問題には役に立つけど、それっきりなんだろうなぁ)。「質の高い、けど今は役に立たない基礎研究」をサポートする姿勢というは、将来起こりうる不測の事態に対する先行投資とも思うのだけれど。税金を使っている限り、それは過ぎたお願いだろうか。参照:「なんの役に立つんですか?の暴力性

※ ベイズの定理の計算には、パラメーターpが任意の値をとるすべての事象の数え上げが必要になる。このため、普通の計算手法では、PCをもってしても不可能なほどの計算が必要になる。

ナダルが全仏で10回優勝したすごさについて

ナダルがテニスの全仏オープンで10回目の優勝を達成しました。すごすぎて、このすごさを表現するのに適切な言葉が見つかりません。なので、確率で表してみたいと思います。

指数分布のFitting

全仏オープンといえば、テニスの4大大会のひとつに数えられる大会。出場する選手は全力でここにスケジュールを調整してくるし、基本的にトップ選手はみんなエントリーする。なので、世界で一番優勝の難しい大会の一つといえる。ナダルはその大会で10回目の優勝を達成した。これが全人未踏なのはもちろんなのだが、どれくらいの確率で起きうるのかが気になる。なんとか推定できないか。そこでウィキペディアをググると、なんと全仏オープンの歴代優勝者があるではないか。優勝経験者のみを対象に、ナダルの偉業がどれほど「偉業」なのか推定してみる。横軸にある個人の優勝回数をとり、縦軸にそれぞれの優勝回数を達成した人数をとる(ヒストグラム)。これに対して、指数分布をあてはめる;

f(x) = θexp(-θx)
(xは優勝回数)

どれくらいの確率でおきるのか

指数分布でFittingすると、θ = 0.55と推定され、以下のような図になる。縦軸は確率密度で、値が高いほど対応するx軸の事象が生じる可能性が高い。

よく聞く名前では、フェデラー、ジョコビッチは1回、ヴィランデルは3回、ボルグが6回優勝。ナダルの「10回」という優勝回数が生じる確率は、0.004、つまり0.4%だった。上の図でいうと、(ほとんど見えないけど)赤で塗りつぶされているところの面積と対応する。ちなみにこれは、「全仏優勝経験者」が10回優勝する確率。なので、1000年全仏オープンを続けたら、確率論的には4人10回優勝する人がでるかもね、っていうくらい(実際にはこれよりもずっと低いと思う;選手生命の長さなども考えれば)。地殻変動レベルの現象。これはすごい。ナダル、おめでとう。

負の二項分布モデルの精度検証:glmer.nb、glmmamb、JAGS

精度はいかほどですか

この前紹介したけれども、lme4パッケージの更新で、ランダム効果入りの負の二項分布モデルが扱えるようになっていた。しかし、他の関数(glmmadmbやJAGS)と比べてどの程度の精度なのか、気になるところ。そこで、簡単に推定の精度比較をしてみる。モデルは極めて単純な線形モデル。

Y ~ NB(pi, n)
pi = n/(mui + n)
log(mui) = a + bXi + εj
εj ~ Normal(0,σR2)

まず、上記を真のモデルとするシミュレーションデータを生成し、それをglmer.nbなどに推定させる、という作業を100回繰り返す。100回分の推定値の頻度分布と真の値を比べるということをした。パラメーターの真の値はそれぞれa=1.0, b=0.5, σR = 0.7とした。

結果

赤い曲線が100回分の推定値の頻度分布。縦の赤い実線は100回分の推定値の中央値、グレーの実線は真の値。左から切片(a)、回帰係数(b)、ランダム効果のばらつき(σR)。

glmer.nb

glmmadmb

JAGS

こうしてみると、推定精度の観点から一番難がありそうなのはglmmadmb。特にランダム効果のばらつきの推定値は、右側の裾が随分と重い。これはランダム効果が二つ以上になったら絶望的ではなかろうか。glmer.nbは比較的安定していたが、回帰係数の推定値のピークが平らになっているのがやや不安。中央値は真の値に近いけれども。JAGSは可もなく不可もなく、推定としては一番マシかな、という気がする。ランダム効果が増えたときが気になるけれども(たぶんJAGSの推定がかなりまマシなものになる)、面倒なのでやらない。でも、最尤推定の負の二項分布モデルで、ランダム効果の入れ過ぎはよくないということは間違いなさそう。確認してないけど。

シミュレーションデータ生成に使ったコード

glmer.nbの精度検証に使ったコード。そのほかは割愛。
apply系を使えばfor構文が要らなくなる気がするけど。。。

for(i in 1:100){
  #simulate 100 data points from 10 groups with one continuous fixed effect variable
  x<-runif(100,0,1)
  f<-gl(n = 10,k = 10)
  data<-data.frame(x=x,f=f)
  modmat<-model.matrix(~x,data)
  
  #the fixed effect coefficient
  p <- NULL
  p[1] <- 1; p[2] <- 0.5
  fixed<-c(p[1],p[2])
  
  #the random effect
  p[3] <- 0.7
  rnd<-rnorm(10,0,p[3])
  
  #the simulated response values
  size <- 3
  mu <- exp(modmat%*%fixed+rnd[f])
  prob <- size/(mu+size)
  data$y<-rnbinom(100,size=size,prob=prob)

  fit <- glmer.nb(y~x+(1|f),data)
  est[[1]][i] <- fit@beta[1]
  est[[2]][i] <- fit@beta[2]
  est[[3]][i] <- fit@theta
}

スーパー堤防

父の実家が陸前高田なので、節目で被災地に訪れている(特にボランティアというわけではないです)。そのときに聞いた興味深い話。

被災以降、津波の「減災」を目的としたスーパー堤防がガンガン作られており、それはもうすごい勢い。生態学の観点からそれが好ましくないことはもちろんそうなのだが、現地に生きる人の感情を考えると極めて難しい問題だなぁ、と感じていた。しかし、それは大きな勘違い(みんな堤防ほしいだろう、という先入観)で、現地からでる「堤防」に対する意見は住民の間でも大きく割れており、そしてそのパターンがどうにも奇妙だった。

奇妙だったというのは、海の近くに住んでいる住人(つまり被災リスクが高い)ほど堤防反対の意見をもつことが多かったこと。逆に、海から離れた場所に住む人(被災リスクは相対的に低いと思われる)ほど堤防建設を熱望しているという。

海際の人は年配の人が多く住んでいる可能性あるので、年齢層との偽相関は拭えないが、海際の住人の共通した意見として、海に住んでいるのだからそういったリスクは承知の上、海を見れないことのほうがよっぽどつらいし、仕事もしにくい、とのこと(もちろん、これが一般的な傾向というつもりはございません。あくまで聞いた範疇です)。

まったく定量的な調査ではないですが、現地の人に聞く限り、上述のような傾向があるそうです。だれかきちんと証明しているのかな、これ。

さまざまな情報量基準

生態学の界隈では、情報量基準なるものが大人気です。そのくせ、(筆者も含め)それぞれの情報量基準なるものがどんな意味をもつのか、よくわかっていない場合がほとんどです。わかっている範疇で整理。

AIC(赤池情報量基準)

情報量基準といったらAICと言わんばかりですが、なんだったら決定係数との違いすらよく分かっていないかったのが博士課程のころの僕です。AICが低いモデルでは、「同じ母集団から同じサンプル数の別のデータセットを取り出してきたときにも、それなりのパフォーマンス(当てはまり)を発揮する」ということになる。だが、これはサンプル数を無限大に増やしたとしても、必ずしも「真のモデル」には収束しないらしい。あくまで評価基準が予測精度だからだろうか。

AICc

サンプルサイズが小さいときはAICc!みたいな使われ方をしている。だが実は、正規分布を仮定した補正をしているので、正規分布以外のモデルに対して適用するのはご法度らしい。Burnhamさんの本を参照のこと。

BIC(ベイズ情報量基準)

AICと異なり、こちらは真のモデルを選ぶらしい(あくまでサンプルが無限大になったとき、ということだろうか)。あまり使ったことはないが、たしかにAICとは全く異なるモデルを選ぶ傾向がある。予測精度ではなく、真に重要な変数を選び出したい場合にはBICを使った方がよいのだろう。式的には、AICの罰則項にサンプルサイズが入っていないのに対し、BICには入っている。そのため、サンプルサイズが大きくなったときに、前者は対数尤度の項の影響が相対的に大きくなり、より複雑なモデルを選びがちになるようだ。

DIC(Deviance Information Criteria)

こちらはベイズモデルなどでよく使われる情報量基準。だが、広く使われている割に、その数学的根拠は随分と希薄らしい(論文でこういった表記をよく見かける)。DICの式の構造としては、モデルの当てはまりに関わる指標(Deviance)に対し、パラメーター数で罰則を与えている。AICでは2k(kはパラメーター数)という罰則項があるが、これには数学的根拠がある。DICの場合は、そういった裏付けがない、ということだろうか。時間のあるときに調べよう。
また、DIC利用の際の注意点として、パラメーターの正規性が担保された状態でしか役に立たないらしい。詳細はよくわからないが、ベイズのご利益がでる階層線形モデルや状態空間モデルでは、この「パラメーターの正規性」がほとんど担保されないらしい。パラメーターの正規性とはなんぞや。。。(おそらく、パラメーターの事後分布が正規分布で近似できる、ということか)。だとしたらDICを使える場面などそうそうないのでは。

wAIC(widely applicable information criteria)

最近考案されたらしいwAICなるものは、階層をもった複雑なモデルにおいても適用できるらしい。これまたAICに引き続き日本人が発見したので驚き。日本人って統計につおいのだろうか。開発者のワタナベさんのページはこちら(ワタナベさんが開発したからwAICかと思った。違うのかorz)。内容をまったく理解していないのですが、浅はかな理解では予測精度を問題にしているような気がするので、モデル選択の結果の解釈としてはAICと同様になるような気がする。階層モデルのようなものにも適用できる点が違いか(あくまでユーザー視点からすれば)。もう少し勉強する。

lme4パッケージの更新

しばらく見ないうちに、最新版のlme4パッケージでは劇的な進化を遂げていたようなので、それについてメモ。

glmer.nbの出現

これまで負の二項分布を仮定したモデルで、ランダム効果をいれるとしたら、glmmADMBパッケージか、もしくはベイズモデル化(WinBUGSなど)する必要があった(ADMBもベイズみたいなもん?中身はよくわからない)。しかし、最新版のlme4パッケージでは、glmer.nbという関数が表れたことにより、負の二項分布でもランダム効果を簡単に入れることができるようになった。仕様については基本的なlmerと変わらない模様。ただし、頻度主義的に扱うには尤度関数が複雑すぎる気がするので、いったいどこまで結果が信用できるかは不明。>3のランダム項をいれるのはきっとご法度。

mcmcsamp()の廃止

昔のlme4には、最尤法ではなく、MCMCでベイズっぽく推定するための関数が用意されていた。しかし、最新版のlme4パッケージでは除外された模様。Rのhelpには以下のように書いてある。
Due to difficulty in constructing a version of mcmcsamp that was reliable even in cases where the estimated random effect variances were near zero (e.g. https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q4/003115.html), mcmcsamp has been withdrawn (or more precisely, not updated to work with lme4 versions >=1.0.0)

どうやらランダム効果のばらつきがゼロ付近のときに、あまりいい推定をしてくれないということで、やめにしたということらしい。個人的には、ランダム効果が多いとき(>3)の最尤推定は不安なので、mcmcsamp()の取り下げは手痛い、、、JAGSとかはコード書くのが面倒だからなぁ、、、

lmerにおけるREMLのデフォルト化

REMLとはなんぞや、と調べてみると、REstricted Maximum Likelihoodの略らしい。最新版のlmerでは、MLではなくREMLがデフォルトになっている。REMLとMLの違いは、(筆者の浅はかな理解によると)、前者は固定効果による自由度減少を考慮したうえでランダム効果の分散を推定し、後者は考慮しないらしい。つまり、後者の場合、ランダム効果のばらつきの推定は偏りが生じているということになる(不偏推定量でない?)。ただ、REMLではモデル選択ができないため(なぜかはよくわからない)、MuMInパッケージなどでモデル選択する場合はREML=Fとする必要がある。正規分布以外ではMLしか使えないので心配の必要なし。

glmerとlmerの(よくわからない)分離

正規分布の場合はlmer、それ以外はglmerということでよさそう。