（旧サイトより移行）連続量を単位あたりに直してGLMする場合の対処：gaussian("log")、Gamma("log")、offsetなどなど）

# 2014.02.20 旧サイトを閉じるため、移植しました。
# なお、下記のようなケースではゴンペルツ曲線などを用いた成長モデルを使った方がいいかなと最近思ってます。

最近、実験系の統計を引き受けていて気になったので、連続量を単位あたりでGLM推定する方法について検証してみました。

単位あたりの量をGLMで解析する時、個体数や頻度などの整数の場合の対処はoffset項の利用で解決するのが常套手段ですね（密度や○×率のような割り算した値ではなく、元の値はいじらずに単位量を係数1として説明変数に加える。当ブログでもかつて紹介：http://nhkuma269.blog77.fc2.com/blog-entry-9.html）。
例えば、同じく個体密度0.1でも、1/10と100/1000とでは意味合いが違うが、割り算すると両者は同一密度として扱われてしまう。前者では1個体の増減が大きな誤差を生むが、後者ではほとんど影響なし。

では、元々が連続量のものを単位あたりの量にする場合も同様にoffset項による対処がよいのだろうか？例えば、サンゴの枝あたりのクロロフィル量とか、成長量（実験後重量 / 実験前重量）のような量は、必ず単位量（この場合は枝）のバラツキによる誤差が出てしまう（なるべく条件は揃えるように努力はしているだろうが）。これを割り算してしまうと誤差が統計結果に悪影響するだろう。

CrawleyのR統計本によると、単位量は共変量として説明変数に加えるか、log(目的の量 / 単位量)というように、割り算した上で対数変換せよ。と書かれていた。前者は要因の1つとしてカウントするならばそれでいいとして、後者はすでに古典的な対処法でしょう…（本自体がすでに一時代前のものなので仕方ないですが）。

ついでに、論文などで時々見掛ける、family=gaussian(link="log")、正規分布でリンク関数が対数という一見すると対数正規分布っぽいけれど違うらしい（対数正規分布の分散は平均と共に増大するようだ→平均や分散の式をチェックされたし）、こいつの性質もチェックしてみます。もう一つおまけに、古典統計で対数変換する際によく使うlog(x + 1)変換もチェックしてみました。ただし、これはlog(0)が計算できない問題を回避する目的に限定します。


比較するモデルは下記の6バージョン。これを応答変数が対数正規分布とガンマ分布の場合のそれぞれについて一次回帰式のパラメータ推定（切片と回帰係数）をします：
（色は下記の図とおおよそ対応）
m11：正規分布モデル、割合にした上でlog(x + 1)の対数変換という古典統計の常套手段（グラフ上では黒ライン）
m12：正規分布モデル、割合にした上でlog(x + 1e-10)変換、1の代わりにごく小さい小数(1.0 x 10^(-10))でlog(0)を回避
m13：正規分布モデル、割合にしているが、変数変換はせず、連結関数を対数に指定（以下、log(0)回避には1e-10を使用）
m14：正規分布モデル、単位量をoffset(log(off))として係数1の説明変数に加え、連結関数を対数に指定
m23：ガンマ分布モデル、割合にしているが、変数変換はせず、連結関数を対数に指定（ガンマ分布は正の値しか取れないので、1と2の変数変換のモデルは作れない）
m24：ガンマ分布モデル、単位量をoffset(log(off))として係数1の説明変数に加え、連結関数を対数に指定


結論から言うと、ケースバイケースな複雑な結果になりました…。実際の利用では、m13とm23、またはm14とm24をAICで比較するのがよいでしょう。
・やはりデータの対数変換は止めた方がよい、連結関数：logを用いるべきです。対数変換したモデルでは推定が大きく乱れる場合があり信頼できない（図1_1, 1_2のm11黒、m12赤）。
・データが対数正規分布の時、割合の正規分布モデル（m13緑）は、offsetの正規分布モデル（m14青）よりも推定の分散が小さかった（図1_1、2_1）。ただし、データがガンマ分布の時は、m14の方がm13よりも推定の分散が小さかった（図1_2、2_2）。両者の違いはさほど大きくなかった。
・ガンマ分布モデルでは、割合（m23水色）と、offset（m24ピンク）とは推定値がまったく同じ！（そのため、m23の水色は完全にマスクされている）
・つまり、この数値実験の限りでは、予想に反して連続量を割合や比率にした量をGLMで解析する際、割り算をした値を用いても大した問題はないことになる。offsetの使用による論文中での説明のややこしさや、記述可能なモデルの可塑性などを考えると…割り算していい気がしてきました。しかし、なぜ割り算にしても大丈夫だったのか、背景にある数学的なロジックは朧気なままです…たぶん、元の値と割算値とで確率分布が変わらないからだと思います（よくある整数値の割り算の場合は、小数点を取るようになるので本来のポアソンや二項分布が使用不可になるが）。例外があるとすれば、単位量の方にも確率誤差（その値が属している確率分布からのズレ）がある場合でしょう。その場合はoffsetすべきということかと。
・ちなみに、常套手段なlog(x + 1)変換（m11黒）の推定は危険です。とくに平均値の小さい推定では危険（図1_1、1_2）。log(0)の回避には1の代わりにごく小さな値（1e-10くらい）を足すのがよいでしょう。ただし、古典統計で整数を対数変換する場合は、平均と分散の関係を調整する意図があるので + 1 のままで）。

なお、ここでやっているのは入り口と出口がちゃんと一致するかを確認するだけの数値実験に過ぎません。参考にする際には自己責任でお願いします。



# 以下、こんな関数でGLM推定を1000回繰り返し、パラメータ推定の精度をチェックしました。
# 平均 exp(alpha + beta*X)という一次回帰を考え、こちらで指定したalphaとbetaを用いて、対数正規分布とガンマ分布のデータを発生させ、各モデルによってalphaとbetaを逆推定し、初めに指定した値を再現できるかどうかチェックします。
# alpha, beta, distributionは1が対数正規分布、2がガンマ分布、n.itrは繰り返し数を表します。

cont.skew <- function(alpha, beta, distribution, n.itr) {
estim <- numeric(0) # 後でデータをくっつけるためのイントロンみたいなもの
set.seed(1)
X <- rep(c(1:10)*0.1, each=10)
off <- runif(length(X), min=1, max=10) # 下記のoffsetで用います。数字に特に意味はなし
mean <- exp(alpha + beta*X)*off # 平均：log(Y/off) ~ alpha + beta*X、なので。
sd <- exp(0.5)
for (n in 1:n.itr) {
Y <- rbind(rlnorm(100, meanlog=log(mean), sdlog=log(sd)), # 対数正規分布
rgamma(100, shape = mean^2/sd^2, scale = sd^2/mean)) # ガンマ分布（cf. rgamma()のhelp）
Y <- Y[distribution, ]
D <- data.frame(X, Y, off)
m11 <- glm(log(Y/off + 1) ~ X, family=gaussian(link="identity"), D) # 正規分布モデル（割合、対数変換1）
m12 <- glm(log(Y/off + 1e-10) ~ X, family=gaussian(link="identity"), D) # 正規分布モデル（割合、対数変換2）
m13 <- glm((Y/off + 1e-10) ~ X, family=gaussian(link="log"), D) # 正規分布モデル（割合、連結関数：対数）
m14 <- glm((Y + 1e-10) ~ X + offset(log(off)), family=gaussian(link="log"), D) # 正規分布モデル（オフセット、連結関数：対数）
m23 <- glm((Y/off + 1e-10) ~ X, family=Gamma(link="log"), D) # ガンマ分布モデル（割合、連結関数：対数）
m24 <- glm((Y + 1e-10) ~ X + offset(log(off)), family=Gamma(link="log"), D) # ガンマ分布モデル（オフセット、連結関数：対数）
estim2 <- c(rbind(coef(m11), coef(m12), coef(m13), coef(m14), coef(m23), coef(m24)))
estim <- rbind(estim, estim2)
} # forループ、ここまで
par(mfcol=c(1,2))
plot(density(estim[,1]), lwd=4, xlim=c(alpha-0.5, alpha+2), ylim=c(0,3), main="alpha (intercept)")
for (i in 2:6) { lines(density(estim[,i]), lwd=4, col=i)
abline(v=alpha) } # alpha値の図示
plot(density(estim[,7]), lwd=4, xlim=c(beta-1, beta+1), ylim=c(0,3), main="beta (coefficient)")
for (j in 2:6) { lines(density(estim[,j+6]), lwd=4, col=j)
abline(v=beta) } # beta値の図示
estim.m <- matrix(apply(estim, 2, mean), ncol=2) # alpha, betaの推定値を計算
colnames(estim.m) <- c("alpha", "beta")
estim.m} # 推定値を表示（cont.skew関数、ここまで）

# 図1_1： 対数正規分布のパラメータ推定、平均値が小さい場合
cont.skew(alpha=-1, beta=1, n.itr=1000, distribution=1)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# m12赤、m23&m24ピンク、m13緑、の順に推定がよかった。m11黒は大きくずれた。
# 推定値（1~6の順にm11～m24）、m23とm24（一番下の2つ）は完全に推定値が一致！
# alpha beta
# [1,] 0.3131626 0.3959152
# [2,] -1.0009638 0.9999089
# [3,] -0.8811454 1.0011294
# [4,] -0.8879904 1.0073653
# [5,] -0.8793282 0.9996438
# [6,] -0.8793282 0.9996438

# 図1_2：ガンマ分布のパラメータ推定、平均値が小さい場合
cont.skew(alpha=-1, beta=1, n.itr=1000, distribution=2)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# m14青、m13緑、m23&m24ピンク、の順に推定がよかった。m11黒とm12赤は大きくずれた。
# なぜかガンマ分布の推定なのにガンマ分布モデル（m23&m24）が最強にならない…

# 図2_1：対数正規分布のパラメータ推定、平均値が大きい場合
cont.skew(alpha=2, beta=1, n.itr=1000, distribution=1)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# 平均値が大きくなると、モデルによる推定の違いは小さくなる（中心極限定理による正規分布への収束）
# m12赤、m23&m24ピンク、m11黒、の順に推定がよかった。

# 図2_2：ガンマ分布のパラメータ推定、平均値が大きい場合
cont.skew(alpha=2, beta=1, n.itr=1000, distribution=2)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# 平均値が大きくなると、モデルによる推定の違いは小さくなる（中心極限定理による正規分布への収束）
# m11黒以外は真の値にほぼ収束