「統計学関連なんでもあり」の過去ログ--- 048

No.22674　GLMによるモデルと信頼区間の図　　【ぐれいぷ】　2019/01/17(Thu) 11:31

　一般化線形モデルによる解析による予測値と，その信頼区間の求め方についてお教え願います。

　7本の樹（A〜G）について，ある測定値（measure）と比率（30枚の葉のうち症状の出た葉（A）の割合）に関して，一般化線形モデル（ロジスティック回帰）で解析しています。元のデータとモデルによる曲線，さらに95％信頼区間を図示したいと考えています。実際のデータによるRでの解析は以下の通りです（信頼区間は下側のみ示します）。

　図示した場合，信頼区間の破線がx軸の値の細かさで変動し，4から5で0.001刻みにすると階段状になってしまうのですが，モデルによる曲線と信頼区間を示す適切な方法をご存知でしたらお教えください。よろしくお願いします。

> a
  tree measure  A  N
1    A    4.50 11 30
2    B    4.33 11 30
3    C    4.33 13 30
4    D    4.50 17 30
5    E    4.67 23 30
6    F    4.33 23 30
7    G    5.00 23 30

> 
> glm <- glm(cbind(a$A,a$N-a$A)~a$measure,family=binomial)

> glm$coefficients
(Intercept)   a$measure 
  -8.270453    1.900956 

> ##############x軸の値を0.1ごと
> x <- seq(4,5,by=0.1)
> pred.y <- 1/(1+exp(-(-8.270453+1.900956*x)))
> plot(a$measure,a$A/a$N,xlim=c(4,5),ylim=c(0,1))
> lines(x,pred.y)
> 
> ##############信頼区間を求める
> mean <- length(a$N)/sum(1/a$N)
> mean
[1] 30
> lower <- qbinom(0.025,mean,pred.y)/mean
> lines(x,lower,lty=2)
> 
> 
> ##############x軸の値を0.001ごと
> x <- seq(4,5,by=0.001)
> pred.y <- 1/(1+exp(-(-8.270453+1.900956*x)))
> plot(a$measure,a$A/a$N,xlim=c(4,5),ylim=c(0,1))
> lines(x,pred.y)
> 
> lower <- qbinom(0.025,mean,pred.y)/mean
> lines(x,lower,lty=2)

No.22675　Re: GLMによるモデルと信頼区間の図　　【青木繁伸】　2019/01/17(Thu) 14:05

qbinom(0.025,mean,pred.y)/mean は，qbinom(0.025,mean,pred.y) 自体が離散値なので，それを定数で割ったところで離散値には変わりない。したがって，「階段状に表示される」のが正確で正常な状態です。

No.22676　Re: GLMによるモデルと信頼区間の図　　【ぐれいぷ】　2019/01/17(Thu) 14:39

青木先生
　お忙しい中早速のご回答ありがとうございました。モデルと信頼区間について勉強しなおします。今後ともよろしくお願いします。