前回は線形回帰モデルを理解する上でのイロハの一つ、線形回帰モデルの目的変数(従属変数)が確率変数であることを説明しました。今回はより実用度の高い話である、「得られたデータを線形回帰モデルに当てはめて考える」ことの意義とその検討手法について説明します。

得られたデータは線形回帰モデルの実現値であるのか？

結論から言うと、得られたデータは線形回帰モデルの実現値ではありません。¹
では何のために「得られたデータが線形回帰モデルの実現値である」と仮定するのでしょうか？
それは、仮定することで、その仮定が妥当であるかを検討するためです。検討の結果、その仮定が実用上問題ないレベルで妥当だと判断できることはあります。検討の結果、その仮定が妥当でないのであれば、他のモデルに当てはめた場合の妥当性の検討に移り・・・検討を繰り返します。²
ここで統計学未学習者によくある勘違いは、「得られたデータが線形回帰モデルの実現値である」という仮定は統計学的手法を用いれば、「ある指標がある値を超えいればその仮定が正しい!」と明確に判断できると思ってたり、「その仮定が正しい確率は約70%である!」というように理解しやすい数値を計算できると思うことです。
できません。
では何ができるかを次に説明します。

雰囲気でしか判定できない

「得られたデータが線形回帰モデルの実現値である」という仮定の妥当性を検討する場合、主に以下のような検討方法があります。³一つずつ説明していきます。いずれの検討方法においても白黒つけられないことがわかるでしょう。

散布図(線形性の検討)

2変数間の関係であれば散布図で線形な関係っぽいかを検討できます。あくまで判断基準は分析者の感覚です。これとあわせて相関係数という指標もあるのですが、いくつ以上であれば線形な相関があるというような誰もが納得する閾値はありません。0.6以上であれば・・・0.8以上であれば・・・と業界ごとに共通の判断基準はあるかもしれませんが、そのような共通の判断基準はないことがほとんどです。分析の目的と照らし合わせて最終的には分析者の感覚で判断するしかないのです。⁴

# Juliaによる散布図作成の例
using Plots, CSV, DataFrames

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

# 説明変数として全384時点(1～384)を使う。
row = size(df_co2)[1]
mat_x = [r*c for r in 1:row, c in 1:1]
df_x = DataFrame(mat_x)

#世界の二酸化炭素濃度月平均の推移(いちおう散布図)をプロット
pyplot(
    titlefont=font("IPAexGothic"),
    guidefont=font("IPAexGothic"),
    tickfont=font("IPAexGothic"),
    legendfont=font("IPAexGothic")
)
plot(scatter(df_x.x1, df_co2.ave_ppm, legend=false), xlabel="時点(月数)", ylabel="CO2月平均濃度(ppm)", title="世界の二酸化炭素濃度月平均の推移(いちおう散布図)")

ヒストグラム(残差の正規性の検討1)

得られた目的変数値と、得られたデータ(目的変数値と説明変数値)から最小二乗推定した線形回帰式に説明変数を当てはめた推定値との差(残差の推定値)は、ヒストグラムにより正規分布っぽいかを目視で検討できます。判断基準は左右対称の山なりっぽいかです。明らかに違うとわかる場合を除いてこれも分析者の感覚に判断が委ねられます。⁵

# Juliaによるヒストグラム作成の例
using Plots, CSV, DataFrames, DataFramesMeta, GLM

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

# 説明変数として全384時点(0～383)を使う。
row = size(df_co2)[1]
mat_x = [r*c for r in 1:row, c in 1:1]
df_x = DataFrame(mat_x)

# データフレーム統一
df_co2 = @transform(df_co2, mat_x = df_x.x1)

# 線形回帰モデルを推定
model = lm(@formula(ave_ppm ~ mat_x), df_co2)
pred = predict(model)

# 残差の推定値算出
res = df_co2[:ave_ppm] - pred

# 残差の推定値のヒストグラム作成
pyplot(
    titlefont=font("IPAexGothic"),
    guidefont=font("IPAexGothic"),
    tickfont=font("IPAexGothic"),
    legendfont=font("IPAexGothic")
)
p = histogram(res
              , nbins=25
              , legend=false
              , xlabel="残差の推定値"
    　　　　　 , ylabel="件数"
    　　　　　 , title="残差の推定値のヒストグラム")

Q-Qプロットによる判定(残差の正規性の検討2)

正規Q-Qプロットは横軸に「標準化して昇順に並べた残差の推定値」、縦軸に「標準化して昇順に並べた残差の推定値の経験分布関数値に対応する標準正規分布の分位点」をプロットしたものです。⁶ 得られた残差の推定値が「ほぼ」全て、以下の例図における赤線($縦軸=横軸$)に「ほぼ」乗っていれば得られた残差の推定値が正規分布の実現値であると言えます。「ほぼ」がどれくらいであればよいのかの判断の按排は最終的には分析者の感覚に委ねられます。明確な基準はありません。

# JuliaによるQ-Qプロット作成の例
using Plots, CSV, DataFrames, DataFramesMeta, GLM, Statistics, CurveFit

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

# 説明変数として全384時点(0～383)を使う。
row = size(df_co2)[1]
mat_x = [r*c for r in 1:row, c in 1:1]
df_x = DataFrame(mat_x)

# データフレーム統一
df_co2 = @transform(df_co2, mat_x = df_x.x1)

# 線形回帰モデルを推定
model = lm(@formula(ave_ppm ~ mat_x), df_co2)
pred = predict(model)

# 残差の推定値算出
res = df_co2[:ave_ppm] - pred

# Q-Qプロット作成
N = length(res)
stdRes = sort(res .- mean(res) .* ones(N)) ./ (std(res)) # 標準化残差推定値を昇順に
list_expdist = [length(filter((x) -> x <= i, stdRes)) / N for i in stdRes] # 標準化残差推定値の経験分布値リスト
qx = quantile.(Normal(0,1), list_expdist) # 標準化残差推定値の経験分布値に対応する標準正規分布の分位点

x = stdRes
y = x
fittedFunction = curve_fit(Poly, x, y, 2) # Q-Qプロットの基準線定義

pyplot(
    titlefont=font("IPAexGothic"),
    guidefont=font("IPAexGothic"),
    tickfont=font("IPAexGothic"),
    legendfont=font("IPAexGothic")
)
plot(scatter(stdRes, qx, legend=false), xlabel="標準化して昇順に並べた残差推定値", ylabel="横軸に対応する標準正規分布の分位点", title="残差推定値のQ-Qプロット")
plot!(x, fittedFunction(x)) # Q-Qプロットの基準線

？？？(得られたデータが独立に得られているかの検討)

得られたデータがそれぞれ互いに独立な確率変数の実現値であるかどうかを検討する手段を私は知りません。ただ今回の場合は残差推定値それぞれが互いに同分布(正規分布)に従う確率変数の実現値かどうかを検討したいので、上記の正規性の検討方法で正規分布に従うと判断できたなら、同時に残差推定値が互いに独立に同じ正規分布に従う確率変数の実現値であると言えると思います。

残差推定値-最小二乗推定値のプロット(残差推定値それぞれが互いに独立に等分散分布からそれぞれ得られているかの検討)

以下は線形回帰モデルにおける残差の等分散の仮定を検討するためのプロットです。最小二乗推定値(横軸:目的変数の推定値)が大きくなるにつれて残差推定値(縦軸:最小二乗推定残差)のバラツキが大きくなっていたり小さくなっていれば等分散の仮定を満たしているとは言えませんし、何らかの関係性を持っている場合は残差推定値が互いに独立な分布からの実現値とは言えないでしょう。これにも明確な判断基準はありません。あくまで線形回帰モデルの推定値としての妥当性を検討していることにも留意する必要があります。

# Juliaによる残差推定値-最小二乗推定値プロット作成の例
using Plots, CSV, DataFrames, DataFramesMeta, GLM

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

# 説明変数として全384時点(0～383)を使う。
row = size(df_co2)[1]
mat_x = [r*c for r in 1:row, c in 1:1]
df_x = DataFrame(mat_x)

# データフレーム統一
df_co2 = @transform(df_co2, mat_x = df_x.x1)

# 線形回帰モデルを推定
model = lm(@formula(ave_ppm ~ mat_x), df_co2)
pred = predict(model)

# 残差の推定値算出
res = df_co2[:ave_ppm] - pred

#世界の二酸化炭素濃度月平均の推移をプロット
pyplot(
    titlefont=font("IPAexGothic"),
    guidefont=font("IPAexGothic"),
    tickfont=font("IPAexGothic"),
    legendfont=font("IPAexGothic")
)
plot(scatter(pred, res, legend=false), xlabel="最小二乗推定値", ylabel="最小二乗推定残差", title="残差推定値-最小二乗推定値プロット")

統計的仮説検定について

以上の検討手法として統計的仮説検定による検討の例を出しませんでした。理由は統計的仮説検定が難解だからです。統計的仮説検定の枠組みを理解するには相当の基礎知識を要します。下手に説明してしまうといらぬ誤解が生じてしまう恐れがあります。また、統計的仮説検定を用いても白黒結論はつけられれません。今後丁寧にわかりやすく一から説明できるようになることが今後の私の課題です。

雰囲気で正規分布っぽいとわかったら次にすべきこと

様々な手法の検討により得られたデータは線形回帰モデルからの実現値っぽいなと判断したとしましょう。次はその雰囲気を共通の認識にすることに尽力しましょう。雰囲気でそれっぽいというのは最初は分析者本人のみの認識です。ほかの人が見ればそうじゃない雰囲気を感じることもあるからです。多くの統計学未学習者はそうであるという絶対的な根拠を欲していて、それを絶対的に判定できる統計学的手法があると思っています。そんな手法がないことは前述したとおりです。なので、それをあらゆる手を使って説得するのです。一番効果的なのは既に出版されている書籍や投稿されている論文で行われている判定方法と同じ方法で判定したことをアピールすることです。そのような権威っぽいものに多くの人は弱いです。そこを突くのです。

まとめ

得られたデータは線形回帰モデルからの実現値ではありません。
得られたデータが線形回帰モデルの実現値であると仮定してその仮定が妥当かどうかを検討します。
様々な検討手法がありますが、どの手法を用いても白黒断定できません。
統計的仮説検定は難しいです。この手法を理解しようとする際は基礎からしっかり学習しましょう。
検討手法をクライアントに納得もらうよう尽力しましょう。

参考にした記事と書籍

https://qiita.com/piccolist/items/9632e3bafd46732ddff6
線形回帰分析 (統計ライブラリー) 蓑谷

次回

次回は、今回検討した線形回帰モデルの前提が妥当だと判断できたして、線形回帰モデルの推定係数値の妥当性を統計的仮説検定を用いて検討するうえで必要な基礎知識を説明します。

以上誤り等ご指摘いただけると助かります。

実用できるレベルで近似できることはあるかもしれませんが、完全一致はないと考えてよいでしょう。 ↩
線形回帰モデル以外のモデルを考えてもよいでしょう。 ↩
ここに並べた手法は私が知っている手法の一部を並べただけです。 ↩
無相関の仮説検定なるものがあるじゃないかと思われた方がいるかもしれませんが、これを使っても根本的な解決につながりません。今後別途説明します。あと、これ本当に実データか？と疑いたくなるような周期性を持った時系列データですね・・・ ↩
ks検定というある分布に従っていることを帰無仮説にして行う検定もありますが、この検定結果を残差の推定値が正規分布に従っているかの判断には使いづらいです。今回の場合は実用上問題ないレベルでおおよそ正規分布に従っていればよく、そもそもks検定の構造上、検定結果からは正規分布に従っているという断定はできません。 ↩
統計未学習者に対してはこれでは説明不足ですね。すみません。 ↩