No.21423 Re: 一般化線形モデルのリンク関数について 【青木繁伸】 2014/10/30(Thu) 11:22
丸付き数字は使わないでと，上の方に書いてあります。
さて，おおもとの理論曲線は y = a*b^x だということでしょうか？
従属変数を対数変換して lm を適用するというのは，近似値というか，log(y) の誤差が正規分布ということを仮定するものです。
glm で family=gaussian(link="log") を指定したときに実際にどのように計算されるのかは私にはフォローできませんが，nls で非線形回帰して直接 y = a*b^x の，a, b を求める場合とほとんど同じ結果になります。また，非線形回帰では，当然ですが，従属変数が 0 である場合も何の問題もありません。
ということで，あなたの場合には非線形回帰でよいのでは？
従属変数が 0 になるデータを除いて，lm, glm, nls を試して比較するとよいでしょう。
No.21424 Re: 一般化線形モデルのリンク関数について 【竹田】 2014/10/30(Thu) 12:29
青木先生
注意事項をよく読んでおらず申し訳ありません。
早急なご回答，ありがとうございます。
先行研究内におおもとの理論曲線に関する記述はないのですが，行われている解析としましては，
先生がおっしゃられた"y = a*b^x"の式のbにあたる部分を推定しています。
早速，従属変数が0になるデータを除いて，lm()，glm()，nls()を試してみました。
結果としましては，glm()，nls()がほとんど同じ結果(減少速度？でいうと0.01%以下の差)，lm()が結構違う結果(同じく減少速度？でいうと10％ほど違う)が得られました。
これはnls()，もしくはglm()を使ったほうがよいという解釈になるのでしょうか？
また，続けての質問で申し訳ないのですが，私の研究では「減少速度のグループ差(3グループ)」もみる予定です。
glm()では，グループのカテゴリーを説明変数に入れる共分散分析の形をとればよいと考えているのですが，
nls()の場合は，グループごとに係数を推定する形になるのでしょうか？
大変無知で申し訳ありませんが，ご回答いただければと存じます。
何卒よろしくお願い申し上げます。
No.21425 Re: 一般化線形モデルのリンク関数について 【青木繁伸】 2014/10/30(Thu) 13:06
> これはnls()，もしくはglm()を使ったほうがよいという解釈になるのでしょうか？
答の分かっているデータを使って分析してみるとよいでしょう。
以下のデータは a=20，b=0.85 で作り，小数点以下7桁まで示したもの。x ynls では
0 20.0000000
1 17.0000000
2 14.4500000
3 12.2825000
4 10.4401250
5 8.8741062
6 7.5429903
7 6.4115418
8 5.4498105
9 4.6323389
10 3.9374881
11 3.3468649
12 2.8448351
13 2.4181099
14 2.0553934
15 1.7470844> ans3 <- nls(y ~ a*b^x, start=list(a=20, b=0.85))となる。lm, glm ではいくつになるか。比較すれば結論は明確。更に，y=0 のデータがあれば1つに決まる。（って，前のコメントに書いた）
> summary(ans3)
Parameters:
Estimate Std. Error t value Pr(>|t|)
a 1.99999999954e+01 1.72868447474e-08 1156949130.26 < 2.22e-16
b 8.50000000006e-01 1.94419718026e-10 4371984532.42 < 2.22e-16
> glm()では，グループのカテゴリーを説明変数に入れる共分散分析の形をとればよいと考えているのですが，
グループ別に推定するほうがよいでしょう
No.21426 Re: 一般化線形モデルのリンク関数について 【竹田】 2014/10/30(Thu) 18:10
青木先生
丁寧なご説明ありがとうございます。
実際に手を動かしてみて，ようやく分かりました。
nls()を使おうと思うのですが，今後のためにもう1点教えていただけないでしょうか？
さきほどいただいたコメントの中で
>グループ別に推定するほうがよいでしょう
とおっしゃられましたが，これはどうしてでしょうか？
3グループの比較なので，共分散分析のみだとグループ個々の差までは言えず，
多重比較をしなければいけないからでしょうか？
それともほかの理由，例えばグループ別にモデルを作成したほうが，係数の推定の精度が上がったりするのでしょうか？
的外れな質問でしたら申し訳ありません。
こんなに質問をしてばっかりでよいものかわかりませんが，ご回答のほど何卒よろしくお願い申し上げます。
No.21427 Re: 一般化線形モデルのリンク関数について 【青木繁伸】 2014/10/30(Thu) 18:31
この質問も，実際に答の明らかなテストデータで実際に分析してみるとよいでしょう。
あなたのやろうとしていることで，「グループ別の係数 b は推定できますか？」
No.21428 Re: 一般化線形モデルのリンク関数について 【竹田】 2014/10/30(Thu) 20:53
青木先生
何度もお手数をおかけして申し訳ありません。
先ほどの青木先生のデータに一部追加して，glm()で推定してみました。
x y group
0 20.0000000 A
1 17.0000000 A
2 14.4500000 A
3 12.2825000 A
4 10.4401250 A
5 8.8741062 A
6 7.5429903 A
7 6.4115418 A
8 5.4498105 A
9 4.6323389 A
10 3.9374881 A
11 3.3468649 A
12 2.8448351 A
13 2.4181099 A
14 2.0553934 A
15 1.7470844 A
0 30.0000000 B
1 22.5000000 B
2 17.3250000 B
3 13.3402500 B
4 10.2719925 B
5 7.9094342 B
6 6.0902644 B
7 4.6895036 B
8 3.6109177 B
9 2.7804067 B
10 2.1409131 B
11 1.6485031 B
12 1.2693474 B
13 0.9773975 B
14 0.7525961 B
15 0.5794990 B
0 40.0000000 C
1 36.4000000 C
2 33.1240000 C
3 30.1428400 C
4 27.4299844 C
5 24.9612858 C
6 22.7147701 C
7 20.6704408 C
8 18.8101011 C
9 17.1171920 C
10 15.5766447 C
11 14.1747467 C
12 12.8990195 C
13 11.7381077 C
14 10.6816781 C
15 9.7203270 C
y = a*b^xとすると，グループAは(a,b)=(20,0.85)，Bは(30,0.77)，Cは(40,0.91)で作ってあります。
Rでは以下のように推定しました。
model<-glm(y~x*group,family=gaussian(link="log"),start=c(20,-0.1,1,1,-0.1,-0.1))
その結果，
> model$coefficient
(Intercept) x groupB groupC x:groupB x:groupC
2.99573227 -0.16251893 0.39648943 0.69314718 -0.10337401 0.06820825
> exp(model$coefficient[2])
x
0.85
> exp(model$coefficient[2]+model$coefficient[5])
x
0.7665212
> exp(model3$coefficient[2]+model3$coefficient[6])
x
0.91
となり，おおよそ近い結果が得られました。
切片に関しても同様です。
データの作り方，glm()の使い方，どこがどう間違っているのか分かりません。
本旨とはかけ離れた内容で申し訳ありませんが，頭が混乱しています。
どの点が間違っているのでしょうか？助けてください。
何卒よろしくお願い申し上げます。
No.21429 Re: 一般化線形モデルのリンク関数について 【青木繁伸】 2014/10/30(Thu) 21:03
完璧でしょう。データが真実を語ったのです。
No.21430 Re: 一般化線形モデルのリンク関数について 【竹田】 2014/11/04(Tue) 10:53
青木先生
全ての質問に対しての丁寧なご返答ありがとうございました。
● 「統計学関連なんでもあり」の過去ログ--- 047 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る