雰囲気で分析している人のために(線形回帰モデル編1)


はじめに

得られたデータが線形回帰モデルの実現値であると仮定して最小二乗法で線形回帰モデルの係数を推定すると以下のような推定結果と検定結果が得られます。

(説明するための前振りです。例としてpythonのライブラリと落ちてたデータ使ってます。)

import pandas as pd
import statsmodels.api as sm

# 世界の二酸化炭素濃度月平均の推移
# (https://www.data.go.jp/data/dataset/mlit_20180523_0032)より
df_co2 = pd.read_csv('co2.csv')

# 年々世界の二酸化炭素濃度は上がってきてるのではないか? #

# 説明変数として全384時点(0~383)を使う。
df_co2['x'] = df_co2.index
X = df_co2.loc[:, ['x']]

# 目的変数として各月の二酸化炭素濃度(ppm)平均を使う。
Y = df_co2.loc[:, ['ave_ppm']]

# 最小二乗法で線形回帰モデルの係数を推定する。(時系列データでこんなことやるのは・・・)
model = sm.OLS(Y,sm.add_constant(X))
results = model.fit()
print(results.summary())
OLS Regression Results                            
==============================================================================
Dep. Variable:                ave_ppm   R-squared:                       0.983
Model:                            OLS   Adj. R-squared:                  0.983
Method:                 Least Squares   F-statistic:                 2.195e+04
Date:                Tue, 24 Dec 2019   Prob (F-statistic):               0.00
Time:                        00:01:54   Log-Likelihood:                -840.53
No. Observations:                 384   AIC:                             1685.
Df Residuals:                     382   BIC:                             1693.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        341.6819      0.221   1549.122      0.000     341.248     342.116
x              0.1477      0.001    148.154      0.000       0.146       0.150
==============================================================================
Omnibus:                       17.898   Durbin-Watson:                   0.198
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               10.180
Skew:                          -0.229   Prob(JB):                      0.00616
Kurtosis:                       2.347   Cond. No.                         442.
==============================================================================

近年分析と称して実際にやってみる人は増えたのではないでしょうか。
しかし、この推定結果の意味を理解している人はそれに対してかなり少ないと思ってます。(偏見)

なんか漠然と「p値が5%より小さいから有意な説明変数だ!仮説通り!」とか言って実は自分がやってることの意味がわかってない人の方が多いのではないでしょうか?
どのような場合でもデータさえあれば必ず統計的有意性を自動的に計算できるとか思ってる人も多そうです。
有意な結果を追い求めてるのに統計的有意性の意味がわかってないとかも・・・。
etc

線形回帰モデルは理解するのに相当な基礎知識を要する難解な統計モデルなのに、データと分析ツールがあれば分析結果は容易に出力できます。このギャップがこのような事態を招いていると思います。
ここで、

  • 「得られたデータが線形回帰モデルの実現値であると仮定して最小二乗法で線形回帰モデルの係数を推定すると以下のような推定結果と検定結果が得られます。」
  • 「線形回帰モデルは理解するのに相当な基礎知識を要する難解な統計モデルなのに、」

に対して、「え、そうなの?」、「どういうこと?」と思われた方がいるのではないでしょうか?
この記事はそういう人向けの記事です。
(もしくは分析ツールで出力してみたもののどう解釈していいのか?という人向けです。)

初回(線形回帰モデル編1)は、
「得られたデータが線形回帰モデルの実現値である。」
の意味を説明します。

「得られたデータが線形回帰モデルの実現値である。」とは?

まず線形回帰モデルが確率モデルであることを理解していない人が多いのではないでしょうか。線形回帰モデルは説明変数が1つの場合、

y_j = {\beta}_0 + {\beta}_{1}{x_{1j}}  + {u_j} \\
u_j  \sim N(0, \sigma^{2}), \quad i.i.d.\\
(j = 1, \cdots , n)\\

と表せます。1 
$y$が目的変数で$x_{1}$が説明変数の線形回帰モデルです。
この線形回帰モデルは、$y$と$x_{1}$の組のデータが$n$件得られたときに当てはめて考えることができるモデルの一つです。実際に得られたデータが必ずこのモデルで説明できるわけではないことにご注意ください。

「$u_j \sim N(0, \sigma^{2}), \quad i.i.d.$」は「$u_j$は各$j$間で互いに独立に平均$0$、分散$\sigma^{2}$、の正規分布に従う確率変数である。」ことを意味します。(線形回帰モデルでは$x_{1}$は確率変数ではありません。)
ここまでの線形回帰モデルの説明を聞いて確率変数とは?独立とは?などと疑問を持たれる場合は、まだ線形回帰モデルを理解するのに必要な基礎知識が足りていません。まずは用語の意味を理解しましょう。めんどい本題からそれるのでここでは説明しません。他のサイトや教科書を読んで理解願います。2

さて、以上の線形回帰モデルを見ればわかる通り、確率変数「$u_j \sim N(0, \sigma^{2}), \quad i.i.d.$」が含まれていることが線形回帰モデルが確率モデルである所以です。確率変数が含まれるモデルということです。この項$u_j$は誤差項と呼ばれます。

理解していない人によくあるのはこの誤差項が目に入らず、$y$が$\beta_i$と$x$の線形和でのみ表されていると勘違いしていることです。3ここだけを見て簡単なモデルだと思ってしまうパターンです。これは線形回帰モデルに従って得られたデータ値が確率変数の実現値であることをよく理解していないことからくる勘違いだと思います。

具体例として$y_j$が従う線形回帰モデルが$y_j = 1 + 2{x_{1j}} + {u_j}$であるとします。
このとき${x_{1j}} = 3$の場合$y_j$はいくつになるでしょうか?
ここでよく理解していない人は$y_j = 7$になると答えるでしょう。
もちろんこれは誤りです。こう答える人は$y_j$が確率変数であることを理解してないということです。
正しくは$y_j = 7 + {u_j}$なので$y_j$の値は${u_j}$の値によって決まります。
つまり$y_j$の値はサイコロの目の値のようにその時々によって変わるということです。
「サイコロをふって出た目の値はサイコロの実現値です。」のようにこの表現を使ってます。
実際にデータ値として得られた$y_j$の値は仮定している確率分布に従って得られた値というわけです。($y_j = 7 + {u_j}$の場合に$y_j$が従う確率分布は$N(7, \sigma^{2})$です。)

以上が「得られたデータが線形回帰モデルの実現値である。」ということの意味の説明でした。
ご質問、誤りのご指摘、などご意見いただけると幸いです。

次回

次回(線形回帰モデル編2)は、
『モデルを「仮定する」しかできない。』
ことを説明します。
よろしくお願いいたします。


  1. 誤差項が正規分布に従うモデルを仮定しないと、最小二乗推定された係数が正規分布に従わないし、残差平方和と$\sigma^2$の比がカイ二乗分布に従わないし、で一番最初にやってた$t$検定ができないしで・・・でも 

  2. 久保川統計(現代数理統計学の基礎)とかがわかりやすい教科書だと思います。測度論を理解しろとか言ってません。私も測度論は理解してません。ただ、確率分布の概念は理解しておきたいところです。 

  3. 私も昔そうでした。