1 Review of multiple linear regressions

本章では、一般化加法モデル(GAM)の説明に入る前に、多くの人に馴染みのある重回帰分析(multiple linear regression)について説明する、なぜなら、GAMは重回帰分析を拡張したものだからである。

1.1 Light levels and size of the human visual system

Pearce and Dunbar (2012) は、ヒトの集団が住んでいる標高と眼窩の容量に正の関連があることから、住んでいる環境の光量がヒトの視覚システムの進化の原動力になっていると結論付けた。本章ではこの論文のデータを用いて重回帰分析について説明を行う。重回帰分析でデータを探索し、モデルを構築し、モデルを当てはめ、モデル選択を行い、モデルの妥当性を確認する方法はGAMでもほとんど同じように適用できる。

1.2 The variables

Pearce and Dunbar (2012) は、オックスフォード大学博物館にある55人の成人の頭蓋骨から、頭蓋の容量(cranial capacity: CC)と眼窩容量(orbital volume)、大後頭孔(foramen magnum: FM)などの測定を行った。

データは以下のとおりである。平均眼窩容量(mean orbital volume)が目的変数であり、それ以外の変数は説明変数である1。AbsoluteLatitudeとMinimum_Tempreture_celsiusはそれぞれ 頭蓋が発見された場所の標高と最低気温、FMarea_intercondyleは体格の大きさを示す指標、Minimum_Illuminanceはlogスケールで表した光の強度、Genderは頭蓋の性別である。

hvs <- read_delim("data/HVS.txt")

datatable(hvs,
          filter = "top",
          options = list(scrollX = 20))


変数名が長いので、以下のように短く変更する。

hvs %>% 
  rename(OrbitalV = MeanOrbitalVolume,
         LatAbs = AbsoluteLatitude,
         CC = CranialCapacity,
         Illuminance = Minimum_Illuminance,
         Temperature = Minimum_Temperature_celsius,
         FM = FMarea_intercondyle) %>% 
  mutate(fPopulation = factor(Population),
         fGender = factor(Gender)) -> hvs

新しい列名は以下の通り。

colnames(hvs)
##  [1] "Museum"      "Findsite"    "Gender"      "Population"  "Latitude"   
##  [6] "LatAbs"      "CC"          "FM"          "OrbitalV"    "Illuminance"
## [11] "Temperature" "fPopulation" "fGender"

1.3 Protocol for the analysis

いかなるデータ分析も、以下の手順に沿って行わなければならない。

  1. Data exploration
    外れ値がないか、多重共線性(説明変数同士の強い相関)がないか、目的変数と説明変数の関係がどうか、ゼロ過剰はないか、サンプリングの収集が時間や場所によってばらついていないか、などをチェックする必要がある。

  2. Model application
    1の作業で分かったことや研究仮説をもとに、適切なモデルを適用する。今回は重回帰分析を行うが、様々なモデルを適用可能である。

  3. Check the result
    モデルを当てはめたら、どの変数が有意な影響を持つかを調べ、そうでなかった変数についてはどうするかを考える。

  4. Model validation
    最後に、作成したモデルが前提を満たしているかをチェックする。満たしていれば結果の解釈や結果の作図を行い、満たしていなければモデルを改善する必要がある(GAM、GLM、GLMを使うなど)。

1.4 Data exploration

まずは手順1のデータ探索を行う。データ探索については Zuur et al. (2010) に詳しい。

1.4.1 欠損値の確認

まず、欠損値がないかを調べる。

colSums(is.na(hvs))
##      Museum    Findsite      Gender  Population    Latitude      LatAbs 
##           0           0           0           0           0           0 
##          CC          FM    OrbitalV Illuminance Temperature fPopulation 
##           0           1           0           0           0           0 
##     fGender 
##           0

FMの列に1つ欠損値があるので取り除く。

hvs_b <- na.omit(hvs)

1.4.2 外れ値の確認

続いて、外れ値がないかを確認する。ここでは、連続値の説明変数(Latitude、CC、FM、Illuminance、Temperature)についてCleveland Dotplotを作成する。Dotplotは縦軸にサンプル番号、横軸に実際の値をとる。

図は以下のとおりである(図1.1)。図を見る限り、外れ値はなさそうだ。

hvs_b %>% 
  select(LatAbs, CC, FM, Illuminance, Temperature) %>% 
  mutate(sample_number = 1:n()) %>% 
  pivot_longer(cols = 1:5,
               names_to = "var",
               values_to = "values") %>% 
  ggplot(aes(x = values, y = sample_number))+
  geom_point(alpha = 1)+
  facet_rep_wrap(~var,
                 scales = "free_x")+
  labs(x = "Values of the variable", y = "Sample number")+
  theme(aspect.ratio = 1)
Cleveland dotplot for covariates.

図1.1: Cleveland dotplot for covariates.

1.4.3 多重共線性の確認

次に、多重共線性(説明変数同士の強い相関)がないかを調べる。もしあると、推定結果にバイアスが生じてしまう。

説明変数同士の関連を調べたところ(図1.2)、LatAbsとIlluminance、IlluminanceとTemperature、TemperatureとLatAbsに強い相関があることが分かる。また、CCは男性で高い傾向があることが分かったので、CCとfGenderを同じモデルに説明変数として入れない方がよさそうである(今回はfGenderを用いない)。

ggpairs(hvs_b %>% select(LatAbs, CC, FM, Illuminance, Temperature, fGender))
correlation between covariates

図1.2: correlation between covariates

1.4.4 目的変数と説明変数の関係の確認

最後に、目的変数と説明変数の関連を調べる(図1.3)。図には局所回帰(LOESS)による回帰曲線を追加している。

図から、LatAbsとOrbitalVの間に線形の関係がありそうだということが分かる(’CC’も?)。IlluminanceやTemperatureにも同様のことがいえるが、これは変数間に強い相関があることを考えれば当然だろう。多重共線性を考慮し、Pearce and Dunbar (2012) にもとづいて解析ではLatAbsを用いてIlluminanceとTemperatureは用いないこととする。

hvs_b %>% 
  select(LatAbs, CC, FM, Illuminance, Temperature, OrbitalV) %>% 
  pivot_longer(cols = 1:5,
               names_to = "var",
               values_to = "values") %>% 
  ggplot(aes(x = values, y = OrbitalV))+
  geom_point(alpha = 1)+
  facet_rep_wrap(~var,
                 scales = "free_x")+
  labs(x = "Explanatory Variables", y = "OrbitalV")+
  theme(aspect.ratio = 1)+
  geom_smooth(method = "loess", se= F, color = "grey32",
              span = 0.9)
Relationship between explanatory variables and orbital volume (OrbitalV). A LOESS smoother was added to the plot.

図1.3: Relationship between explanatory variables and orbital volume (OrbitalV). A LOESS smoother was added to the plot.

1.5 Multiple linear regression

1.5.1 Underlying statistical theory

それでは、重回帰分析を行う。まずは説明のために説明変数が1つだけのモデル(= 単回帰)を考えよう。

標高のみを説明変数とする単回帰モデルは以下のように実行できる。

M1 <- lm(OrbitalV ~ LatAbs, data = hvs_b)

summary(M1)
## 
## Call:
## lm(formula = OrbitalV ~ LatAbs, data = hvs_b)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.9893 -1.4166 -0.1616  1.5037  3.8887 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 22.91064    0.51412  44.563  < 2e-16 ***
## LatAbs       0.06598    0.01411   4.677 2.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.253 on 52 degrees of freedom
## Multiple R-squared:  0.2961, Adjusted R-squared:  0.2825 
## F-statistic: 21.87 on 1 and 52 DF,  p-value: 2.111e-05

これは、実際には何をやっているのだろうか?
lm関数で短回帰を実行するとき、私たちは下記のモデルを実行している。なお、\(i\)はサンプル番号(今回は55個の頭蓋がある)を、2行目は\(\epsilon_i\)が平均0、分散\(\sigma^2\)の正規分布に従うことを表す。

\[ \begin{aligned} OrvitalV_i &= \alpha + \beta\times LatAbs_i + \epsilon_i \;\; (i = 1,2,\dots, 55)\\ \epsilon_i &\sim N(0,\sigma^2) \end{aligned} \]

このモデルで推定するパラメータは\(\alpha\)、\(\beta\)、\(\sigma\)である。\(\beta\)は標高(LatAbs)が眼窩容量(OrbitalV)に与える影響の大きさを表し、式からわかるように標高が1増えると眼窩容量が\(\beta\)増えることを表す。先ほどlm関数で実行した結果でEstimateの下に書かれていた数字(\(22.91...と0.065...\))は\(\alphaと\beta\)の推定値を示しているのである。

もし以下のように行列を定義すると、
\[ \begin{aligned} \mathbf{y} = \begin{pmatrix} OrbitalV_1\\ OrbitalV_2\\ \vdots\\ OrvitalV_{55} \end{pmatrix} , \mathbf{X} = \begin{pmatrix} 1 & LatAbs_1\\ 1 & LatAbs_2\\ \vdots & \vdots\\ 1 & LatAbs_{55} \end{pmatrix}, \mathbf{\beta} = \begin{pmatrix} \alpha\\ \beta \end{pmatrix}, \mathbf{\epsilon} = \begin{pmatrix} \epsilon_1\\ \epsilon_2\\ \vdots \\ \epsilon_{55} \end{pmatrix} \end{aligned} \]

モデル式は以下のように書ける。
\[ \mathbf{y = X\times\beta + \epsilon} \]

\(\mathbf{X}\)はRで以下のように求められる。

head(model.matrix(M1))
##   (Intercept) LatAbs
## 1           1   1.33
## 2           1   5.42
## 3           1   5.42
## 4           1  28.51
## 5           1  28.51
## 6           1  28.51

パラメータは最小二乗法(ordinary least square)で求められる。最小二乗法は、残差の二乗和(実際の測定値と推定されたモデルによる予測値の差、ここでは\(\sum_i^{55}( y_i - \alpha + \beta \times OrbitalV_i)\))が最小になるようにパラメータを推定する方法である。

\(\mathbf{\beta}\)の推定値は数学的に以下のように求められる。なお、\(\mathbf{X^t}\)は\(\mathbf{X}\)の転置行列を、\(\mathbf{X^{-1}}\)はは\(\mathbf{X}\)の逆行列を表す。

\[ \mathbf{\hat{\beta}} = (\mathbf{X^t}\times \mathbf{X})^{-1} \times \mathbf{X^t} \times \mathbf{y} \tag{1.1} \]

なお、\(\mathbf{\hat{\beta}}\)は\(\mathbf{\beta}\)の推定値であることを表し、便宜的にここでは\(\mathbf{\hat{\beta}} = \begin{pmatrix}a\\b \end{pmatrix}\)とする。

Rでは\(\mathbf{\hat{\beta}}\)を以下のように求められ、lm関数で推定した場合と同じ推定値が得られることが分かる。

X <- model.matrix(M1)
solve(t(X) %*% X) %*% t(X) %*% hvs_b$OrbitalV
##                   [,1]
## (Intercept) 22.9106427
## LatAbs       0.0659754

よって、モデルによって推定された眼窩容量(\(\hat{y_i}\))は以下のように表せる。

\[ Fitted OrbitalV_i = 22.930 + 0.066 \times LatAbs_i \]

すなわち、モデルの推定値から得られた残差\(e_i\)は以下のように表せる。

\[ e_i = OrbitalV_i - a + b\times LatAbs_i \]

\(\mathbf{\hat{y}} = \mathbf{X}\times \mathbf{\beta}\)と書けるので、残差の行列\(\mathbf{e}\)は\(\mathbf{H} = (\mathbf{X^t}\times \mathbf{X})^{-1} \times \mathbf{X^t}\)とするとき以下のように書ける。

\[ \mathbf{e} = \mathbf{y} - \mathbf{\hat{y}} = \mathbf{y} - \mathbf{H\times y} = \mathbf{(1-H)\times y} \tag{1.2} \]

モデル式より、\(\mathbf{y}\)の分散共分散行列は\(\sigma^2 \times \mathbf{I}\)と表せるので(\(\mathbf{I}\)は単位行列)、\(\mathbf{e}\)の分散共分散行列は以下のように表せる。

\[ cov(\mathbf{e}) = \sigma^2 \times (\mathbf{1 - H}) \]

また、以下の値を標準化残差(standardized residuals)という。なお、\(H_{ii}\)は\(\mathbf{H}\)の\(i\)番目の対角成分を表す。

\[ e_i ^* = \frac{e_i}{\hat{\sigma} \sqrt{(1-H_{ii})}} \]

もしモデルが正しいとき、標準化残差は標準正規分布に従うので、その値のほとんどは-2から2の間に収まるはずである。

Rでは、残差と標準化残差を以下のように求められる。

e <- resid(M1)
estd <- rstandard(M1)

実際、1つのデータを除いて-2から2の範囲に収まっている(図1.4)。

data.frame(estd = estd) %>% 
  ggplot(aes(x = estd))+
  geom_histogram(binwidth = 0.3)+
  theme_bw()+
  theme(aspect.ratio = 1)+
  scale_y_continuous(breaks = seq(0,10,1))
標準化残差の分布

図1.4: 標準化残差の分布


1.5.1.1 外れ値のチェック

モデルに極端な外れ値がないかを調べるときには、LeverageとCook’s distanceが用いられることが多い。

leverageは\(\mathbf{H}\)の対角成分で0から1の値をとり、大きいほどそのデータポイントが結果に大きな影響を与える外れ値であることを示す。特に極端な値はないよう(図1.5)。

data.frame(h = diag(H),
           n = 1:nrow(hvs_b)) %>% 
  ggplot(aes(x = n, y = h))+
  geom_col(width = 0.3)+
  theme_bw()+
  theme(aspect.ratio = 1)+
  labs(x = "sample number", y = "Leverage")
Leveragevalue for each observation

図1.5: Leveragevalue for each observation

式(1.2)より、以下のように書ける。すなわち、\(i\)番目の観察に対してモデルから予測される値(fitted value)は元データの観測値の加重平均であり、重みは行列\(\mathbf{H}\)により与えられる。

\[ \hat{y_i} = H_{i1} \times y_1 + H_{i2} \times y_2 + \dots + H_{i55} \times y_55　\tag{1.3} \]

特に\(H_{ii}\)は観測値\(y_i\)がモデルからの予測値\(\hat{y_i}\)に与える影響の大きさを表している。そのため、\(H_{ii}\)が高いことは、その観測値がモデルの推定に大きな影響を及ぼしていることを示しているのである。

ある観測値がモデルの推定に影響を与えている度合いを表すのがCook’s distanceで、以下の式で表せる。\(\hat{y_{(i)}}\)は式(1.3)から\(H_{ii}\times y_i\)を除いたものである。また、\(p\)はモデルの回帰式の中のパラメータ数である。

\[ D_i = \frac{||\hat{y_i} - \hat{y_{(i)}}||}{p \times \hat{\sigma^2}} \]

Cook’s distanceはRで以下のように取得できる。

D <- cooks.distance(M1)

data.frame(cookD = D,
           n = 1:nrow(hvs_b)) %>% 
  ggplot(aes(x = n, y = cookD))+
  geom_col(width = 0.2)+
  theme_bw()+
  theme(aspect.ratio = 1)+
  labs(x = "sample number", y = "Cook's distance")
Cook's distance for each observation

図1.6: Cook’s distance for each observation

1.5.1.2 95%信頼区間と予測区間の算出

\(\mathbf{\beta}\)の推定値の分散共分散行列は\(\mathbf{y}\)の分散共分散行列が\(\mathbf{\sigma^2 \times I}\)であることを考えると、式(1.1)より以下のようになる。
\[ \begin{aligned} cov(\mathbf{\hat{\beta}}) &= \mathbf(X^t \times X)^{-1} \times \mathbf{X^t} \times cov(\mathbf{y}) \times \mathbf{X} \times (\mathbf{X^t} \times \mathbf{X})^{-1} \\ &= \sigma^2 \times (\mathbf{X^t} \times \mathbf{X})^{-1} \end{aligned} \]

\(\hat{\beta}\)の標準偏差は上式から得られる行列の対角成分の1/2乗である。この値は、lm関数を利用した結果(Std. Errorの下の数値)と一致する。

SE <-summary(M1)$sigma * sqrt(diag(solve(t(X) %*% X)))

SE
## (Intercept)      LatAbs 
##  0.51411731  0.01410749

95%信頼区間は、モデル式より\(\beta\)が正規分布に従うので以下のように求められる。

Z <- rbind(coef(M1) + 1.96*SE,
           coef(M1) - 1.96*SE) %>% 
  data.frame() 

rownames(Z) <- c("Upper bound", "Lower bound")

Z

ただし、これはサンプルサイズが十分に大きいときのみ成り立つ。実際は、\(\hat{\beta}\)は自由度\(55-2\)(サンプル数 - パラメータ数)のt分布に従うので、より正確に95%信頼区間を求めるには、1.96ではなく2.005746…を用いる必要がある。

qt(1- 0.05/2, df = 55-2)
## [1] 2.005746

95％信頼区間は、もし100回同様の方法で実験/観察を行ってそれぞれについて95%信頼区間を算出するとき、そのうち95個には真の値が含まれていることを表す。95%信頼区間に0が含まれていないとき、\(\hat{\beta}\)が有意に0とは違うということができる。

同様に、モデルに基づいた予測値の分散共分散行列は以下のように求められ、その対角成分の1/2乗が予測値の標準誤差(SE)になる。　　

\[ cov(\mathbf{\hat{y}}) = \mathbf{X} \times cov(\mathbf{\hat{\beta}}) \times \mathbf{X^t} \tag{1.4} \]

予測値の95%信頼区間はある標高(LatAbs)に対して100回データをサンプリングすれば95個のデータが含まれる範囲を表し、予測値が自由度55-2のt分布に従うので、\(予測値 ± 2.0057 \times SE\)で求められる。

一方で95%予測区間は新たにデータをサンプリングしたときにデータの95%が収まる範囲を指す。予測区間を求める際の標準誤差には式(1.4)に\(\hat{\sigma^2}\)を足したものを用いればよい。

Rでは予測値と95%信頼区間、95%予測区間は以下のように求められる。

## 係数の分散共分散行列  
covbeta <- vcov(M1)
data <- data.frame(LatAbs = seq(0,65,length.out = 100))
X <- model.matrix(~LatAbs, data = data)

t <- qt(1-0.05/2, df = 55-2)
  
data %>% 
  ## 予測値
  mutate(p = X %*% coef(M1)) %>% 
  ## se(信頼区間)
  mutate(se.ci = sqrt(diag(X %*% covbeta %*% t(X)))) %>% 
  ## se(予測区間)
  mutate(se.pi = sqrt(diag(X %*% covbeta %*% t(X)) + summary(M1)$sigma^2)) %>% 
  mutate(ci_upper = p + t*se.ci,
         ci_lower = p - t*se.ci,
         pi_upper = p + t*se.pi,
         pi_lower = p - t*se.pi) -> pred

図示すると以下のようになる(図1.7)。

pred %>% 
  ggplot(aes(x = LatAbs, y = p))+
  geom_line()+
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper),
              fill = "grey21", alpha = 0.5)+
  geom_ribbon(aes(ymin = pi_lower, ymax = pi_upper),
              fill = "grey72", alpha = 0.5)+
  geom_point(data = hvs_b,
             aes(y = OrbitalV))+
  theme(aspect.ratio = 1)+
  theme_bw()+
  labs(y = "OrbitalV")
95%信頼区間と予測区間

図1.7: 95%信頼区間と予測区間

predict関数を用いて簡単に求めることができる。

predict(M1,
        newdata = data.frame(LatAbs = seq(0,65,0.2)),
        interval = "confidence") %>% 
  data.frame() %>% 
  bind_cols(data.frame(LatAbs = seq(0,65,0.2))) -> conf_M1

predict(M1,
        newdata = data.frame(LatAbs = seq(0,65,0.2)),
        interval = "prediction") %>% 
  data.frame() %>% 
  bind_cols(data.frame(LatAbs = seq(0,65,0.2))) -> pred_M1

hvs_b %>% 
  ggplot(aes(x = LatAbs, y = OrbitalV))+
  geom_line(data = conf_M1,
              aes(y = fit))+
  geom_ribbon(data = conf_M1,
              aes(y = fit, ymin = lwr, ymax = upr),
              fill = "grey21", alpha = 0.5)+
  geom_ribbon(data = pred_M1,
              aes(y = fit, ymin = lwr, ymax = upr),
              fill = "grey72", alpha = 0.5)+
  geom_point()+
  theme(aspect.ratio = 1)+
  theme_bw()
95%信頼区間と予測区間(predict関数を使用)

図1.8: 95%信頼区間と予測区間(predict関数を使用)

1.5.2 Multiple linear regression

それでは、2つ以上の変数を含めたモデリングを行う。データ探索や先行研究の知見から、モデルにはLatAbs、CC、FMを説明変数として入れ、2変数の交互作用を全ての組み合わせについて含めた。

\[ \begin{aligned} OrbitalV_i &= \alpha + \beta_1 \times LatAbs_i + \beta_2 \times CC_i + \beta_3 \times FM_i \\ & + \beta_4 \times LatAbs_i \times CC_i\\ & + \beta_5 \times LatAbs_i \times FM_i \\ & + \beta_6 \times CC_i \times FM_i + \epsilon_i\\ \epsilon_i &\sim N(0,\sigma^2) \end{aligned} \]

1.5.3 Fitting the model in R and estimate parameters

モデルのパラメータは、先ほどと同様にlm関数で推定できる。推定は前節で行ったのと全く同じ方法(最小二乗法)で行う。

M2 <- lm(OrbitalV ~ LatAbs + CC + FM + LatAbs:CC + LatAbs:FM + CC:FM,
         data = hvs_b)

結果は以下の通り。Estimateはパラメータの推定値を、t valueは5%水準でパラメータが有意に0と異なっているかを判断する際に用いられる。P値(Pr(>|t|))を見ると、有意な変数は一つもなかった。

print(summary(M2), digits = 3, signif.stars = FALSE)
## 
## Call:
## lm(formula = OrbitalV ~ LatAbs + CC + FM + LatAbs:CC + LatAbs:FM + 
##     CC:FM, data = hvs_b)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -3.79  -1.50  -0.23   1.40   4.63 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.19e+01   2.24e+01    0.97    0.335
## LatAbs      -2.51e-01   1.60e-01   -1.56    0.124
## CC           9.37e-04   1.65e-02    0.06    0.955
## FM           2.16e-03   3.50e-02    0.06    0.951
## LatAbs:CC    2.37e-04   1.30e-04    1.83    0.074
## LatAbs:FM   -4.49e-05   1.57e-04   -0.29    0.776
## CC:FM       -1.56e-06   2.50e-05   -0.06    0.951
## 
## Residual standard error: 2.14 on 47 degrees of freedom
## Multiple R-squared:  0.427,  Adjusted R-squared:  0.354 
## F-statistic: 5.84 on 6 and 47 DF,  p-value: 0.000131

1.6 Finding the optimal model

さて、先ほどのモデルでは1つも有意に0と異なるパラメータ(\(\alpha, \beta_1 \sim \beta_6\))はなかった(= 目的変数に有意な影響を持つ説明変数変数がなかった)。このようなとき、選択肢がいくつかある。

  1. そのままのモデルを採用し、全ての交互作用が5%水準では有意ではなかったと報告する。
  2. AICを利用して古典的なモデル選択を行う。
  3. 仮説検定の結果に基づいて変数選択を行う(効果のなさそうな変数を外す)。
  4. 交互作用についてのみモデル選択を行う。
  5. 情報理論的アプローチを用い、モデル平均化などを行う(c.f., Burnham et al. (2011) )。

ここでは、 Pearce and Dunbar (2012) に従い、AICを用いたモデル選択を行うことにする。

Rでは、step関数を用いることでステップワイズ法(AICが最も低くなるまで説明変数を1つずつ増減させる方法)を用いた変数選択を行うことができる。分析の結果、LatAbs、CC、これらの交互作用のみを含むモデルが最もAICが低い(= 予測精度が高い)と判断された。

step(M2)
## Start:  AIC=88.54
## OrbitalV ~ LatAbs + CC + FM + LatAbs:CC + LatAbs:FM + CC:FM
## 
##             Df Sum of Sq    RSS    AIC
## - CC:FM      1    0.0177 214.75 86.547
## - LatAbs:FM  1    0.3756 215.11 86.637
## <none>                   214.74 88.543
## - LatAbs:CC  1   15.2596 230.00 90.250
## 
## Step:  AIC=86.55
## OrbitalV ~ LatAbs + CC + FM + LatAbs:CC + LatAbs:FM
## 
##             Df Sum of Sq    RSS    AIC
## - LatAbs:FM  1    0.4172 215.17 84.652
## <none>                   214.75 86.547
## - LatAbs:CC  1   15.2501 230.00 88.252
## 
## Step:  AIC=84.65
## OrbitalV ~ LatAbs + CC + FM + LatAbs:CC
## 
##             Df Sum of Sq    RSS    AIC
## - FM         1    0.5409 215.71 82.788
## <none>                   215.17 84.652
## - LatAbs:CC  1   16.8268 232.00 86.718
## 
## Step:  AIC=82.79
## OrbitalV ~ LatAbs + CC + LatAbs:CC
## 
##             Df Sum of Sq    RSS    AIC
## <none>                   215.71 82.788
## - LatAbs:CC  1     17.24 232.95 84.940
## 
## Call:
## lm(formula = OrbitalV ~ LatAbs + CC + LatAbs:CC, data = hvs_b)
## 
## Coefficients:
## (Intercept)       LatAbs           CC    LatAbs:CC  
##   2.323e+01   -2.552e-01   -5.829e-05    2.201e-04

そこで、選ばれた変数のみを用いたモデルを作成し、分析を行う。

M3 <- lm(OrbitalV ~ LatAbs*CC, data = hvs_b)

分析の結果は以下のとおりである。交互作用項の係数の推定値のP値が0.051であり、わずかに交互作用項の影響があることが示唆された?

print(summary(M3), digits = 3, signif.stars = FALSE)
## 
## Call:
## lm(formula = OrbitalV ~ LatAbs * CC, data = hvs_b)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.759 -1.504 -0.098  1.513  4.588 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.32e+01   5.23e+00    4.44  4.9e-05
## LatAbs      -2.55e-01   1.53e-01   -1.67    0.101
## CC          -5.83e-05   3.94e-03   -0.01    0.988
## LatAbs:CC    2.20e-04   1.10e-04    2.00    0.051
## 
## Residual standard error: 2.08 on 50 degrees of freedom
## Multiple R-squared:  0.425,  Adjusted R-squared:  0.39 
## F-statistic: 12.3 on 3 and 50 DF,  p-value: 3.81e-06

1.7 Degree of freedom

重回帰分析では、　モデルの自由度は回帰式のパラメータ数なので4である。しかし、このような自由度の求め方はGAMではできない。一般に、パラメータ数は\(\mathbf{H}\)の対角成分の和(= trace)で求めることができる。

\[ p = \sum_{i =1} ^{55} \mathbf{H_{ii}} \]

確かに4になっている。

X <- model.matrix(M3)
H <- X %*% solve(t(X) %*% X) %*% t(X)

sum(diag(H))
## [1] 4

1.8 Model validation

最後に、作成したモデルが前提を満たしているかをチェックする。残差には普通の残差のほかに、標準化残差、スチューデント化残差などがあるが、ここでは標準化残差を用いる。

  1. 等分散性が成り立つか (→ モデルに基づく予測値と残差をプロットする)
  2. モデルがデータに当てはまっているか(線形性が成り立つか)、残差の独立性があるか (→ 残差とモデルに含まれる説明変数、モデルに含まれない説明変数の関係をプロットする)
  3. データの独立性があるか(→ 自己相関があるかを確認する)
  4. 残差が仮定した分布(今回は正規分布)に従うかを確認する(→ QQプロットを書く)
  5. モデルへの影響が大きいデータがないか確認する(→ Cook’s distanceやLeverage)

まず、モデルによる予測値と標準化残差の関係をプロットする(図1.9)。このときプロットにパターンがあってはいけない(e.g., 広がっていく、ばらつきが不均等など)。もし前提を満たしていれば、0を中心に均等にばらつくはずである。図からは、明確なパターンは見られない(横軸が26~27で負の残差がほとんどないなどを除けば)。

data.frame(fitted = fitted(M3),
           resstd = rstandard(M3)) %>% 
  ggplot(aes(x = fitted, y = resstd))+
  geom_point()+
  geom_hline(yintercept = 0)+
  theme_bw()+
  theme(aspect.ratio = 1)+
  labs(x = "Fitted values", y = "Standardized residuals")
Standardized residuals versus fitted values to assess homogeneity

図1.9: Standardized residuals versus fitted values to assess homogeneity


続いて、Cook`s distanceを用いて影響の大きいデータがないか確認する。Cook’s distanceは通常1を超えると影響が大きいと判断される2。ほんもでるではそのような値をとる観測値はない(図1.10)。

data.frame(cook = cooks.distance(M3),
           n = 1:nrow(hvs_b)) %>% 
  ggplot(aes(x = n, y = cook))+
  geom_col(width = 0.2)+
  geom_hline(yintercept = 1, linetype = "dotted")+
  theme_bw()+
  theme(aspect.ratio = 1)+
  labs(x = "Observations", y = "Cook's distance")
Cook's disance for each observation

図1.10: Cook’s disance for each observation

残差の正規性は、ヒストグラムとQQプロットを基に判断される。QQプロットとは、「得られたデータ(今回の場合は標準化残差)と理論分布(今回の場合は標準正規分布)を比較し、その類似度を調べるためのグラフ」である。詳細についてはこちら。もし類似度が高ければ、点が直線上に乗る。そこまで大きくは外れていなさそう。

qqnorm(rstandard(M3))
qqline(rstandard(M3))
QQplot for M3

図1.11: QQplot for M3

続いて、(モデルに入れていない説明変数を含めた)各連続変数と標準化残差の関連を図示する(図1.12)。いずれも関連はないようで、これはモデルによって説明できない部分(= 残差)はこうした変数とも関連がないことを示している。

hvs_b %>% 
  mutate(resstd = rstandard(M3)) %>% 
  select(CC, LatAbs, FM, Illuminance, Temperature, resstd) %>% 
  pivot_longer(cols = 1:5, names_to = "var", values_to = "values") %>% 
  ggplot(aes(x = values, y = resstd))+
  geom_point()+
  geom_smooth(method = "lm", color = "black")+
  facet_rep_wrap(~var, repeat.tick.labels = TRUE, scales = "free_x")+
  theme_bw()+
  theme(aspect.ratio = 1)+
  labs(x = "Explanatory variables", y = "Standardized residuals")
Multiple scatterplots of the standardized residuals versus each continuous variables

図1.12: Multiple scatterplots of the standardized residuals versus each continuous variables


これは、離散的な変数についても同様のようである(図1.13)。

hvs_b %>% 
  mutate(resstd = rstandard(M3)) %>% 
  select(fPopulation, fGender, resstd) %>% 
  pivot_longer(1:2, names_to = "var", values_to = "values") %>% 
  ggplot(aes(x = values, y = resstd))+
  geom_boxplot()+
  facet_rep_wrap(~var, scales = "free_x")+
  theme_bw()+
  theme(aspect.ratio = 1)+
  labs(x = "Explanatory variables", y = "Standardized residuals")
Boxplot of standardized residuals versus discrete variable

図1.13: Boxplot of standardized residuals versus discrete variable


自己相関も大きくなさそうである。

acf(hvs_b$OrbitalV)

1.9 Model interpretation

モデルの推定値に基づくと、眼窩容量の予測値は以下のように書ける。

\[ OrbitalV_i = 23.23 -0.25 \times LatAbs_i -0.000058 \times CC_i + 0.00022 \times LatAbs_i \times CC_i \]

モデルに基づく回帰平面を描くと以下のようになる(図1.14)。

data_M3 <- crossing(LatAbs = seq(0.02,65, length.out = 100),
                    CC = seq(1100,1700, length.out = 100))

pred <- predict(M3, newdata = data_M3) %>% 
  data.frame() %>% 
  rename(pred = 1) %>% 
  bind_cols(data_M3) %>% 
  pivot_wider(names_from = CC, values_from = pred) %>% 
  select(-1) 

plot_ly(hvs_b,
        x = ~LatAbs,
        y = ~CC,
        z = ~OrbitalV,
        type = "scatter3d",
        size = 2) %>% 
  add_trace(z = as.matrix(pred),
            x = seq(0.02,65, length.out = 100),
            y = seq(1100,1700, length.out =100),
            type = "surface") -> p

p

図1.14: 3d plot of the fitted surface

1.10 What to do if things go wrong

もしモデルが前提を満たさなかったらどうすればよいだろうか?本節では特にどのようなときにGAMを適用すべきかを解説する。

等分散性の仮定が満たされないとき、以下の選択肢を検討する必要がある。

  1. 目的変数に変数変換を施す
  2. 一般化最小二乗法(GLS)を用いる
  3. さらに説明変数や交互作用を加える
  4. 目的変数の分布として他の分布を用いる(ガンマ分布など)

3つ目の選択肢において、非線形のパターンも許容するように拡張すると、GAMが使えるようになる。これは、残差と説明変数が何らかのパターンを示した時にも有効である(= 関係が線形でない可能性があるため)。

References

Burnham KP, Anderson DR, Huyvaert KP (2011) AIC model selection and multimodel inference in behavioral ecology: Some background, observations, and comparisons. Behav Ecol Sociobiol 65:23–35
Dunn PK, Smyth GK (2018) Generalized linear models with examples in R. Springer New York
Pearce E, Dunbar R (2012) Latitudinal variation in light levels drives human visual system size. Biol Lett 8:90–93
Zuur AF, Ieno EN, Elphick CS (2010) A protocol for data exploration to avoid common statistical problems. Methods Ecol Evol 1:3–14

  1. 説明変数(独立変数)とは、物事の原因となっている変数のこと、目的変数(応答変数)とは説明変数の影響を受けて発生した結果となっている変数のことである。今回の場合は、様々な標高など(= 説明変数)が眼窩容量(= 目的変数)に与える影響をモデリングする。↩︎

  2. 　より正確には、F分布の50パーセンタイルが基準になるが、多くの場合その値は1に近い(Dunn and Smyth 2018)。↩︎