一般化線形モデル入門 第6章 解答例

https://zenn.dev/honwakasan/articles/b9ae08fd8a694e の第6章の解答例です.

演習問題6.1

(a)

df = pd.DataFrame(
    {
        "時期":["1936-39", "1946-49", "1956-59", "1966-69", "1976-79", "1986-89"],
        "精糖の消費量":[32, 31.2, 27, 21, 14.9, 8.8],
        "加工食品に含まれる砂糖の消費量":[16.3, 23.1, 23.6, 27.7, 34.6, 33.9]
    }
)

plt.title("精糖及び加工食品に含まれる砂糖の消費量")
plt.scatter(df["時期"], df["精糖の消費量"], label="精糖")
plt.scatter(df["時期"], df["加工食品に含まれる砂糖の消費量"], label = "砂糖")
plt.legend()
plt.show()

により砂糖の消費量と期間の散布図は以下になる. 期間を経るに従い, 加工食品に含まれる砂糖の消費量は上昇するが, 精糖の消費量は減少していることがわかる.

df_concat = pd.concat(
    [
        df,pd.get_dummies(df["時期"]).astype(int)
    ],
    axis = 1
)
df_concat = sm.add_constant(df_concat)
df_concat["期間"] = list(range(1, len(df_concat) + 1))

食品に含まれる砂糖の消費量の場合

以下のコードで食品に含まれる砂糖の消費量と機関との関係を回帰モデルにより推定する.

from scipy.stats import t

alpha = 0.05
n, p = X.shape
t_crit = t.ppf(1 - alpha/2, df = n - p)

X = np.stack([np.ones(6), df_concat["期間"].values], axis = 1)
y = df_concat["加工食品に含まれる砂糖の消費量"].values

# 係数
beta = np.linalg.inv(X.T@X)@X.T@y

# 標準誤差の計算
sigma_hat_sq = ((y - X@beta).T@(y - X@beta)) / (n - p)
std_error = np.diag(np.sqrt(np.linalg.inv(X.T@X)*sigma_hat_sq))

# 表に出力
table = pd.DataFrame(beta, index = ["const", "係数"], columns = ["推定値"])
table["標準誤差"] = std_error
table["95%信頼区間下限"] = beta - t_crit * table["標準誤差"]
table["95%信頼区間上限"] = beta + t_crit * table["標準誤差"]
print(table)

以下を得る. 期間を経るごとに食品に含まれる砂糖の消費量が増加していることがわかる. また, 信頼区間が0を含んでいないことから, 5%で有意であることもわかる.

	推定値	標準誤差	95%信頼区間下限	95%信頼区間上限
const	13.873333	1.910613	8.568622	19.178045
係数	3.617143	0.490600	2.255019	4.979267

精糖の消費量

また, 加工食品に含まれる砂糖の消費量ではなく, 精糖の場合の結果も同様に係数を推定できる.

from scipy.stats import t

alpha = 0.05
n, p = X.shape
t_crit = t.ppf(1 - alpha/2, df = n - p)

X = np.stack([np.ones(6), df_concat["期間"].values], axis = 1)
y = df_concat["精糖の消費量"].values

# 係数
beta = np.linalg.inv(X.T@X)@X.T@y

# 標準誤差の計算
sigma_hat_sq = ((y - X@beta).T@(y - X@beta)) / (n - p)
std_error = np.diag(np.sqrt(np.linalg.inv(X.T@X)*sigma_hat_sq))

# 表に出力
table = pd.DataFrame(beta, index = ["const", "係数"], columns = ["推定値"])
table["標準誤差"] = std_error
table["95%信頼区間下限"] = beta - t_crit * table["標準誤差"]
table["95%信頼区間上限"] = beta + t_crit * table["標準誤差"]
print(table)

以下を得る. 期間を経るごとに精糖の消費量が減少していることがわかる. また, 信頼区間が0を含んでいないことから, 5%で有意であることもわかる.

	推定値	標準誤差	95%信頼区間下限	95%信頼区間上限
const	39.573333	1.899239	34.300201	44.846466
係数	-4.882857	0.487679	-6.236872	-3.528842

(b)

以下のコードにより総消費量を計算し, 散布図を作成した.

df_concat["総消費量"] = df_concat["精糖の消費量"] + df_concat["加工食品に含まれる砂糖の消費量"]
plt.title("総消費量")
plt.scatter(df_concat["時期"], df_concat["総消費量"], label="総消費量")
plt.legend()
plt.show()

得られた結果は以下の図である. 変動が大きく, 1946年から減少傾向であることがわかる.

また, 総消費量が各期間において一定（= 回帰係数の値が0）という仮説を検定する. つまり, 帰無仮説 ${H_0}$ を「期間の回帰係数 = 0」として設定する.

alpha = 0.05
n, p = X.shape
t_crit = t.ppf(1 - alpha/2, df = n - p)

X = np.stack([np.ones(6), df_concat["期間"].values], axis = 1)
y = df_concat["総消費量"].values

# 係数
beta = np.linalg.inv(X.T@X)@X.T@y

# 標準誤差の計算
sigma_hat_sq = ((y - X@beta).T@(y - X@beta)) / (n - p)
std_error = np.diag(np.sqrt(np.linalg.inv(X.T@X)*sigma_hat_sq))

# t値
t_values = beta / std_error

# 両側p値
p_values = 2 * (1 - t.cdf(np.abs(t_values), df=n - p))

# 表に出力
table = pd.DataFrame(beta, index = ["const", "係数"], columns = ["推定値"])
table["標準誤差"] = std_error
table["95%信頼区間下限"] = beta - t_crit * table["標準誤差"]
table["95%信頼区間上限"] = beta + t_crit * table["標準誤差"]
table["p値"] = p_values
table

により以下を得る. 係数値 = 11.0681 であり, p-value = 0.182 > 0.05 である^[1]ことから, 帰無仮説は棄却されない.^[2] よって, 総消費量が各期間において一定かどうかをデータから判断することを留保する.^[3]

	推定値	標準誤差	95%信頼区間下限	95%信頼区間上限	p値
const	53.446667	3.057030	44.958991	61.934342	0.000063
係数	-1.265714	0.784973	-3.445148	0.913720	0.182163

脚注

1. 本書には, p値の正確な定義が記載されていないため計算するのが難しい. p値の定義については, https://www.kyoritsu-pub.co.jp/book/b10003681.html などを参照

2. 帰無仮説を検証するのにp値を計算する必要は特にないが, なんとなく計算してみた. 実際には, t値とt分布のパーセンタイル点が計算できれば事足りる.

3. 有意水準は前段の値である5%に設定している. また, 実際の実証分析では, 5%有意かどうかだけで仮説を検討するのは危ういがここでは深入りしない.