ガンマ分布のglm()モデリング[R] - 井出草平の研究ノート

Kazuki Yoshidaさんによって作成されたものらしい。

farawayパッケージに含まれる半導体ウェハのデータを用いる。

library(faraway)
data(wafer)
plot(density(wafer$resist))

f:id:iDES:20210701224253p:plain

結果は連続的なものだが、右に傾いており、常に正の値を示している。このようなデータを分析する最も一般的な方法は、結果を対数変換することである。しかし、一般化線形モデルの枠組みの中で、より良い代替手段がある。

"resitivity of wafer in semiconductor experiment"

Description:
     A full factorial experiment with four two-level predictors.
Format:
     A data frame with 16 observations on the following 5 variables.
     x1 a factor with levels ‘-’ ‘+’
     x2 a factor with levels ‘-’ ‘+’
     x3 a factor with levels ‘-’ ‘+’
     x4 a factor with levels ‘-’ ‘+’
     resist Resistivity of the wafer
Source:
     Myers, R. and Montgomery D. (1997) A tutorial on generalized
     linear models, Journal of Quality Technology, 29, 274-291.

用語解説算術平均とは通常の平均である。 geometric mean（幾何平均）は、算術平均を対数スケールで計算したものである。

$\sum_{i = 1}^n log(Y_i)$

ここでは，オッズ比やレート比と同じように，2つの平均値を比較する比率を表すために「平均比」を使っている。

通常の方法 - 結果を対数変換した線形モデル（乗法幾何平均モデル）

res.lm.logY <- lm(log(resist) ~ x1 + x2 + x3 + x4, data = wafer)
summary(res.lm.logY)

Call:
lm(formula = log(resist) ~ x1 + x2 + x3 + x4, data = wafer)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.1757 -0.0622  0.0175  0.0877  0.1084 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.4405     0.0598   90.95  < 2e-16 ***
x1+           0.1228     0.0535    2.29  0.04243 *  
x2+          -0.2999     0.0535   -5.60  0.00016 ***
x3+           0.1784     0.0535    3.34  0.00665 ** 
x4+          -0.0561     0.0535   -1.05  0.31652    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.107 on 11 degrees of freedom
Multiple R-squared:  0.816, Adjusted R-squared:  0.75 
F-statistic: 12.2 on 4 and 11 DF,  p-value: 0.000491

exp(coef(res.lm.logY))

(Intercept)         x1+         x2+         x3+         x4+ 
   230.5536      1.1306      0.7409      1.1953      0.9454

結果が歪んでいて常に正である場合は、変換することでモデル化することができる。このモデルは次のようにモデル化している。

$\ E[log(Y_i)] = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4$

このように、対数スケールでの（算術的）平均変化、つまり、元のスケールでの幾何学的な「平均比」という観点からのみ解釈される。

例えば、x1＝＋とすると、対数結果の期待値（算術平均）は0.12277変化する。両辺の指数をとった共変量 $e^{0.12277} = 1.130624$ は、元の結果の幾何平均を掛け合わせた比率（元の尺度での幾何平均比）である。

切片を指数化したもの $\ e^{5.44} = 230.6$ は、すべての予測因子について-を持つウェハの幾何平均の期待値である。

したがって，このモデルは，予測変数による元の結果への乗法的効果を仮定している。

代替案1 - ガンマファミリとログリンクを用いた一般化線形モデル（乗算算術平均モデル）

res.glm.Gamma.log <- glm(formula = resist ~ x1 + x2 + x3 + x4,
                         family  = Gamma(link = "log"),
                         data    = wafer)
summary(res.glm.Gamma.log)

Call:
glm(formula = resist ~ x1 + x2 + x3 + x4, family = Gamma(link = "log"), 
    data = wafer)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.17548  -0.06486   0.01423   0.08399   0.10898  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.44552    0.05856  92.983  < 2e-16 ***
x1+          0.12115    0.05238   2.313 0.041090 *  
x2+         -0.30049    0.05238  -5.736 0.000131 ***
x3+          0.17979    0.05238   3.432 0.005601 ** 
x4+         -0.05757    0.05238  -1.099 0.295248    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.01097542)

    Null deviance: 0.69784  on 15  degrees of freedom
Residual deviance: 0.12418  on 11  degrees of freedom
AIC: 152.91

Number of Fisher Scoring iterations: 4

exp(coef(res.glm.Gamma.log))

(Intercept)         x1+         x2+         x3+         x4+ 
   231.7186      1.1288      0.7405      1.1970      0.9441

結果が歪んでいて常に正である場合は、ガンマ分布を用いてモデル化することができる。

リンク関数は、先の線形モデルとの整合性を考慮してlog()とし、以下のようにモデル化している。

$log(E[Y_i]) = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4$

となり、以下も成立する。

$E[Y_i] = e^{\beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4} \$

このように、元の尺度での算術平均 $E[Y_i]$ を推論することができる。

例えば、x1 = +とすると、対数算術平均結果が0.12115増加する。両辺の指数を取った共変量 $e^{0.12115} = 1.128794$ は、元の尺度での算術平均の結果に乗じる係数で、x1 = +の場合、他の変数のレベル内では、元の尺度での算術平均がx1 = -に比べて1.13倍になる。

また、指数化された切片 $e^{5.45} = 231.7$ は、すべての予測変数で-となったウェハの算術平均結果である。

したがって、このモデルは、予測変数による元の結果への乗算効果を仮定している。

相違点と類似点

これらの2つのモデルは非常によく似ており、今回のケースでも似たような係数を出している、指数化された係数は異なる解釈をしている（幾何学的「平均比」と算術的「平均比」）。

これは

$E[log(Y_i)|X] \ne log(E[Y_i|X])$

(対数スケールの平均は、元のスケールの平均の対数とは等しくない)

このように、線形回帰モデルに入る前に結果が対数変換されていれば、幾何平均に関する推論が可能となります。一方、一般化線形モデルのアプローチでは、元の尺度での算術平均についての推論が可能である。

いずれの場合も，各予測変数の効果は乗法的である（平均の変化率）。

代替案2 ガンマ・ファミリと同一性リンクを持つ一般化線形モデル（加法的算術平均モデル）

res.glm.Gamma.identity <- glm(formula = resist ~ x1 + x2 + x3 + x4,
                         family  = Gamma(link = "identity"),
                         data    = wafer)
summary(res.glm.Gamma.identity)

Call:
glm(formula = resist ~ x1 + x2 + x3 + x4, family = Gamma(link = "identity"), 
    data = wafer)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.190394  -0.066533   0.005538   0.091549   0.126838  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   235.43      14.21  16.572 3.97e-09 ***
x1+            27.83      12.24   2.274 0.044017 *  
x2+           -65.74      12.68  -5.184 0.000302 ***
x3+            36.94      12.32   2.999 0.012102 *  
x4+           -11.88      12.15  -0.978 0.349294    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.0124071)

    Null deviance: 0.69784  on 15  degrees of freedom
Residual deviance: 0.14009  on 11  degrees of freedom
AIC: 154.84

Number of Fisher Scoring iterations: 6

各予測変数の効果が元の尺度上で相加的であると考えられる場合，同一性リンク関数を持つ一般化線形モデルを使用できる。
この場合、生の係数は元の尺度になる。x1 = +を持つことは、算術平均の結果に27.83を加える（加法的効果）。
生の切片（235.43）は、すべての予測変数について-を有するウェーハの算術平均結果である。