はじめに

　本講義では、先週学んだ重回帰分析について発展的な話題について説明する。

独立変数が非線形な線形回帰

　線形回帰の「線形」とはパラメータについてのことであり、予測変数は「線形」でなくてかまわない。例えば、 \[ \frac 1 {y_i} = \beta_0 + \frac {\beta_1} {x_i} + \epsilon_i, \:\: i = 1, ..., n \] や \[ y_i = \epsilon_i \times \exp(\beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i}), \:\: i = 1,..., n \] も、線形回帰である。後述するように、\(x\)または\(y\)あるいは両方に対する変換を行って、線形の関係を実現することもできる。ただし、多くの場合、誤差の正規性や等分散性について何らかの妥協をしなければならなくなることを心にとどめておく必要がある。

　なお、誤差が変数変換の影響を受けない例の一つは、\(x\)と\(y\)の多項式の関係である。 \[ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \cdots + \beta_k x_i^k + \epsilon_i,\:\: i = 1,...,n \]

　次は、2次の多項式の回帰モデルの例である。

例　コムギにおける収量とタンパク質含量の関係

　以下は、コムギの収量(\(X\))と穀粒のタンパク質含量(\(Y\))のデータ（Snedecor and Cockran, 1989より引用）である。

x <- c(5, 8, 10, 11, 14, 16, 17, 17, 18, 20, 22, 24, 26, 30, 32, 34, 36, 38, 43)
y <- c(16.2, 14.2, 14.6, 18.3, 13.2, 13.0, 13.0, 13.4, 10.6, 12.8, 12.6, 11.6, 11.0, 9.8, 10.4, 10.9, 12.2, 9.8, 10.7)
df <- data.frame(y, x)

　まず、両者の関係を散布図を描いて確認してる。

plot(y ~ x, data = df)

　収量\(x\)が増えると、タンパク質含量\(y\)が減る傾向がある。 　 そこで、タンパク質含量\(y\)を収量\(x\)の上に回帰してみる。

model1 <- lm(y ~ x, df)
plot(y ~ x, data = df)
abline(model1)

anova(model1)
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## x          1 53.760  53.760  26.655 7.804e-05 ***
## Residuals 17 34.286   2.017                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

割とよくあてはまっているように見えますが、収量\(x\)の低いところと、高いところで\(y\)が過小に推定されていることが分かる。

　そこで、2次の項までの多項式回帰を行ってみる。y ~ x + I(x^2)というモデルの記述の仕方に注意する。

model2 <- lm(y ~ x + I(x^2), data = df) 
summary(model2)
## 
## Call:
## lm(formula = y ~ x + I(x^2), data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1170 -0.4762 -0.1914  0.3699  3.7172 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.676720   1.372098  13.612 3.25e-10 ***
## x           -0.436731   0.129363  -3.376  0.00385 ** 
## I(x^2)       0.005869   0.002665   2.202  0.04267 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.282 on 16 degrees of freedom
## Multiple R-squared:  0.7012, Adjusted R-squared:  0.6638 
## F-statistic: 18.77 on 2 and 16 DF,  p-value: 6.361e-05

2次の項も有意である。

　多項式回帰モデルは、次のようなコマンドでもあてはめることができる。y ~ poly(x, 2)というモデルの記述の仕方に気をつける。

model2.2 <- lm(y ~ poly(x, 2), data = df) # どちらの書き方でも構いません
summary(model2.2)
## 
## Call:
## lm(formula = y ~ poly(x, 2), data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1170 -0.4762 -0.1914  0.3699  3.7172 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  12.5421     0.2942  42.632  < 2e-16 ***
## poly(x, 2)1  -7.3321     1.2824  -5.718 3.18e-05 ***
## poly(x, 2)2   2.8239     1.2824   2.202   0.0427 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.282 on 16 degrees of freedom
## Multiple R-squared:  0.7012, Adjusted R-squared:  0.6638 
## F-statistic: 18.77 on 2 and 16 DF,  p-value: 6.361e-05

　2次の多項式回帰モデルと単回帰モデルを比較してみる。

anova(model1, model2)
## Analysis of Variance Table
## 
## Model 1: y ~ x
## Model 2: y ~ x + I(x^2)
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     17 34.286                              
## 2     16 26.312  1    7.9745 4.8492 0.04267 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

この場合は、2次の項の回帰係数の\(t\)検定の結果と一致する。

　最後に、推定値\(\hat y\)と予測値\(\tilde y\)の信頼区間を求めて図示してみる。

plot(y ~ x, data = df, ylim = c(7, 20))
xp <- seq(min(x), max(x), 0.5)
y.est <- predict(model2, data.frame(x = xp), interval = "confidence")
matlines(xp, y.est, lty = c(1, 2, 2), col = c("black", "orange", "orange"))
y.pred <- predict(model2, data.frame(x = xp), interval = "prediction")
matlines(xp, y.pred, lty = c(1, 2, 2), col = c("black", "green", "green"))

例　木の胸高直径と樹高から体積を推定する

伐採された31本のブラックチェリーの木の胸高直径\(D\)、樹高\(H\)、体積\(V\)のデータをもとに、胸高直径と樹高から体積を予測するモデルを作る。これら変数は、木の主幹の円柱状の形から \[ V = c H D^2 \] というモデルに従うのではないかと考えられる。なお、\(c\)は変換のための定数である。したがって、これら変数の対数をとると、 \[ \log V = \beta_0 + \beta_1 (\log H) + \beta_2 (\log D) \] というような関係がみられると考えられ、さらに\(\beta_1\)や\(\beta_2\)は、それぞれ、1, 2に近い値をとると期待される。

　このことを、Rのデータtreesを使って確認してみる。treesはirisと同様に、いつでも呼び出して使うことができる。まず、データの準備する。

head(trees)
##   Girth Height Volume
## 1   8.3     70   10.3
## 2   8.6     65   10.3
## 3   8.8     63   10.2
## 4  10.5     72   16.4
## 5  10.7     81   18.8
## 6  10.8     83   19.7
log.trees <- log(trees)
colnames(log.trees) <- c("log.D", "log.H", "log.V") # 分かりやすいように名前を変えた
head(log.trees)
##      log.D    log.H    log.V
## 1 2.116256 4.248495 2.332144
## 2 2.151762 4.174387 2.332144
## 3 2.174752 4.143135 2.322388
## 4 2.351375 4.276666 2.797281
## 5 2.370244 4.394449 2.933857
## 6 2.379546 4.418841 2.980619

　回帰モデルを当てはめる。

model <- lm(log.V ~ log.H + log.D, data = log.trees)
summary(model)
## 
## Call:
## lm(formula = log.V ~ log.H + log.D, data = log.trees)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.168561 -0.048488  0.002431  0.063637  0.129223 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.63162    0.79979  -8.292 5.06e-09 ***
## log.H        1.11712    0.20444   5.464 7.81e-06 ***
## log.D        1.98265    0.07501  26.432  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08139 on 28 degrees of freedom
## Multiple R-squared:  0.9777, Adjusted R-squared:  0.9761 
## F-statistic: 613.2 on 2 and 28 DF,  p-value: < 2.2e-16

　期待された通り、\(b_1, b_2\)は、1, 2に近い値となる。

回帰係数\(\beta_i\)の95%信頼区間を示す。

confint(model)
##                 2.5 %    97.5 %
## (Intercept) -8.269912 -4.993322
## log.H        0.698353  1.535894
## log.D        1.828998  2.136302

\(\beta_1=1\)と\(\beta_0=2\)は、信頼区間内に入っている。

　最後に、胸高直径\(D = 16\)、樹高\(H = 80\)のときの木の体積\(V\)の推定値とその95%信頼区間を計算してみる。データの変換をしなければならないことに注意する。

df <- data.frame(log.D = log(16), log.H = log(80))
(log.V.pred <- predict(model, df, interval = "confidence"))
##       fit      lwr      upr
## 1 3.76072 3.719341 3.802099
(V.pred <- exp(log.V.pred))
##        fit     lwr      upr
## 1 42.97935 41.2372 44.79511
hist(trees$Volume)
abline(v = V.pred, col = c("black", "orange", "orange"))

信頼区間の範囲も狭く、信頼度の高い推定ができていることが分かる。胸高直径、樹高は伐採前にも計測ができるため、こうした回帰モデルは有用と考えられる。

共分散分析

　共分散分析（analysis of covariance: ANCOVA）は、2つのタイプの予測変数、回帰分析のように量的な変量と分散分析のようにカテゴリカルな変量、をもつ線形モデルである。これは、より一般的な形で定式化することもできるが、ここでは、それぞれのタイプの予測変数が一つずつある場合を例にして説明していく。

　回帰と分散分析を融合する理論的根拠は、様々な実験においては、問題を回帰としてモデル化するだけでも、分散分析としてモデル化するだけでも、不十分な場合があるためである。また、分散分析に量的な共変量を導入するか、あるいは、回帰に処理やグループの効果を導入することにより、データ内の変動をよく説明でき、よりよいモデリングや予測が可能になる場合もある。

　以下の例では、2つの降圧剤の効果を、薬剤摂取後の血圧低下を分析することによって比較している。ただし、例えば、初期圧力が180の場合と90の場合に、50低下することの意味は同じではないと考えられるため、降圧剤を摂取する前の初期血圧を考慮に入れたモデルを作るべきだと考えた。

　すなわち、モデルを、 \[ y_{ij} = \mu + \alpha_i + \beta(x_{ij} - \bar x) + \epsilon_{ij}, \:\: i = 1,...,a; j = 1,...,n \] と仮定した。ここで、\(a\)は処理・水準の数であり、\(n\)はそれぞれの水準での標本数である。また、\(\bar x = \frac 1 {an} \sum_{i,j} x_{ij}\)であり、\(x\)の全体平均である。なお、これまでのモデルと同様に、誤差\(\epsilon_i\)は、独立に平均0、分散\(\sigma^2\)の正規分布に従うと仮定する。回帰式をより簡単にするために、共変量\(x_{ij}\)について、\(x_{ij} - \bar x\)として平均を0に基準化してある。 　 なお、ここで、\(\bar x_i = \frac 1 n \sum_j x_{ij}\)を\(i\)番目の処理における\(x_{ij}\)の平均とする。また、平均\(\bar y\)、\(\bar y_i\)についても、\(\bar x\), \(\bar x_i\)と同様に定義する。

　共分散分析では、以下のように平方和と積和を計算することができる。 \[ S_{xx} = \sum_{i=1}^a \sum_{j=1}^n (x_{ij} - \bar x)^2 \\ S_{xy} = \sum_{i=1}^a \sum_{j=1}^n (x_{ij} - \bar x)(y_{ij} - \bar y) \\ S_{yy} = \sum_{i=1}^a \sum_{j=1}^n (y_{ij} - \bar y)^2 \\ T_{xx} = \sum_{i=1}^a (\bar x_i - \bar x)^2 \\ T_{xy} = \sum_{i=1}^a (\bar x_i - \bar x)(\bar y_i - \bar y) \\ T_{yy} = \sum_{i=1}^a (\bar y_i - \bar y)^2 \\ Q_{xx} = \sum_{i=1}^a \sum_{j=1}^n (x_{ij} - \bar x_i)^2 = S_{xx} - T_{xx} \\ Q_{xy} = \sum_{i=1}^a \sum_{j=1}^n (x_{ij} - \bar x_i)(y_{ij} - \bar y_i) = S_{xy} - T_{xy} \\ Q_{yy} = \sum_{i=1}^a \sum_{j=1}^n (y_{ij} - \bar y_i)^2 = S_{yy} - T_{yy} \\ \]

　これら平方和と積和を使うと、上述したモデルのパラメータの推定値を計算できる。 \[ \hat \mu = \bar y \\ b = \hat \beta = Q_{xy} / Q_{xx} \\ \hat \alpha_i = \bar y_i - \bar y - b(\bar x_i - \bar x) \]

　また、分散の推定量は、 \[ {\hat \sigma}^2 = s^2 = MSE = \frac {SSE} {a(n-1)-1} \] として計算できる。ここで、\(SSE = Q_{yy} - Q_{xy}^2 / Q_{xx}\)である。

　もし、処理の効果がなければ、すなわち、全ての\(\alpha_i = 0\)となり、上述したモデルは通常の単回帰となり、 \[ y_{ij} = \mu + \beta(x_{ij} - \bar x) + \epsilon_{ij}, \:\: i = 1,...,a; j = 1,...,n \] と表せる。

　この場合、誤差の平方和は、\(SSE' = S_{yy} - S_{xy}^2 / S_{xx}\)で、\(an-2\)の自由度をもつ。 　 　したがって、\(H_0: \alpha_i = 0\)の検定は、\(F\)統計量 \[ F = \frac {(SSE' - SSE) / (a-1)}{SSE / (a(n-1)-1)} \] をもとに行えばよいことがわかる。この統計量は、帰無仮説のもとで自由度\(a-1\)と\(a(n-1)-1\)をもつ\(F\)分布に従う。

例　降圧剤の効果と摂取前の初期血圧の関係

　Kodlin（1951）は、血圧を下げるための2つの降圧剤（AとB）を比較した実験を報告した（Vidakovic B (2011）より引用）。2つの動物群を2つの物質に無作為に割り付け、血圧の減少を記録した。前述のように初期血圧を記録した。血圧実験のデータは以下の表にある。初期圧力が圧力の低下に与える影響を考慮して、2つの物質の比較を行う。 　共分散分析の結果を検討し、それらを初期圧力の潜在的な影響を無視して得られた結果と比較する。検定には有意水準\(\alpha = 0.05\)を用いる。

まずデータの準備をする。

# data preparation
x <- c(135, 125, 125, 130, 105, 130, 140, 93, 110, 100, 90, 135, 130, 115, 110, 140, 130, 95, 90, 105)
y <- c(45, 45, 20, 50, 25, 37, 50, 20, 25, 15, 34, 55, 50, 45, 30, 45, 45, 23, 40, 35)
g <- as.factor(c(rep("A", 10), rep("B", 10)))
a <- 2
n <- 10

　次に、平方和と積和の計算をする。

# preparation for calculation
x1 <- x[g == "A"]; x2 <- x[g == "B"]
y1 <- y[g == "A"]; y2 <- y[g == "B"]
# calculate SS
x1b <- mean(x1)
x2b <- mean(x2)
y1b <- mean(y1)
y2b <- mean(y2)
xb <- mean(x)
yb <- mean(y)
Sxx <- sum((x - xb)^2)
Sxy <- sum((x - xb) * (y - yb))
Syy <- sum((y - yb)^2)

Qxx <- sum((x1 - x1b)^2 + (x2 - x2b)^2)
Qxy <- sum((x1 - x1b) * (y1 - y1b) + (x2 - x2b) * (y2 - y2b))
Qyy <- sum((y1 - y1b)^2 + (y2 - y2b)^2)

　次に、母数の推定値を計算する。

# estimators of model parameters mu, alpha_i, beta
(mu <- yb)
## [1] 36.7
(b <- Qxy / Qxx)
## [1] 0.5174961
(alpha <- c(y1b, y2b) - c(yb, yb) - b * (c(x1b, x2b) - c(xb, xb)))
## [1] -4.871365  4.871365

分散の推定値を計算する。

SSE <- Qyy - Qxy^2 / Qxx
(MSE <- SSE / (a * (n - 1) - 1))
## [1] 60.65422

　処理の効果（降圧剤の違いによる効果）を検定する。

SSEp <- Syy - Sxy^2 / Sxx
# F-test for tesing H_0: alpha_i = 0 (all i)
F1 <- ((SSEp - SSE) / (a - 1)) / MSE
F1
## [1] 7.632095
pvalF1 <- 1 - pf(F1, a - 1, a * (n - 1) - 1)
pvalF1
## [1] 0.01331085

　検定の結果、帰無仮説が有意であることが分かる。したがって、降圧剤AとBの効果は、\(\alpha = 0.0133\)で有意に異なっていると結論付けられる。

　共変量として加えた初期血圧の効果（傾き）の検定も行うことができる。

# F-tet for testing H_0: beta = 0
F2 <- (Qxy^2 / Qxx) / MSE
F2
## [1] 24.56677
pvalF2 <- 1 - pf(F2, 1, a * (n - 1) - 1)
pvalF2
## [1] 0.0001200347

こちらも高度に有意である。

　最後に、降圧剤AとBを分けて、初期値\(x\)に対する血圧の減少\(y\)について散布図を描いてみる。散布図の2つの回帰式\(\hat y = 36.7 \pm 4.8714 + 0.5175\cdot (x - 116.65)\)も重ね書きする。

plot(x, y, col = as.numeric(g))
abline(mu + alpha[1] - xb * b, b, col = "black")
abline(mu + alpha[2] - xb * b, b, col = "red")

　なお、共変量とした血圧の初期値を無視して分散分析を行うと、降圧剤AとB間の差は有意ではない（\(p\)値は0.2025）。

# ANOVA
model <- lm(y ~ g)
anova(model)
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value Pr(>F)
## g          1  245.0  245.00  1.7492 0.2025
## Residuals 18 2521.2  140.07

また、同じように、2つの処理を仮に1つにして応答\(y\)を\(x\)に回帰した場合は、回帰は有意だが、\(F\)検定の\(p\)値は共分散分析のときに比べて8倍以上大きなものとなる。

# regression
model <- lm(y ~ x)
summary(model)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.643  -5.650   0.815   6.996  15.885 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -18.3865    14.2175  -1.293  0.21228   
## x             0.4722     0.1206   3.915  0.00102 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.111 on 18 degrees of freedom
## Multiple R-squared:  0.4599, Adjusted R-squared:  0.4299 
## F-statistic: 15.33 on 1 and 18 DF,  p-value: 0.001015

これらの結果からも、共分散分析の有効性が示唆される。

　なお、上の一連の解析は、以下のようにするとRの関数lmを用いて実行できる。

xs <- x - xb
model <- lm(y ~ g + xs)
summary(model)
## 
## Call:
## lm(formula = y ~ g + xs)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.1497  -4.4518  -0.0662   4.7338  12.2199 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  31.8286     2.4783  12.843 3.54e-10 ***
## gB            9.7427     3.5266   2.763  0.01331 *  
## xs            0.5175     0.1044   4.956  0.00012 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.788 on 17 degrees of freedom
## Multiple R-squared:  0.6272, Adjusted R-squared:  0.5834 
## F-statistic:  14.3 on 2 and 17 DF,  p-value: 0.0002276

　 　 ただし、回帰モデルの分散分析には、通常のanova関数ではなく、パッケージcarに含まれるAnova関数を用いて行う必要がある。これは、Type IIの平方和とよばれる平方和を計算する必要なためである。

# 分散分析には、carパッケージのAnova関数をつかう。
require(car)
## Loading required package: car
## Loading required package: carData
Anova(model, type="II")
## Anova Table (Type II tests)
## 
## Response: y
##            Sum Sq Df F value  Pr(>F)    
## g          462.92  1  7.6321 0.01331 *  
## xs        1490.08  1 24.5668 0.00012 ***
## Residuals 1031.12 17                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

比較のために、通常のanova関数を用いた分散分析を行う。通常のanova関数では、Type Iの平方和が計算される。

anova(model)
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## g          1  245.0  245.00  4.0393 0.06060 .  
## xs         1 1490.1 1490.08 24.5668 0.00012 ***
## Residuals 17 1031.1   60.65                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model2 <- lm(y ~ xs + g)
anova(model2)
## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## xs         1 1272.16 1272.16 20.9740 0.0002663 ***
## g          1  462.92  462.92  7.6321 0.0133108 *  
## Residuals 17 1031.12   60.65                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type Iの平方和に基づく分散分析では、変数としてモデルに含める順番に依存して結果が変わってしまうことに注意する。

参考文献

以下は、本テキストを作成するにあたり参考にした書籍である。