No.13988 non-selectiveになる場合のGLM  【GLM躓き者】 2010/12/17(Fri) 00:56

餌A,Bの二者択一で,ある動物の成長に伴う,餌選好性の変化を評価したいと思っています。
結果が,小さいサイズのうちはAを極端に好むが,成熟すると選好性がなくなり,A,Bそれぞれの摂食割合が50%くらいに漸近するとします。
この場合,個体重と各個体の摂食割合の関係をプロットすると,“重量の増減に伴って 1 または 0 に近づく漸近線”ではなく,“ 1 または 0.5 に近づく漸近線”のようなグラフになると思います。
そのようなグラフに対してもGLM(あるいはGLMMなど)をそのまま当てはめてしまってもいいものでしょうか?
実際にGLMをそのまま適用してみると,そのグラフは,成長に伴って0.5に漸近はせず,0方向にテイリングするような感じになり,当てはまりが良くないようです。

よろしくお願いします。

No.13989 Re: non-selectiveになる場合のGLM  【青木繁伸】 2010/12/17(Fri) 08:50

GLM,GLMと仰いますが,family は何を指定したのですか。また,簡単なテストデータとその分析結果も示した方が,ほかの人に分かりやすいと思います。

No.13995 Re: non-selectiveになる場合のGLM  【GLM躓き者】 2010/12/18(Sat) 03:56

大変,失礼しました。
以下のような生データです。
weight	preyA	preyB	sum
10 10 0 10
15 12 0 12
19 20 1 21
21 21 2 23
25 30 4 34
30 32 7 39
32 29 10 39
33 31 15 46
35 28 21 49
38 27 24 51
41 30 25 55
42 29 27 56
44 30 28 58
50 32 30 62
餌はカウントデータで,family=binomialを当てています。
> model <- glm(cbind(preyA, sum-preyA) ~ weight, family=binomial, data=data01)
> summary(model)

Call:
glm(formula = cbind(preyA, sum - preyA) ~ weight, family = binomial,
data = data01)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.6861 -0.4585 0.3545 0.9439 2.0676

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.68739 0.45945 8.026 1.01e-15 ***
weight -0.08300 0.01183 -7.018 2.26e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 75.874 on 13 degrees of freedom
Residual deviance: 15.882 on 12 degrees of freedom
AIC: 65.908
グラフは以下のように当てはまりがよくありません。
おそらく二項分布(binomial)を指定しているのがまずいと思うのですが,代わりに何を指定すればいいかわかりません。

よろしくお願いします。


No.13996 Re: non-selectiveになる場合のGLM  【GLM躓き者】 2010/12/18(Sat) 04:02

生データについて説明不足でした。
weightが個体重,preyがその個体が食べた餌の数,sumが食べた餌の総数
実験設定としては,preyA, Bそれぞれ同数づつ与えたとき,どちらをどれだけ食べたかを測定した,ということです。
weight preyA preyB sum
10 10 0 10
15 12 0 12
19 20 1 21
21 21 2 23
25 30 4 34
30 32 7 39
32 29 10 39
33 31 15 46
35 28 21 49
38 27 24 51
41 30 25 55
42 29 27 56
44 30 28 58
50 32 30 62

No.13998 Re: non-selectiveになる場合のGLM  【青木繁伸】 2010/12/18(Sat) 10:47

まず第一に,あなたのモデルに family=binomial を指定するのは不適切でしょう。binomial は,生存率曲線の場合のように,独立変数によって対象中の何%が故障(死亡など)するかというような現象を記述するものです。あなたのデータで,「餌の個 数」という記述があるのですが,たぶん餌がペレット状になっているから「個数」で測定したのでしょうけど,その1個のペレットの重量は採用した餌では一定 なのでしょうけど,例えば別のメーカーのペレットは1.5倍の重量だった等という場合には,研究者によって「個数」は異なってくるでしょう。もし,餌が粉 末状で,どちらの餌も摂取量がグラムで測られたとするとどうなりますか。もはや「個数ではない」ので,だれも binomial を指定しようとは思わないでしょう。しかるに,どのような餌の形態であろうと,測定によって得られる分析結果(生長につれ,餌の好みはなくなる)というこ とに変わりないでしょう。
つまり,餌の数は離散変数(計数値)のようには見えたけど,飛び飛びの値をとるものの,本当は連続変数なのです。

では,どうすればよいか。
まずはモデル式を立てる必要があるでしょう。示されたデータに合うような数学曲線はいくつかありますけど,このデータがどのような理論曲線に従うかは私には不明です。数式が決まれば,非線形最小二乗法でパラメータを求めることで解決するでしょう。

No.13999 Re: non-selectiveになる場合のGLM  【青木繁伸】 2010/12/18(Sat) 11:19

ただ単に似ているからと言うだけでロジスティック曲線を当てはめてみると以下のようになった。
> d <- read.table("data01.R", header=TRUE)
> d$y <- d$preyA/d$sum
> ans <- nls(y ~ 0.5/(1+a*exp(-b*weight))+c, d, start=list(a=-0.5, b=-0.06, c=0.5))
> p <- ans$m$getPar()
> (a <- p[1])
a
0.000284451
> (b <- p[2])
b
-0.262103
> (c <- p[3])
c
0.4885675
> x <- seq(0, 60, length=1000)
> y <- 0.5/(1+a*exp(-b*x))+c
> plot(x, y, type="l")
> points(d$weight, d$y, pch=19, col="maroon2", xpd=TRUE)


No.14000 Re: non-selectiveになる場合のGLM  【GLM躓き者】 2010/12/18(Sat) 15:35

明快なご回答,ありがとうございました。
餌の場合,重量換算と考えてbinomialを使わない手がある,ということ,よく理解しました。

あえてカウントデータとして扱うことは,難しいでしょうか?
餌の選択が栄養価や,重量によるものだと想定すれば,連続変数ですが,preyA が赤色,preyB が黒色のペレットで,動物は餌の色をひとつひとつ見て摂食判断をする場合(栄養価に関係なく),やはりカウントデータとして扱うのがふさわしいように思えます。

重ね重ねの質問になりますが,もし何かご教授いただければ幸いです。

No.14001 Re: non-selectiveになる場合のGLM  【竹澤】 2010/12/20(Mon) 09:04

ノンパラメトリック一般化線形回帰を使うと,例えば以下のようになります。
function() {
library(mgcv)
weight <- c(10,15,19,21,25,30,32,33,35,38,41,42,44,50)
preyA <- c(10,12,20,21,30,32,29,31,28,27,30,29,30,32)
preyB <- c(0,0,1,2,4,7,10,15,21,24,25,27,28,30)
sum <- c(10,12,21,23,34,39,39,46,49,51,55,56,58,62)
gam123 <- gam(cbind(preyA, preyB)~s(weight), family=binomial)
print(summary(gam123))
ey <- predict(gam123, type="response")
plot(weight, preyA/sum, ylim=c(0.4,1.1))
lines(weight, ey)
}
回帰曲線は,正規分布を仮定してロジスティック曲線をあてはめた場合とほとんど違いません。しかし,信頼区間などにはかなりの違いが現れるだろうと思います。

No.14005 Re: non-selectiveになる場合のGLM  【知ったかぶり】 2010/12/20(Mon) 16:31

link関数に手を加えれば,
>“ 1 または 0.5 に近づく漸近線”
となるようにすることができます.glm()のlink関数は,自分で定義できるらしいのですが,よくわからないので,optim()を使用.
LL2 <- function(a){ # 対数尤度関数を定義
lb <- 1 / (1 + exp(a[3])) # preyAの比率の下限
q <- lb / (1 + exp( - (a[1] + a[2] * weight))) + (1 - lb)
sum(lchoose(sum, preyA) + preyA * log(q) + (sum - preyA) * log(1 - q))
}
model2 <- optim(c(1, -0.1, -0.1), LL2, control=list(fnscale=-1)) # 対数尤度を最大化するパラメータの推定
pa2 <- model2$par # パラメータ推定値
-2 * model2$value + 2 * length(pa2) # AIC
xx <- seq(5, 55, 1)
plot(weight, preyA / sum, ylim=c(0, 1), xlim=c(0, 55))
lines(xx, (1 / (1 + exp(pa2[3]))) / (1 + exp( - (pa2[1] + pa2[2] * xx))) + (1 - (1 / (1 + exp(pa2[3])))), col=2)
パラメータの数は1つ増えますが,AICは普通にlogitリンクを使った場合より小さくなります.

No.14013 Re: non-selectiveになる場合のGLM  【GLM躓き者】 2010/12/21(Tue) 16:52

みなさん,ご回答ありがとうございます。

竹澤 様
ノンパラメトリックは確かによさそうな気がします。
今まで耳にはしておきながら手を出していませんでしたが,勉強します。

知ったかぶり 様
比率の下限を設定する方法があるとは知りませんでした。
下限の設定をどう決めるかの根拠が難しそうですが,場合によってはこれも有効に思えます。

改めまして,貴重なお時間いただきありがとうございました。

● 「統計学関連なんでもあり」の過去ログ--- 044 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る