lm 関数を利用して、複数の独立変数で 1 つの従属変数を説明するモデル

R による重回帰分析

回帰分析は多変量解析法の一つである。1 つの従属変数 y を、複数の独立変数 xi (i = 1, 2, ..., n) でモデル化し、それらの関係や傾向を分析する方法である。そのモデル式は、次のように書くことができる。

\[ y = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + \cdots + \beta_{n}x_{n} \]

lm 関数を利用した重回帰分析

R の lm を利用して、重回帰分析を行なってみる。ここではサンプルデータとして、哺乳類の睡眠時間を調査したデータを利用する。このデータには、動物の名前(Species)、体重(BodyWt)、睡眠時間(TotalSleep)、寿命(LifeSpan)などが記載されている。ここでは、寿命を従属変数として、体重と睡眠時間を独立変数として、重回帰分析を行う。その結果をもとに、寿命は、体重あるいは睡眠時間のどちらかに強く影響されるかを考えてみる。

x <- read.table("https://stats.biopapyrus.jp/data/sleep_in_mammals.txt", header = TRUE)
head(x)
##                  Species   BodyWt BrainWt NonDreaming Dreaming TotalSleep LifeSpan Gestation Predation Exposure Danger
## 1        Africanelephant 6654.000  5712.0          NA       NA        3.3     38.6       645         3        5      3
## 2 Africangiantpouchedrat    1.000     6.6         6.3      2.0        8.3      4.5        42         3        1      3
## 3              ArcticFox    3.385    44.5          NA       NA       12.5     14.0        60         1        1      1
## 4   Arcticgroundsquirrel    0.920     5.7          NA       NA       16.5       NA        25         5        2      3
## 5          Asianelephant 2547.000  4603.0         2.1      1.8        3.9     69.0       624         3        5      4
## 6                 Baboon   10.550   179.5         9.1      0.7        9.8     27.0       180         4        4      4

LifeSpan を従属変数とし、BodyWt と TotalSleep を独立変数として、回帰モデルの各係数を lm 関数で計算する。

\[ LifeSpan = \beta_{0} + \beta_{1}BodyWt + \beta_{2}TotalSleep \]
LifeSpan <- x$LifeSpan
BodyWt <- x$BodyWt
TotalSleep <- x$TotalSleep

解析を行う前に、まず、この 3 つのデータの分布を確認してみる。もっとも簡単な方法として、3 つデータに対して、ヒストグラムを描くことである。

par(mfrow = c(1, 3))
hist(LifeSpan)
hist(BodyWt)
hist(TotalSleep)
睡眠時間、寿命、体重のデータの分布

このデータの分布を確認すると、LifeSpan と BodyWt が、正規分布から著しく離れていることがわかる。回帰分析は、基本的にデータが正規分布に従う場合のみに利用できる。そこで、ここで対数変換を行なって、LifeSpan と BodyWt を正規分布に近づけてみる。(※ 対数変換を行わずに、GLM を利用して解析した方がいいかも。)

LifeSpan.log <- log10(LifeSpan)
BodyWt.log <- log10(BodyWt)

par(mfrow = c(1, 3))
hist(LifeSpan.log)
hist(BodyWt.log)
hist(TotalSleep)
対数変換化された睡眠時間、寿命、体重のデータの分布

LifeSpan と BodyWt のデータ分布がだいぶ正規分布に近づいた。ここで、回帰分析を行なってみる。体重と睡眠時間は、単位が異なるので、解析する前に両者が同じ比重となるように標準化する。

BodyWt.log.std <- scale(BodyWt.log)
TotalSleep.std <- scale(TotalSleep)

m <- lm(LifeSpan.log ~ TotalSleep.std + BodyWt.log.std)
summary(m)
##
## Call:
## lm(formula = LifeSpan.log ~ TotalSleep.std + BodyWt.log.std)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max
## -0.57995 -0.19894 -0.05503  0.16127  0.84891
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)     1.11941    0.04121  27.167  < 2e-16 ***
## TotalSleep.std -0.02846    0.04829  -0.589    0.558
## BodyWt.log.std  0.27866    0.04873   5.718 5.65e-07 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 0.3024 on 51 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.5033,	Adjusted R-squared:  0.4838
## F-statistic: 25.84 on 2 and 51 DF,  p-value: 1.781e-08

Cofficients 項の Estimate 列に、回帰モデルの係数が書かれている。この場合、この回帰モデルは、LifeSpan = -0.02846×TotalSleep + 0.27866×BodyWt + 1.11941 となる。この解析結果から、体重が増えるほど寿命が長くなり、また、睡眠時間が長いほど寿命が短くなる傾向がみられる。また、睡眠時間よりも、体重が寿命に与える影響が強いと読み取ることができる。

寿命と体重、寿命と睡眠時間の散布図を描くと、次のようになる。回帰分析で得られた結果と一致することが確認できる。

par(mfrow = c(1, 2))
plot(LifeSpan.log, BodyWt.log.std, ylim = c(-3, 3))
plot(LifeSpan.log, TotalSleep.std, ylim = c(-3, 3))
睡眠時間、寿命、体重のデータの相関関係