Murakami's Memorandum

Home Page > My Memo               災害科学ホワイトボード   災害科学コース   理学部

GLM入門編1

(工事中)

一般化線形モデル

 一般化線形モデル(Generalized Linear Model)のお勉強中のメモ書きです。決して鵜呑みにしないでください。

 モチベーションとしては,カウントデータ(日数)にExcelで直線を当てはめて議論する人がちらほらいるのだが,それでよいのだろうかという漠然とした疑問をもったのですっきりしたいということです。例えば,真夏日の経年変化(横軸:年−縦軸:真夏日)に対して直線を当てはめ1年当たり何日増えているので温暖化の傾向が見えるというような話なのだが,縦軸のデータはそもそも単位は日なので連続データというわけではない。また,データとしてはそもそもゼロ以下の値をとらない。しかし,回帰直線をモデルという観点(そもそもそう見ていないのか)から見るとセロを越えてマイナスになる領域が出てくる。卒論などでこの手が話が出てくる時には,専門の指導教員の目を経ているのであろうから問題がないのかもしれないしのだが...。

 データに正規分布を仮定しないで解析をするために,データの確率分布を見極める必要がある。
確率分布 特徴
ポアソン分布 データが離散値,ゼロ以上の範囲,上限とくになし,標本平均=標本分散
負の二項分布 データが離散値,ゼロ以上の範囲,上限とくになし,標本平均<標本分散
二項分布 データが離散値,ゼロ以上で有限の範囲(0 1 2 ... N)
正規分布 データが連続値,範囲(-∞〜∞)
対数正規分布 データが連続値,範囲(0〜∞)
ガンマ分布 データが連続値,範囲(0〜∞)
ベータ分布 データが連続値,範囲(0〜1)

ref.) http://hosho.ees.hokudai.ac.jp/~kubo/stat/2008/b/kubostat2008b.pdf

データ特性

 真夏日(hday)の経年変化にモデルを当てはめる演習問題をRのGLMを使いながら試みる。最初に,matatu.datというデータファイルを読み込み,ヒストグラム,正規確率プロットをしてみる。

> D <- read.table("C:/Users/murakami/Documents/manatu.dat", header=TRUE, 
+   sep="", na.strings="NA", dec=".", strip.white=TRUE)
> Hist(D$hday, scale="frequency", breaks="Sturges", col="darkgray")
> qq.plot(D$hday, dist= "norm", labels=FALSE)

GLM ポアソン回帰

 真夏日データは(見た目では)正規分布しないことがわかる。ではポアソン分布かというと,平均=分散は成り立っていない。

> mean(D$hday)
[1] 19.08929
> var(D$hday)
[1] 108.8464

 いきなり困ってしまうがとりあえずデータはポアソン分布にしたがうものとして解析を実行してみる。Rの一般化線形モデルglmで,family=poissonを選択すると標準では対数変換を指定することになる。対数変換するとデータはどうなるかを見てみる。

  • 生のデータ
  • 対数変換したデータ

 生データではバラツキが年毎に大きくなる傾向があるが,対数変換した結果バラツキの大きさがほぼ一様になった(?)。一般化線形モデルで,データのポアソン分布を指定すると計算途中で次の線形の式を当てはめることになる。

 ポアソン回帰  log(μ) = α+βX  <---> 期待値 μ=exp(α+βX)

> GLM.1 <- glm(hday ~ year, family=poisson(log), data=D)
> summary(GLM.1)

Call:
glm(formula = hday ~ year, family = poisson(log), data = D)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-4.50925  -1.28531  -0.03577   1.22616   4.03139  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -26.918753   3.821829  -7.043 1.88e-12 ***
year          0.015066   0.001926   7.823 5.16e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 308.59  on 55  degrees of freedom
Residual deviance: 246.31  on 54  degrees of freedom
AIC: 509.91

Number of Fisher Scoring iterations: 4

 ポアソン回帰  log(μ) = α+βX  <---> 期待値 μ=exp(α+βX)

 ここで得られたモデルは, hday = exp( -26.9 + 0.015*year) ということになる。このモデルならhdayが負になることはない。しかし,最初からポアソン分布から期待される分散よりも大きいことがわかっていたが, 「Residual deviance/degree of freedom=246.1/54 >1」となり,過分散(overdispersion)なので最適なモデルとは言えない。

GLM 負の二項分布

 ポアソン分布と仮定して解析した結果が過分散の場合には,ものの本によれば「負の二項分布」として解析すればよいなどと紹介されている。負の二項分布とは何かはとりあえず置いておいて,データが負の二項分布に従うものとして解析してみる。

>library(MASS)
> NGLM <- glm.nb(hday ~ year,data=D)
> summary(NGLM)

Call:
glm.nb(formula = hday ~ year, data = D, init.theta = 5.25104821567019, 
    link = log)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.43375  -0.66780  -0.02390   0.57497   1.57584  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -25.470269   8.131666  -3.132 0.001735 ** 
year          0.014335   0.004104   3.493 0.000478 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for Negative Binomial(5.251) family taken to be 1)

    Null deviance: 71.164  on 55  degrees of freedom
Residual deviance: 58.288  on 54  degrees of freedom
AIC: 406.03

Number of Fisher Scoring iterations: 1


              Theta:  5.25 
          Std. Err.:  1.29 

 2 x log-likelihood:  -400.026 

負の二項分布 log(ρ/1+ρ) = log(μ/μ+k) = α+βX ;期待値 μ=ρθ

 今度のモデルでは,過分散の程度はかなり改善されている。切片,傾斜の値は少し変わった。

モデルの比較

 いわゆる直線回帰,一般化線形モデル(ポアソン回帰),一般化線形モデル(負の二項分布)の比較をしてみる。

 その前に従来どおりの回帰直線をglm関数で求める。

> GLM.2 <- glm(hday ~ year, family=gaussian(identity), data=D)
> summary(GLM.2)

Call:
glm(formula = hday ~ year, family = gaussian(identity), data = D)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-17.9060   -5.3898   -0.2948    4.8927   22.6522  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -543.85735  154.45265  -3.521 0.000882 ***
year           0.28424    0.07798   3.645 0.000602 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for gaussian family taken to be 88.97263)

    Null deviance: 5986.6  on 55  degrees of freedom
Residual deviance: 4804.5  on 54  degrees of freedom
AIC: 414.23

Number of Fisher Scoring iterations: 2

回帰直線  期待値 μ = α+βX 

 では,データに各モデルの予測式(予測値)を重ねてみると

>plot(D$year,D$hday)                               # 散布図
>lines(D$year,GLM.1$fitted.values,col=3)   # 一般化線形モデル(ポアソン回帰)
>lines(D$year,NGLM$fitted.values,col=2)    # 一般化線形モデル(負の二項分布)
>lines(D$year,GLM.2$fitted.values,col=1)   # 回帰直線[一般化線形モデル(ガウス分布)]
>legend("topleft",pch=c("-","-","-"),col=c(3,2,1),c("GLM(Poisson)","GLM(NegativeBinominal)","GLM(Gauss)"))

 データのある範囲内では3つの予測値の差はほとんどないのでどれでもいいじゃんと言われそうではあるが,元来のデータの持つ性質を考えてモデルを構築するという基本的な立場をつらぬけば,「負の二項分布」を仮定したモデルが最適モデルとなる(AIC最小)。

ポアソン分布


負の二項分布


参考文献・URL

compGLM.jpg hist_plot.jpeg qq_plot.jpeg y-h.jpeg y-logh.jpeg

最終更新時間:2009年09月03日 09時19分25秒


HomePage > FrontPage