GLMの基礎：オフセット項（割算値の回避）

「二項分布を使った GLM」からの続き。

「6.6 割算値の統計モデリングはやめよう」をもとにした。

ここでの割算値とは、打率、人口密度などの「割合」のこと。なんらかの割合同士を比較した分析は、確かに頻繁にみかける。「単純な平均」も「割算値」の例と考えるが、それも多用＆誤用されていると思う（このブログでも「普通の平均ではダメ」等で取り上げた）。

観測データを何らかの「割合」に変換すると情報が失われる

というのがポイントだ。そのような「割算値」を回避する方法が、ここで取り上げる「オフセット項」。


観測データ

以下の観測データは、植物の個体数と観測した地点の明るさと面積で、全 100 箇所。解析の目的は

調査値 i における植物個体の「人口密度」が「明るさ」x_i にどう影響されているか知りたい P.132 

「人口密度」が不自然に思えたので、ここでは「個体密度」とした。

> d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/binomial/data4b.csv")
> summary(d)
       y               x                A        
 Min.   :13.00   Min.   :0.0500   Min.   : 3.50  
 1st Qu.:36.00   1st Qu.:0.4000   1st Qu.: 8.80  
 Median :46.00   Median :0.5150   Median :10.30  
 Mean   :48.09   Mean   :0.5141   Mean   :10.47  
 3rd Qu.:59.50   3rd Qu.:0.6425   3rd Qu.:12.15  
 Max.   :95.00   Max.   :0.9900   Max.   :17.40  

y：植物の個体数、x：観測地点の明るさ、A：観測地点の面積

> library(ggplot2)
> theme_set(theme_gray(base_size=12, base_family="HiraKakuProN-W3"))
> ggplot(d,aes(x=A,y=y,colour=x))+geom_point(size=6)+scale_color_gradient(low = "black",high = "white")+xlab("面積 A")+ylab("個体数 y")+ theme(legend.position=c(.94,.13))+ labs(colour="明るさ x")

このプロット図（本書図 6.10A に相当）から

• 同じ面積なら、明るい（白色に近い点）ほど、個体数が多い

という傾向がうかがえる。


分析

個体密度 = 個体数／面積より、面積が A_i である観測地点 i における個体密度は

平均個体数 λ_i / A_i = 個体密度

個体密度は正の数なので、指数関数と明るさ x_i の依存性を組み合わせて

　λ_i = A_i × 個体密度
   　= A_i exp(β₁ + β₂x_i)
      = exp(β₁ + β₂x_i + logA_i)

線形予測子 z_i = β₁ + β₂x_i + logA_i 、対数リンク関数・ポアソン分布の GLM となる。この logA_i のように、線形予測子の中でパラメーターがつかない項を「オフセット項」。

> glm(data=d, y~x, offset = log(A), family = poisson)

Call:  glm(formula = y ~ x, family = poisson, data = d, offset = log(A))

Coefficients:
(Intercept)            x
     0.9731       1.0383

以下は、glm の結果から predict したもの（末尾 Rコード 参照）。

左は本書図 6.10B に相当し、直線は明るさ x で {0.1, 0.3, 0.5, 0.7, 0.9} 。右は本書にはないが、面積 A を {4, 8, 12, 16, 20} で、明るさを x 軸にしたもの。先の図 6.10A のプロットとの関係では、左の方が分かり易いが、明るさと個数密度の関係は右の曲線の方が分かりやすいかもしれない。

推定パラメーターより平均個体数は

　λ_i = exp(0.9731 + 1.0383x_i)
   　= exp(0.9731) × exp(1.0383x_i)

> exp(1.0383)
[1] 2.824411
> exp(1.0383*0.1)
[1] 1.109412

よって、明るさが 10% 増えるごとに平均個体数は約 1.1 倍増加と予測。

このように、「観測データの割算」をすることなく、面積 A に比例、明るさ x の効果を推定できた。

最後に、先のモデルと「切片だけのモデル」で AIC を比較する。

> glm(data=d, y~x, offset = log(A), family = poisson)
Null Deviance:    261.5 
Residual Deviance: 81.61 AIC: 650.3
> glm(data=d, y~1, offset = log(A), family = poisson)
Null Deviance:    261.5 
Residual Deviance: 261.5 AIC: 828.2

明るさ x を含んだモデルの方が「切片だけのモデル」より、予測の良さを示している。

「目次」へ続く。

R コード

> fit <- glm(data=d, y~x, offset = log(A), family = poisson)
# x = Square A
> A=seq(0,20,length=100)
> y1<-predict(fit, data.frame(A=A,x=0.1),type = "response")
> y3<-predict(fit, data.frame(A=A,x=0.3),type = "response")
> y5<-predict(fit, data.frame(A=A,x=0.5),type = "response")
> y7<-predict(fit, data.frame(A=A,x=0.7),type = "response")
> y9<-predict(fit, data.frame(A=A,x=0.9),type = "response")
> df<-data.frame(A=A,y=y1,x=as.factor(0.1))
> df <- rbind(df,data.frame(A=A,y=y3,x=as.factor(0.3)))
> df <- rbind(df,data.frame(A=A,y=y5,x=as.factor(0.5)))
> df <- rbind(df,data.frame(A=A,y=y7,x=as.factor(0.7)))
> df <- rbind(df,data.frame(A=A,y=y9,x=as.factor(0.9)))
> pB<-ggplot(df,aes(x=A,y=y,color=x))+geom_line()+xlab("面積 A")+ylab("個体数 y")+ theme(legend.position=c(.1,.8))+ labs(colour="明るさ x")
# x = Brightness x

> x=seq(0,1,length=100)
> y4<-predict(fit, data.frame(A=4,x=x),type = "response")
> y8<-predict(fit, data.frame(A=8,x=x),type = "response")
> y12<-predict(fit, data.frame(A=12,x=x),type = "response")
> y16<-predict(fit, data.frame(A=16,x=x),type = "response")
> y20<-predict(fit, data.frame(A=20,x=x),type = "response")
> df<-data.frame(x=x,y=y4,A=as.factor(4))
> df <- rbind(df,data.frame(x=x,y=y8,A=as.factor(8)))
> df <- rbind(df,data.frame(x=x,y=y12,A=as.factor(12)))
> df <- rbind(df,data.frame(x=x,y=y16,A=as.factor(16)))
> df <- rbind(df,data.frame(x=x,y=y20,A=as.factor(20)))
> pC<-ggplot(df,aes(x=x,y=y,color=A))+geom_line()+xlab("明るさ x")+ylab("個体数 y")+ theme(legend.position=c(.1,.8))+ labs(colour="面積 A")

# Putting two charts together

> library(grid)
> grid.newpage()
> pushViewport(viewport(layout = grid.layout(1,2)))
> print(pB,vp=viewport(layout.pos.row=1,layout.pos.col = 1))
> print(pC,vp=viewport(layout.pos.row=1,layout.pos.col = 2))