10  回帰分析の基礎

10.1 線形回帰 – 散布図への直線の当てはめ

前章では2変数の関係を表すために散布図を作り，相関係数を求めましたが，この章では，背後にモデルとして原因と結果の関係をもつと仮定し，そのモデルを表すために直線を当てはめる回帰分析について学びます．

分析に使うデータとして，2006年から2022年の東京証券取引所のプライム市場に上場している会社の決算データを使います。 売上高と広告宣伝費や研究開発費の関係を分析するために，研究開発費が1億円以上，広告宣伝費が1000万円以上の会社に限定し，決算月数が12ヶ月の会社を対象にします。

df <- read_csv("data/adv_2023.csv")
df <- df %>%
  filter(決算月数 == 12 & 研究開発費 >= 100 & 広告宣伝費 > 10)
glimpse(df)
Rows: 4,442
Columns: 14
$ 日経会社コード <chr> "0000001", "0000001", "0000003", "0000003", "0000003", …
$ 企業名称       <chr> "極洋", "極洋", "日本水産", "日本水産", "日本水産", "日…
$ 決算期         <chr> "2006/03", "2007/03", "2006/03", "2007/03", "2008/03", …
$ 決算種別       <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,…
$ 連結基準       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ 決算月数       <dbl> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,…
$ 業種           <dbl> 235341, 235341, 235341, 235341, 235341, 235341, 235341,…
$ 資産合計       <dbl> 65049, 66459, 384819, 404173, 396739, 385462, 383924, 4…
$ 売上高         <dbl> 152899, 157088, 539653, 552871, 533970, 505250, 481574,…
$ 販管費         <dbl> 13702, 14455, 95566, 98200, 100394, 98413, 99938, 10490…
$ 広告宣伝費     <dbl> 304, 279, 2699, 2569, 2953, 2568, 2636, 3160, 3009, 288…
$ 拡販費         <dbl> 111, 158, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ 研究開発費     <dbl> 193, 188, 3083, 3377, 3718, 3803, 3994, 4499, 4809, 361…
$ 設備投資額     <dbl> 897, 1841, 17186, 16031, 19105, 28872, 21121, 18633, 16…

では売上高と広告宣伝費の散布図を書いてみます。

g <- ggplot(df) + aes(x = 広告宣伝費, y = 売上高) + geom_point() + mystyle
print(g)

散布図を見ると，売上高も広告宣伝費もばらつきが大きいので，とりあえず対数変換してみます。

df <- df %>%
  mutate(
    log_広告宣伝費 = log(広告宣伝費),
    log_売上高 = log(売上高),
    log_研究開発費 = log(研究開発費)
  )
glimpse(df)
Rows: 4,442
Columns: 17
$ 日経会社コード <chr> "0000001", "0000001", "0000003", "0000003", "0000003", …
$ 企業名称       <chr> "極洋", "極洋", "日本水産", "日本水産", "日本水産", "日…
$ 決算期         <chr> "2006/03", "2007/03", "2006/03", "2007/03", "2008/03", …
$ 決算種別       <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,…
$ 連結基準       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ 決算月数       <dbl> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,…
$ 業種           <dbl> 235341, 235341, 235341, 235341, 235341, 235341, 235341,…
$ 資産合計       <dbl> 65049, 66459, 384819, 404173, 396739, 385462, 383924, 4…
$ 売上高         <dbl> 152899, 157088, 539653, 552871, 533970, 505250, 481574,…
$ 販管費         <dbl> 13702, 14455, 95566, 98200, 100394, 98413, 99938, 10490…
$ 広告宣伝費     <dbl> 304, 279, 2699, 2569, 2953, 2568, 2636, 3160, 3009, 288…
$ 拡販費         <dbl> 111, 158, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ 研究開発費     <dbl> 193, 188, 3083, 3377, 3718, 3803, 3994, 4499, 4809, 361…
$ 設備投資額     <dbl> 897, 1841, 17186, 16031, 19105, 28872, 21121, 18633, 16…
$ log_広告宣伝費 <dbl> 5.717028, 5.631212, 7.900637, 7.851272, 7.990577, 7.850…
$ log_売上高     <dbl> 11.93753, 11.96456, 13.19868, 13.22288, 13.18809, 13.13…
$ log_研究開発費 <dbl> 5.262690, 5.236442, 8.033658, 8.124743, 8.220941, 8.243…

対数変換した売上高と広告宣伝費の関係を示す、散布図を作成します。

g <- ggplot(df) + aes(x = log_広告宣伝費, y = log_売上高) + geom_point() + mystyle
print(g)

キレイに正の相関があるようにみえる散布図ができあがりました。 相関係数を計算してみると，

cor(df$log_広告宣伝費, df$log_売上高)
[1] 0.7530297

となり，非常に高い正の相関があることが分かりました。

この散布図に回帰直線を追加するには，geom_smoothを使います。

g <- g + geom_smooth(method = "lm", se = FALSE, color = "red") + xlab("広告宣伝費の対数") + ylab("売上高の対数") + mystyle
print(g)

この赤い直線を回帰直線(regression line)といい，原因となる変数y(ここでは広告宣伝費)と結果となる変数x(ここでは売上高)の線形関係を表す次のようなモデルとなります。

y = a + bx

先の広告宣伝費と売上高の散布図に引いた直線だと，

res01 <- lm(log_売上高 ~ log_広告宣伝費, data = df)
summary(res01)

\text{売上高の対数} = 7.770869 + 0.550978 \times 広告宣伝費の対数

10.2 最小二乗法

先ほどの求めた回帰直線の切片や傾きの値を求める方法の1つに最小二乗法(least squares method)があります。

最小二乗法は、観測値y_iとモデル上の予測値\hat{y}_iの差を残差(residual)と定義し、この残差の二乗の和を最小にするような回帰直線を求める方法です。

サイズNの標本(y_i, x_i)_{i \in N}をデータとして持っているとしましょう。 ここでは、広告宣伝費と売上高のデータがN組あるということです。

このデータの関係を線形で表すモデルとして、\hat y = \hat a + \hat b xを考えます。 観測値はy_i = a + bx_i + e_iです。 したがって、残差e_iはy_i - \hat y_iとなります。

data <- data.frame(
  x = c(1, 2, 3, 4, 5),
  y = c(3, 10, 9, 12, 18)
)
model <- lm(y ~ x, data = data)

# 回帰直線と残差を含むグラフを作成する
ggplot(data, aes(x, y)) +
  geom_point(size = 2) +  # 散布図のプロット
  geom_smooth(method = "lm", se = FALSE, fullrange = TRUE) +  # 回帰直線のプロット
  geom_segment(aes(xend = x, yend = predict(model), color = "Residuals")) +  # 残差のプロット
#  geom_hline(yintercept = 0, linetype = "dashed") +  # 残差0の水平線
  labs(title = "回帰直線と残差", x = "x", y = "y") +  # グラフのタイトルと軸ラベル
  scale_color_manual(values = c("Residuals" = "blue")) + # 残差の色を指定する
  xlim(0, 6) + ylim(0, 20) + mystyle # x軸とy軸の範囲を指定する

この観測値を表す点から回帰直線に向けて書かれた垂線の長さが残差eとなります。 この二乗和を最小にするような回帰直線を求めることが最小二乗法です。 残差の二乗和は、次のように計算できます。

\begin{aligned} \sum_{i \in N} e_i^2 &= \sum_{i \in N} (y_i - \hat y_i)^2 \\ &= \sum_{i \in N} (y_i - \hat a - \hat b x_i)^2 \\ &= \sum_{i \in N} (y_i - \bar y + \bar y - \hat a - \hat b x_i)^2 \\ &= \sum_{i \in N} (y_i - \bar y)^2 + \sum_{i \in N} (\bar y - \hat a - \hat b x_i)^2 + 2 \sum_{i \in N} (y_i - \bar y)(\bar y - \hat a - \hat b x_i) \\ &= \sum_{i \in N} (y_i - \bar y)^2 + \sum_{i \in N} (\bar y - \hat a - \hat b x_i)^2 + 0 \\ \end{aligned}

この最小二乗法による回帰直線の傾きと切片を推定するには、lm関数を使います。 lm()関数は引数にformulaとdataを取ります。

  • formula : 回帰式を指定します。y ~ xとすると、yをxで回帰します。
  • data : データフレームを指定します。
res01 <- lm(log_売上高 ~ log_広告宣伝費, data = df)

推定結果を代入したresの中を詳細に見るには、summary関数を使います。

summary(res01)

Call:
lm(formula = log_売上高 ~ log_広告宣伝費, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6125 -0.7040 -0.0514  0.7580  3.5226 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    7.770869   0.056878  136.62   <2e-16 ***
log_広告宣伝費 0.550978   0.007225   76.26   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.077 on 4440 degrees of freedom
Multiple R-squared:  0.5671,    Adjusted R-squared:  0.567 
F-statistic:  5815 on 1 and 4440 DF,  p-value: < 2.2e-16

重要な情報が表示されますが、いろいろ表示されすぎて見づらいので、modelsummary関数を使って見やすくします。

library(modelsummary)
modelsummary(res01)
 (1)
(Intercept) 7.771
(0.057)
log_広告宣伝費 0.551
(0.007)
Num.Obs. 4442
R2 0.567
R2 Adj. 0.567
AIC 13272.0
BIC 13291.2
Log.Lik. -6632.994
F 5815.313
RMSE 1.08

いろいろオプションを追加して、もっと見やすい表にします。

cm  <-  c(
    "(Intercept)" = "切片", # 切片のラベルを変更する
    "log_広告宣伝費" = "対数広告宣伝費" # 傾きのラベルを変更する
)
msummary(res01,
  stars = TRUE, # p値の有意性を星で表示する
  fmt = '%.2f',
  coef_map = cm,
  gof_omit = "AIC|BIC|Log.Lik.",
  output = 'html'
)
 (1)
切片 7.77***
(0.06)
対数広告宣伝費 0.55***
(0.01)
Num.Obs. 4442
R2 0.567
R2 Adj. 0.567
F 5815.313
RMSE 1.08
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

となり、切片が7.77、傾きが0.55と推定されました。 独立変数も従属変数も対数変換をしているため、解釈としては、 広告宣伝費が1\%上昇すると、売上が0.55\%増加する、という統計的に正の関係があるといえます。

10.3 単回帰と重回帰

今までは説明変数が1つだけでしたが、複数の要因が応答変数に影響を与える場合もあります。その場合は、複数の説明変数を使った重回帰分析を行います。

y_i = \beta_0 + \beta _1 x_1 + \beta_2x_2 + e_i

という回帰式を考えます。 この式は、x_1とx_2の2つの説明変数が応答変数yに影響を与えることを表す回帰モデルとなっています。 先の例で用いたデータを使って、広告宣伝費と研究開発費が与える影響についてみてみましょう。

まずは、単回帰と同様に、散布図と回帰平面を描いてみます。 単回帰では、説明変数が1つであったため、応答変数との関係を2次元平面上にプロットすることができましたが、今回の重回帰では説明変数2個と応答変数の関係となるため、3次元空間にプロットする必要があります。 そのため、回帰直線では無く、回帰平面を描くことになります。

# install.packages("rgl")
library(rgl) # 3Dグラフのパッケージ
library(knitr) # グラフの表示を操作
knitr::knit_hooks$set(webgl = hook_webgl) # おまじない

# 回帰モデルを推定する
model <- lm(log_売上高 ~ log_広告宣伝費 + log_研究開発費, data = df)

# プロットする範囲の値を指定する
x1_range <- range(df$log_広告宣伝費)
x2_range <- range(df$log_研究開発費)

# # メッシュグリッドを作成する
grid <- expand.grid(
  log_広告宣伝費 = seq(x1_range[1], x1_range[2], length.out = 50),
  log_研究開発費 = seq(x2_range[1], x2_range[2], length.out = 50)
  )
# メッシュグリッド上での予測値を計算する
grid$log_売上高 <- predict(model, newdata = grid)

# メッシュグリッドの座標を行列に変換する
x1_matrix <- matrix(grid$log_広告宣伝費, nrow = 50)
x2_matrix <- matrix(grid$log_研究開発費, nrow = 50)
y_matrix <-  matrix(grid$log_売上高,    nrow = 50)

# グラフを出力
plot3d(df$log_広告宣伝費, df$log_研究開発費, df$log_売上高,
 xlab = "対数広告宣伝費", ylab = "対数研究開発費", zlab = "対数売上高"
 ) # 3D散布図を描く
surface3d(x1_matrix, x2_matrix, y_matrix, alpha = 0.5, color = "red") # 回帰平面を描く

10.4 会計データを使った重回帰

では実際に会計データを使った重回帰分析を行ってみましょう。 売上高yを広告宣伝費x_1と研究開発費x_2で説明する、つまり

\log y_i = \beta_0 + \beta_1 \log x_{1,i} + \beta_2 \log x_{2,i} + e_i

という線形回帰モデルを考えます。 重回帰式の推定は、単回帰式と同様にlm関数を使い、複数の説明変数を+でつなげて指定します。

res02 <- lm(log_売上高 ~ log_広告宣伝費 + log_研究開発費, data = df)
summary(res02)

Call:
lm(formula = log_売上高 ~ log_広告宣伝費 + log_研究開発費, 
    data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9369 -0.5689 -0.0217  0.5057  3.5677 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    6.298486   0.055753  112.97   <2e-16 ***
log_広告宣伝費 0.323990   0.007583   42.73   <2e-16 ***
log_研究開発費 0.402482   0.008479   47.47   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8775 on 4439 degrees of freedom
Multiple R-squared:  0.7128,    Adjusted R-squared:  0.7127 
F-statistic:  5509 on 2 and 4439 DF,  p-value: < 2.2e-16

回帰係数だけを知りたいときは、coef()関数を使うと便利です。 ここでは単回帰の例と同様に、modelsummaryパッケージを使って結果をまとめてみましょう。

results <- list() # 結果を格納するリストを作成
results[["単回帰"]] <- res01 # 単回帰の結果
results[["重回帰"]] <- res02 # 重回帰の結果

modelsummary(results, # 結果を表にして出力
  stars = TRUE, # p値の有意性を星で表示する
  fmt = '%.2f', # 小数点以下2桁で表示する
  gof_omit = "AIC|BIC|Log.Lik.", # 不要な指標を除去
  output = 'html' # 出力形式をhtmlにする
  )
単回帰 重回帰
(Intercept) 7.77*** 6.30***
(0.06) (0.06)
log_広告宣伝費 0.55*** 0.32***
(0.01) (0.01)
log_研究開発費 0.40***
(0.01)
Num.Obs. 4442 4442
R2 0.567 0.713
R2 Adj. 0.567 0.713
F 5815.313 5509.365
RMSE 1.08 0.88
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

広告宣伝費の対数の回帰係数が、単回帰のときの0.55から0.32に減少しています。 これは、研究開発費の対数が追加されたことで、広告宣伝費が売上高を説明する部分が減少したことを意味します。 そして、研究開発費の対数の回帰係数は0.40となっており、広告宣伝費とともに統計的に有意な値となっています。 つまり、研究開発費も売上高を説明する要因となっていることがわかります。

10.5 単回帰と重回帰の違い

単回帰と重回帰の違いを、実際にデータを使って確認してみましょう。 つまり、以下の3つの回帰モデルの関係について考えてみます。

\begin{aligned} y &= \beta_0 + \beta _1 x_1 + \beta _2 x_2 + e\\ y &= \alpha_0 + \alpha _1 x_1 + e_1\\ y &= \gamma_0 + \gamma _1 x_2 + e_2 \end{aligned}

このとき、\beta_1 = \alpha_1と\beta_2 = \gamma_1でとなるのでしょうか？ やってみましょう。

まず、研究開発費の対数を説明変数として、売上高の対数を目的変数として回帰分析を行います。

res03 <- lm(log_売上高 ~ log_研究開発費, data = df)

これをmodelsummaryパッケージを使ってまとめてみましょう。

results <- list()
results[["単回帰1"]] <- res01 # 単回帰の結果
results[["単回帰2"]] <- res03 # 単回帰の結果
results[["重回帰"]] <- res02 # 単回帰の結果
modelsummary(results, # 結果を表にして出力
  stars = TRUE, # p値の有意性を星で表示する
  fmt = '%.3f', # 小数点以下3桁で表示する
  gof_omit = "AIC|BIC|Log.Lik.", # 不要な指標を除去
  output = 'html' # 出力形式をhtmlにする
  )
単回帰1  単回帰2 重回帰
(Intercept) 7.771*** 6.936*** 6.298***
(0.057) (0.064) (0.056)
log_広告宣伝費 0.551*** 0.324***
(0.007) (0.008)
log_研究開発費 0.631*** 0.402***
(0.008) (0.008)
Num.Obs. 4442 4442 4442
R2 0.567 0.595 0.713
R2 Adj. 0.567 0.595 0.713
F 5815.313 6515.603 5509.365
RMSE 1.08 1.04 0.88
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

\beta_1 = \alpha_1と\beta_2 = \gamma_1にはなっていないようです。 なぜでしょう？ それは、重回帰分析の係数は、他の変数を一定としたときのある説明変数が応答変数に与える影響の強さ、なので、単回帰ごとの結果と重回帰の結果は異なります。 そのことを確認するために、研究開発費と広告宣伝費の散布図を描いてみましょう。

ggplot(df) + aes(x = log_広告宣伝費, y = log_研究開発費) + geom_point() + mystyle

研究開発費と広告宣伝費の間には、やや強い相関があるようです。相関係数は、0.6306058となります。

この関係を前提として，単回帰分析の結果から重回帰分析の結果を導出してみます。 この方法を回帰解剖(regression anatomy)といいます。

まず，売上高の対数を広告宣伝費の対数で回帰した結果から回帰誤差を得ます。回帰残差を取り出すには、residuals()関数を使います。 この回帰誤差は，売上高の変動のうち，広告宣伝費とは関係ない部分といえます。

res <- lm(log_売上高 ~ log_広告宣伝費, data = df)
resid_sale <- residuals(res)

売上高の自然対数を広告宣伝費の自然対数で回帰した残差をresid_saleに代入しました。 次に，広告宣伝費を研究開発費で回帰します。

res_adv = lm(log_研究開発費 ~ log_広告宣伝費 , data = df)
coef(res_adv)
   (Intercept) log_広告宣伝費 
     3.6582599      0.5639725

この結果の残差は、研究開発費の対数の変動のうち、広告宣伝費の対数が説明できていない部分です。 つまり広告宣伝費のうち、研究開発費とは関係のない部分となります。 回帰残差を取り出すには、residuals()関数を使います。

resid_adv <- residuals(res_adv)
plot(resid_adv)

売上高の変動のうち広告宣伝費とは関係ない部分である残差resid_saleを，研究開発費の変動の内広告宣伝費とは関係のない部分である残差resid_advで回帰します。

res_ee <- lm(resid_sale ~ resid_adv)
round(coef(res_ee), digits = 3)
(Intercept)   resid_adv 
      0.000       0.402

このresid_advの回帰係数をみると，0.402となっており，前に計算した重回帰分析の研究開発費の自然対数の回帰係数と一致しております。

つまり，重回帰分析の係数とは，複数の説明変数の中で，他の説明変数とは関係の無い部分が，応答変数と他の説明変数との間の関係のない部分に与える影響の強さを表しているということです。

ここでいうと，重回帰分析における研究開発費の回帰係数の意味するところは，研究開発費の変動のうち広告宣伝費とは関係ない部分が，応答変数の変動のうち広告宣伝費とは関係の無い部分に与える影響，を意味しています。

10.6 決定係数

重回帰分析の結果を見ると、R^2という値が表示されています。 これは、決定係数と呼ばれる値で、回帰モデルの当てはまりの良さを表す指標の1つです。

決定係数について，以下の点くらいは覚えておくとよいでしょう。

  1. 決定係数は、0から1の値をとります。
  2. 決定係数は、説明変数の数が増えると必ず増加するので，重回帰分析では修正決定係数を使うことが多いです。
  3. 決定係数が1に近いほど、回帰モデルがデータによく当てはまっていることを表しますが，どれだけ高い数値だと良いかの目安はなく、研究分野やリサーチデザインによっても異なります。
  4. 決定係数は検定統計量ではないため，すべての回帰係数がゼロであるという帰無仮説を検定するF統計量でモデルの有意性を判断します。
  5. 決定係数以外にも，赤池情報量規準やベイズ情報量規準などの指標もありますが，今回は後者の2つは扱いません。