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

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

演習問題2.1

次のコードでサンプルデータを作成する.

data = {
    "Treatment": [4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.98, 4.90, 5.75, 5.36, 3.48, 4.69, 4.44, 4.89, 4.71, 5.48, 4.32, 5.15, 6.34],
    "Control":   [4.17, 3.05, 5.18, 4.01, 6.11, 4.10, 5.17, 3.57, 5.33, 5.59, 4.66, 5.58, 3.66, 4.50, 3.90, 4.61, 5.62, 4.53, 6.05, 5.14]
}
df_2_8 = pd.DataFrame(data)

(a)

散布図を描くと次のような感じになり, 一見すると両群に違いはみられない.

plt.figure()
plt.scatter(range(len(df_2_8)), df_2_8["Treatment"], label="Treatment")
plt.scatter(range(len(df_2_8)), df_2_8["Control"], label="Control")
plt.xlabel("Index")
plt.ylabel("Dried weight")
plt.title("Scatter plot of Treatment and Control groups")
plt.legend()
plt.savefig("test.png")
plt.show()

統計量は以下のようになり, Treatmentの方が植物の乾重量が重い.

Statistic	Treatment	Control
count	20.000000	20.000000
mean	4.860000	4.726500
std	0.790922	0.863526
min	3.480000	3.050000
25%	4.387500	4.077500
50%	4.850000	4.635000
75%	5.390000	5.392500
max	6.340000	6.110000

(b)

帰無仮説を $H_0 : \mu_1 = \mu_2$ とし, 対立仮説を $H_0 : \mu_1 \neq \mu_2$ として設定する. 次のコードで対応のない $t$ 検定を行う. なお, 検定の際には, 等分散を仮定した.

from scipy.stats import t

# 各群における平均と分散
mean_treatment = df_2_8.mean().values[0]
mean_control = df_2_8.mean().values[1]
var_treatment = df_2_8.var(ddof = 1).values[0]
var_control = df_2_8.var(ddof = 1).values[1]

# poolした分散
pooled_variance = ((len(df_2_8) - 1) * var_treatment + (len(df_2_8) - 1) * var_control) / (len(df_2_8) + len(df_2_8) - 2)

# t値
t_value = (mean_treatment - mean_control)/np.sqrt(pooled_variance * (1/len(df_2_8) + 1/len(df_2_8)))

# 自由度
freedom = len(df_2_8) + len(df_2_8) - 2

# 95%信頼区間なので
alpha = 0.025
t_crit = t.ppf(1 - alpha, freedom)

# 群の平均の差
diff = mean_treatment - mean_control

# 標準誤差
se = np.sqrt(pooled_variance * (1/len(df_2_8) + 1/len(df_2_8)))

# 信頼区間
ci_low = diff - t_crit * se
ci_high = diff + t_crit * se

信頼区間は以下である. p値は出していないが, 信頼区間が0を含んでいるので, 帰無仮説は棄却されない.

指標	値
t値	0.5098476266071734
自由度	38
両群の平均の差	0.1335000000000015
95%信頼区間（下限）	-0.39657331363038195
95%信頼区間（上限）	0.663573313630385

(c)

群の数を $J$, サンプルサイズを $K$ とおく. 先に答えを言うと, 本問の状況下では, 最尤推定量と最小2乗推定量は一致する.

最尤推定量について

$\mu_1 \neq \mu_2$のとき, 対数尤度関数は

$$l(\mu_1, \mu_2;\ y_{11}, \cdots, y_{JK}) = -\sum_{j = 1}^J \sum_{k = 1}^K \left(\frac{1}{2} \log(2\pi \sigma^2) + \frac{1}{2\sigma^2}(y_{jk} - \mu_j)^2 \right)$$

となる. 対数尤度関数を $\mu_j$ で偏微分し, その値を0とおくことで

$$\hat{\mu}_{j,\ MLE} = \frac{1}{K}\sum_{k = 1}^K y_{jk}$$

を得る.

$\mu_1 = \mu_2$のとき, 同様の議論により

$$\hat{\mu}_{MLE} = \frac{1}{JK} \sum_{j = 1}^J \sum_{k = 1}^K y_{jk}$$

を得る.

最小二乗推定量について

$\mu_1 \neq \mu_2$のとき, $Y_{jk}$ を

$$Y_{jk} = \mu_j + \epsilon_{jk}, \, \epsilon_{jk} \sim \mathcal{N}(0, \sigma^2)$$

として定式化する. このとき, 最小二乗法により

$$\begin{aligned} S_1 := \sum_{j = 1}^J \sum_{k = 1}^K (y_{jk} - \mu_j)^2 \end{aligned}$$

を最小化すれば良い. $\mu_j$で偏微分することで

$$\begin{aligned} \frac{\partial S_1}{\partial \mu_j} = -2 \sum_{k = 1}^K (y_{jk} - \mu_j) \end{aligned}$$

を得る. よって

$$\hat{\mu}_{j,\ OLS} = \frac{1}{K} \sum_{k = 1}^K y_{jk}$$

を得る.

なお, $\mu_1 = \mu_2 = \mu$ のとき, つまり

$$Y_{jk} = \mu + \epsilon_{jk}, \, \epsilon_{jk} \sim \mathcal{N}(0, \sigma^2)$$

は同様の議論を行うことで

$$\begin{aligned} S_0 := \sum_{j = 1}^J \sum_{k = 1}^K (y_{jk} - \mu)^2 \end{aligned}$$

$$\hat{\mu}_{OLS} = \frac{1}{JK} \sum_{k = 1}^K y_{jk}$$

を得る.

最尤推定値および最小二乗推定量を求める

上記の方法で求めた推定量にサンプルデータの値を代入すれば良い.

記号	推定値
$\mu$	4.79325
$\mu_1$	4.86
$\mu_2$	4.7265

(d)

(c)の最小二乗推定量の導出結果をもとに, $J = 2$, $K = 20$ を代入すれば良い.

$$\begin{aligned} \hat{S}_0 &= \sum_{j = 1}^J \sum_{k = 1}^K (y_{jk} - \hat{\mu}_{OLS})^2 \\ &= \sum_{j = 1}^J \sum_{k = 1}^K (y_{jk} - \bar{Y})^2 \end{aligned}$$

$$\begin{aligned} \hat{S}_1 &= \sum_{j = 1}^J \sum_{k = 1}^K (y_{jk} - \hat{\mu}_{j,\ OLS})^2 \\ &= \sum_{j = 1}^J \sum_{k = 1}^K (y_{jk} - \bar{Y}_j)^2 \end{aligned}$$

を得る.

(e)

$H_1$ が真のとき

$\hat{S}_1$を変形すると

$$\begin{aligned} \hat{S}_1 &= \sum_{j = 1}^J \sum_{k = 1}^K (Y_{jk} - \bar{Y}_j)^2 \\ &= \sum_{j = 1}^J \sum_{k = 1}^K (Y_{jk} - \mu_j + \mu_j - \bar{Y}_j)^2 \\ &= \sum_{j = 1}^J \sum_{k = 1}^K (Y_{jk} - \mu_j)^2 -2 \sum_{j = 1}^J \sum_{k = 1}^K (Y_{jk} - \mu_j)(\mu_j - \bar{Y}_j) + \sum_{j = 1}^J \sum_{k = 1}^K (\bar{Y}_j - \mu_j)^2 \\ \end{aligned}$$

を得る. 第二項に着目すると

$$\begin{aligned} \sum_{j = 1}^J \sum_{k = 1}^K (Y_{jk} - \mu_j)(\mu_j - \bar{Y}_j) &= \sum_{j = 1}^J (\mu_j - \bar{Y}_j) \sum_{k = 1}^K (Y_{jk} - \mu_j) \\ &= \sum_{j = 1}^J (\mu_j - \bar{Y}_j) K(\bar{Y}_j - \mu_j) \\ \end{aligned}$$

となる. 元の式に代入すると

$$\begin{aligned} \hat{S}_1 = \sum_{j = 1}^J \sum_{k = 1}^K (Y_{jk} - \mu_j)^2 - \sum_{j = 1}^J \sum_{k = 1}^K (\bar{Y}_j - \mu_j)^2 \\ \end{aligned}$$

を得る. これを $\sigma^2$ で割り, $J = 2, K = 40$ を代入することで求める式

$$\begin{aligned} \frac{\hat{S}_1}{\sigma^2} = \sum_{j = 1}^2 \sum_{k = 1}^{20} \left(\frac{Y_{jk} - \mu_j}{\sigma} \right)^2 - 20 \sum_{j = 1}^2 \left(\frac{\bar{Y}_j - \mu_j}{\sigma} \right)^2 \\ \end{aligned}$$

が導出できる. さらに, $Y_{jk} \sim \mathcal(\mu_j, \sigma^2)$ であったので

$$\begin{aligned} \frac{Y_{jk} - \mu_j}{\sigma} &\sim \mathcal{N}(0, 1) \\ \frac{\bar{Y}_j - \mu_j}{\sqrt{\sigma^2/K}} &\sim \mathcal{N}(0, 1) \\ \end{aligned}$$

が成り立つ. 本文1.4.2の3より, 標準正規分布の二乗和はカイ二乗分布に従うことに注意すれば

$$\begin{aligned} \sum_{j = 1}^2 \sum_{k = 1}^{20} \left(\frac{Y_{jk} - \mu_j}{\sigma} \right)^2 &\sim \mathcal{\chi}^2(40) \\ 20 \sum_{j = 1}^2 \left(\frac{\bar{Y}_j - \mu_j}{\sigma} \right)^2 = \sum_{j = 1}^2 \left(\frac{\bar{Y}_j - \mu_j}{\sqrt{\sigma^2/20}} \right)^2 &\sim \mathcal{\chi}^2(2) \\ \end{aligned}$$

が成り立つ. 以上より, 本文1.5.6を適用すると

$$\frac{\hat{S}_1}{\sigma^2} \sim \mathcal{\chi}^2(40 - 2) = \mathcal{\chi}^2(38)$$

を得る.

$H_0$ が真のとき

$H_1$の時と同様にして示せるので省略する.

(f)

帰無仮説$H_0$が真のとき, (e)の結果および本文1.5の6を用いることで

$$\begin{aligned} \frac{\hat{S}_0 - \hat{S}_1}{\sigma^2} &\sim \chi^2(39 - 38) = \chi^2(1) \\ \frac{\hat{S}_1}{\sigma^2} &\sim \mathcal{\chi}^2(38) \end{aligned}$$

が成り立つ. 故に, F分布の定義により

$$\begin{aligned} \frac{\left( \frac{\hat{S}_0 - \hat{S}_1}{\sigma^2} \right)}{\left(\frac{\hat{S}_1}{38\sigma^2} \right)} = F \sim F(1, 38) \end{aligned}$$

となる.

$H_0$が真でない場合

$$\begin{aligned} \frac{\hat{S}_0 - \hat{S}_1}{\sigma^2} &\sim \chi^2(1, \lambda) \\ \frac{\hat{S}_1}{\sigma^2} &\sim \mathcal{\chi}^2(38) \end{aligned}$$

が成り立つ. 統計量Fの分子が非心カイ二乗分布に従い, 分母がカイ二乗分布に従う. 故に, 非心F分布の定義により, Fは非心F分布に従う.

(g)

以下のコードでF値を算出する.

S_0 = ((df_2_8 - df_2_8.sum().sum() / (df_2_8.shape[0]*df_2_8.shape[1]))**2).sum().sum()
S_1 = ((df_2_8 - df_2_8.mean())**2).sum().sum()
F = (S_0 - S_1) / (S_1/38)

結果は以下となった.

F = 0.2599446023569676

$F(1, 38)$に従うと仮定し, F検定を有意水準5%で実行すると, $F_{0.025}​(1,38) \approx 5.21 > F \approx 0.26$ なので, 帰無仮説$H_0$は棄却されない.

(h)

$F = t^2$ が成り立っている. 実際, (b)および(f)より

t = 0.5098476266071734
F = 0.2599446023569676

であるが, $t^2 = 0.25994460235696776$となり, 等しくなる. したがって, 両者の結論は同一(= 帰無仮説が棄却されない)である.

演習問題2.2

本問のデータを以下のようにデータフレームで作成する.

import numpy as np
import pandas as pd
from scipy.stats import norm
from scipy.stats import t

# 元のデータセット
df = pd.DataFrame({
    "id": list(range(1, 21)),
    "before": [
        100.8, 102.0, 105.9, 108.0, 92.0,
        116.7, 110.2, 135.0, 123.5, 95.0,
        105.0, 85.0, 107.2, 80.0, 115.1,
        103.5, 82.0, 101.5, 103.5, 93.0
    ],
    "after": [
        97.0, 107.5, 97.0, 108.0, 84.0,
        111.5, 102.5, 127.5, 118.5, 94.2,
        105.0, 82.4, 98.2, 83.6, 115.0,
        103.0, 80.0, 101.5, 102.6, 93.0
    ]
})

(a)

以下のコードで, 対応のないt検定を実施する. なお, 本問では, 2標本問題において, 標本平均の差の検定を実施しており, 母分散は未知だが等しいと仮定している.

$\hat{\sigma}_j^2$ を $j$ 回目の計測値を元に算出した不偏分散とする. プールされた分散を

$$s^2 = \frac{(m - 1)\hat \sigma^2_1+(n - 1) \hat\sigma^2_2}{m + n - 2}$$

とするとき, 帰無仮説 $H_0 : \mu_1 = \mu_2$ のもとで, 検定統計量 $T$ は

$$T = \frac{\bar{Y}_{1k} - \bar{Y}_{2k}}{\displaystyle \sqrt{s^2\left(\frac{1}{m} + \frac{1}{n} \right)}}$$

であることに注意する.

# 各群における不偏分散
unbiased_var_before = df["before"].var(ddof = 1)
unbiased_var_after = df["after"].var(ddof = 1)

# 自由度
dof = (len(df) - 1) + (len(df) - 1)

# poolされた分散
pooled_var = ((len(df) - 1) * unbiased_var_before + (len(df) - 1) * unbiased_var_after) / dof

# 対応のないt検定統計量
mean_diff = df["before"].mean() - df["after"].mean()
t_test_stats = (mean_diff) / (np.sqrt(pooled_var * (1/len(df) + 1/len(df))))

# 上側/下側2.5%点
t_upper = t.ppf(0.975, dof)
t_lower = t.ppf(0.025, dof)

# 信頼限界
lower_conf_limit = float(mean_diff - t_upper * np.sqrt(pooled_var *  (1/ len(df) + 1/ len(df))))
upper_conf_limit = float(mean_diff + t_upper * np.sqrt(pooled_var * (1/ len(df) + 1/ len(df))))

# p-value
p_value = 2 * (1 - t.cdf(abs(t_test_stats), dof))

print(f"自由度 : {dof}")
print(f"t検定統計量 : {t_test_stats}")
print(f"95%信頼区間 : {[lower_conf_limit, upper_conf_limit]}")
print(f"p値 : {p_value}")
print(f"上側2.5%点 : {t_upper}")
print(f"下側2.5%点 : {t_lower}")

以下の結果を得る. 95%信頼区間が0を含んでしまっているので, 帰無仮説は棄却できない.

自由度 : 38
t検定統計量 : 0.6430870454738202
95%信頼区間 : [-5.68127962456002, 10.97127962456004]
p値 : 0.524029521709561
上側2.5%点 : 2.024394163911969
下側2.5%点 : -2.0243941639119694

(b)

$D_k = Y_{1k} - Y_{2k}$ として新たに $D_k$ を導入すると, $D_k \sim \mathcal{N}(\mu_D, \sigma_D^2)$ として, (a) の仮説検定は次のように定式化できる.

$$\begin{aligned} H_0 &: \mu_D = 0 \\ H_1 &: \mu_D \neq 0 \end{aligned}$$

D_k = df["before"] - df["after"]

n = len(D_k)
dof = n - 1

# 標本平均と標準偏差
mean_D = D_k.mean()
s_D = D_k.std(ddof=1)

# t統計量
t_stat = mean_D / (s_D / np.sqrt(n))

# 臨界値
t_upper = t.ppf(0.975, dof)

# 信頼区間
lower_conf = float(mean_D - t_upper * (s_D / np.sqrt(n)))
upper_conf = float(mean_D + t_upper * (s_D / np.sqrt(n)))

# p値（両側）
p_value = 2 * (1 - t.cdf(abs(t_stat), dof))

print(f"自由度 : {dof}")
print(f"t検定統計量 : {t_stat}")
print(f"95%信頼区間 : {[lower_conf, upper_conf]}")
print(f"p値 : {p_value}")
print(f"上側2.5%点 : {t_upper}")
print(f"下側2.5%点 : {t_lower}")

出力結果は以下である. 信頼区間が0を含んでおらず, t検定統計量が臨界値よりも大きい値をとっていることから, 有意水準5%で帰無仮説は棄却される. つまり, 減量プログラムの参加前後で体重の増減はなかったことが示唆される.^[1]

自由度 : 19
t検定統計量 : 2.873403522553378
95%信頼区間 : [0.7183480068297599, 4.571651993170241]
p値 : 0.009730462990029443
上側2.5%点 : 2.093024054408263
下側2.5%点 : -2.0243941639119694

(c)

(a) と (b) で結論は異なる.

(d)

(a)においては, $Y_{jk} \sim \mathcal{N}(\mu_j, \sigma^2)$ であり, 体重の分布は計測時点に依存すると仮定していた. 一方, (b) については, $D_k \sim \mathcal{N}(\mu_D, \sigma^2_D)$ であり, 体重は参加前後の時点に依存していないと仮定されている. 更に, (b) では, 同一被験者のプログラム参加前後の体重変化($D_k$)を用いて仮説検定されているが, 通常体重変化には個人差があると思われるので, 同一個人の差を計測した方が減量維持できたか判断する上では, 効果的と考えられる. 従って, (b) の方が適切な手法であるとと考えられる.

演習問題2.3

$\hat{S}_1$ の分解について

$$\begin{aligned} Y_{jk} - a_j - b_j x_{jk} &= Y_{jk} - \alpha_j - \beta_j x_{jk} \\ &\quad + (\beta_j - b_j) (x_{jk} - \bar{x}_j) \\ &\quad - \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right) \\ &\quad + \left(\bar{Y}_j - a_j -b_j \bar{x}_j \right) \\ &= \left(Y_{jk} - \alpha_j - \beta_j x_{jk} \right) \\ &\quad + \left(\beta_j - b_j \right) \left(x_{jk} - \bar{x}_j \right) \\ &\quad - \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right) \\ \end{aligned}$$

と変形できる. ただし, $a_j, b_j$ が最尤推定量であることにより, $\bar{Y}_j - a_j - b_j \bar{x}_j = 0$ となることに注意する. 両辺を2乗することにより

$$\begin{align} \left( Y_{jk} - a_j - b_j x_{jk} \right)^2 &= \left( Y_{jk} - \alpha_j - \beta_j x_{jk} \right)^2 \\ &\quad + (\beta_j - b_j)^2 (x_{jk} - \bar{x}_j)^2 \\ &\quad + \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right)^2 \\ &\quad +2\left( Y_{jk} - \alpha_j - \beta_j x_{jk} \right) \left(\beta_j - b_j \right) \left(x_{jk} - \bar{x}_j \right) \\ &\quad -2 \left(\beta_j - b_j \right) \left(x_{jk} - \bar{x}_j \right) \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right) \\ &\quad -2 \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right) \left( Y_{jk} - \alpha_j - \beta_j x_{jk} \right) \end{align}$$

となる.

(3)式と(6)式の和

(3)式と(6)式の和は

$$\begin{aligned} (3) + (6) &= \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right) \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j -2 (Y_{jk} -\alpha_j -\beta_j x_{jk}) \right) \\ &= \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right) \left((\bar{Y}_j -Y_{jk}) - \beta_j (\bar{x}_j - x_{jk}) \right) \\ &- \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right) \left(Y_{jk} -\alpha_j -\beta_j x_{jk}) \right) \\ \end{aligned}$$

となる. 上式の両辺に対し $k$ で和をとると, $\sum_{k = 1}^K (Y_{jk} - \bar{Y}_j) = 0,\, \sum_{k = 1}^K (X_{jk} - \bar{X}_j) = 0$ であることに注意すれば 2 番目の等式の右辺第1項は 0 となる. よって

$$\begin{aligned} \sum_{j = 1}^J \sum_{k = 1}^K \left( (3) + (6) \right) &= - \sum_{j = 1}^J \sum_{k = 1}^K \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right) \left(Y_{jk} -\alpha_j -\beta_j x_{jk}) \right) \\ &= -K \sum_{j = 1}^J \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right)^2 \end{aligned}$$

(4)式と（５）式の和

(4)式と(5)式の和は

$$\begin{aligned} (4) + (5) &= 2\left(\beta_j - b_j \right) \left(x_{jk} - \bar{x}_j \right) \left (Y_{jk} - \bar{Y}_j - \beta_j (x_{jk} - \bar{x}_j) \right) \\ &= - 2\left(\beta_j - b_j \right) \left(x_{jk} - \bar{x}_j \right) \left (Y_{jk} - (a_j + b_j \bar{x}_j ) - \beta_j (x_{jk} - \bar{x}_j) \right) \\ &= 2\left(\beta_j - b_j \right) \left(x_{jk} - \bar{x}_j \right) \left (Y_{jk} - (a_j + b_j \bar{x}_j ) - \beta_j (x_{jk} - \bar{x}_j) \right) \\ &= 2\left(\beta_j - b_j \right) \left(x_{jk} - \bar{x}_j \right) \left ((b_j - \beta_j)(x_{jk} - \bar{x}_j) + Y_{jk} - a_j - b_j x_{jk} \right) \\ &= -2 \left(b_j - \beta_j \right)^2 \left(x_{jk} - \bar{x}_j \right)^2 + 2\left(\beta_j - b_j \right) \left(x_{jk} - \bar{x}_j \right)(Y_{jk} - a_j - b_j x_{jk}) \end{aligned}$$

となる. $a, b$ は最尤推定量なので

$$\sum_{k = 1}^K (Y_{jk} - a_j - b_j x_{jk}) = 0 , \sum_{k = 1}^K x_{jk}(Y_{jk} - a_j - b_j x_{jk}) = 0$$

を満たす. 故に, $j, k$ について和をとると

$$\begin{aligned} \sum_{j = 1}^J \sum_{k = 1}^K \left( (4) + (5) \right) &= -2 \sum_{j = 1}^J \sum_{k = 1}^K \left(b_j - \beta_j \right)^2 \left(x_{jk} - \bar{x}_j \right)^2 \\ &= -2 \sum_{j = 1}^J \left(b_j - \beta_j \right)^2 \left( \sum_{k = 1}^K x_{jk}^2 - K\bar{x}_j^2 \right) \\ \end{aligned}$$

を得る. なお, (2)式はここで行った変形と同様の処理により, (2)式 = (-1/2)((4)式 + (5)式)であることがわかる.

(1)式 ~ (6)式の合計

以上より, (1) ~ (6)を足し合わせ, $j = 1,...,J, k = 1, ..., K$ で足し合わせると

$$\begin{aligned} \hat{S}_1 &= \sum_{j = 1}^J \sum_{k = 1}^K \left( Y_{jk} - \alpha_j - \beta_j x_{jk} \right)^2 \\ &- K \sum_{j = 1}^J \left( \bar{Y}_j - \alpha_j - \beta_j \bar{x}_j \right) \\ &- \sum_{j = 1}^J \left(b_j - \beta_j \right)^2 \left( \sum_{k = 1}^K x_{jk}^2 - K\bar{x}_j^2 \right) \end{aligned}$$

を得る.

$Y_{jk}, \bar{Y}_j, b_j$ の従う分布について

$Y_{jk} \sim N(\alpha_j + \beta_j x_{jk},\sigma^2)$ の証明

モデル(2.7)より, $Y_{jk} \sim N(\mu_{jk}, \sigma^2)$ であり,

$$\begin{aligned} E[Y_{jk}] = \mu_{jk} = \alpha_j + \beta_j x_{jk} \end{aligned}$$

なので ^[2], $Y_{jk} \sim N(\alpha_j + \beta_j x_{jk}, \sigma^2)$ となる.

$\bar{Y}_j \sim N(\alpha_j + \beta_j \bar{x}_j, \sigma^2)$ の証明

正規分布の再生性および上記の証明により, $\bar{Y}_j \sim N(\alpha_j + \beta_j \bar{x}_j, \sigma^2)$ となる.

$b_j \sim N(\beta_j, \sigma^2/(\sum_{k = 1}^K x_{jk}^2 - K \bar{x}_j^2))$ の証明

愚直に計算する. モデル(2.7)より, 誤差項を含めて記載すると

$$Y_{jk} = \alpha_j + \beta_j x_{jk} + \epsilon_{jk}, \epsilon_{jk} \sim N(0, \sigma^2)$$

として表せる.^[3] ここで, $b_j$を展開し, 上記の誤差項を用いると

$$\begin{aligned} b_j &= \frac{K \sum_{k = 1}^K x_{jk}y_{jk} -(\sum_k x_{jk})(\sum_k x_{jk})}{K \sum_{k} x_{jk}^2 - (\sum_k x_{jk})^2} \\ &= \frac{\sum_{k = 1}^K (x_{jk} - \bar{x}_j) (y_{jk} - \bar{y}_j)} {\sum_{k} (x_{jk} - \bar{x}_j)^2} \\ &= \frac{\sum_{k = 1}^K (x_{jk} - \bar{x}_j) (\beta_j(x_{jk} - x_{j})+\epsilon_{jk} - \bar{\epsilon}_j)} {\sum_{k} (x_{jk} - \bar{x}_j)^2} \\ &= \beta_j + \frac{\sum_{k = 1}^K (x_{jk} - \bar{x}_j)(\epsilon_{jk} - \bar{\epsilon}_j)} {\sum_{k} (x_{jk} - \bar{x}_j)^2} \\ &= \beta_j + \frac{\sum_{k = 1}^K (x_{jk} - \bar{x}_j)\epsilon_{jk}} {\sum_{k = 1}^K (x_{jk} - \bar{x}_j)^2} \\ \end{aligned}$$

を得る. また, $\epsilon_{jk}$ は独立に $N(0, \sigma^2)$ に従うので $E[b_j] = \beta_j$ が成立する. 故に

$$\begin{aligned} Var(b_j) &= \frac{\sum_{k = 1}^K (x_{jk} - \bar{x}_j)^2}{\left( \sum_{k = 1}^K(x_{jk} - \bar{x}_j)^2 \right)^2}Var(\epsilon_{jk}) &= \frac{\sigma^2}{\sum_{k = 1}^K(x_{jk} - \bar{x}_j)^2} \end{aligned}$$

である. 確率変数の和の分散を計算する際には, 各 $\epsilon_{jk}$ が独立であることを用いた.

演習問題2.4

$$\begin{aligned} E[Y] &= \log( \beta_0 + \beta_1 x + \beta_2 x^2) \Rightarrow e^{E[Y]} = \beta_0 + \beta_1 x + \beta_2 x^2 \end{aligned}$$

より, $g(X) = e^X$ である. (2.13)の形で書くと

$$\exp \left[ \begin{array}{ccc} E[Y_1] \\ \vdots \\ E[Y_n] \end{array} \right] = \left[ \begin{array}{ccc} 1 & x_{11} & x_{11}^2 \\ \vdots & \vdots & \vdots \\ 1 & x_{n1} & x_{n1}^2 \end{array} \right] \left[ \begin{array}{ccc} \beta_0 \\ \beta_1 \\ \beta_2 \end{array} \right]$$

と表せるので, 与えられたデータを代入することで

$$X = \left[ \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 1.2 & 1.44 \\ 1 & 1.4 & 1.96 \\ 1 & 1.6 & 2.56 \\ 1 & 1.8 & 3.24 \\ 1 & 2 & 4 \end{array} \right], \quad y = \left[ \begin{array}{ccc} 3.15 \\ 4.85 \\ 6.50 \\ 7.20 \\ 6.25 \\ 16.50 \end{array} \right], \quad \beta = \left[ \begin{array}{ccc} \beta_0 \\ \beta_1 \\ \beta_2 \end{array} \right]$$

である.

演習問題2.5

零和制約より, $\alpha_2 = - \alpha_1,\, \beta_3 = - \beta_1 - \beta_2$ なので

$$\begin{aligned} E[Y_{11}] &= \mu + \alpha_1 + \beta_1 \\ E[Y_{12}] &= \mu + \alpha_1 + \beta_2 \\ E[Y_{13}] &= \mu + \alpha_1 - \beta_1 - \beta_2 \\ E[Y_{21}] &= \mu - \alpha_1 + \beta_1 \\ E[Y_{22}] &= \mu - \alpha_1 + \beta_2 \\ E[Y_{23}] &= \mu - \alpha_1 - \beta_1 - \beta_2 \\ \end{aligned}$$

となる. 行列表示で書くと

$$\begin{aligned} \begin{pmatrix} E[Y_{11}] \\ E[Y_{12}] \\ E[Y_{13}] \\ E[Y_{21}] \\ E[Y_{22}] \\ E[Y_{23}] \\ \end{pmatrix} &= \begin{pmatrix} \mu + \alpha_1 + \beta_1 \\ \mu + \alpha_1 + \beta_2 \\ \mu + \alpha_1 - \beta_1 - \beta_2 \\ \mu - \alpha_1 + \beta_1 \\ \mu - \alpha_1 + \beta_2 \\ \mu - \alpha_1 - \beta_1 - \beta_2 \\ \end{pmatrix} \\ &= \begin{pmatrix} 1 & 1 & 1 & 0 \\ 1 & 1 & 0 & 1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & 1 & 0 \\ 1 & -1 & 0 & 1 \\ 1 & -1 & -1 & -1 \\ \end{pmatrix} \begin{pmatrix} \mu \\ \alpha_1 \\ \beta_1 \\ \beta_2 \\ \end{pmatrix} \\ \end{aligned}$$

となる.

脚注

1. 解答例を見ると, F検定統計量が算出されているが, ここまでの文脈で, 突然出さなくてもt検定統計量で統一した方がわかりやすいと思う.

2. 本文p27のモデル(2.7)より引用

3. 本書には誤差項を記載せずに線形モデルを論じており, その点でこの定義はやや天下りかもしれない. ただ誤差項を $\epsilon_{jk} := Y_{jk} - (\alpha_j + \beta_j x_{jk})$ として定義することで, Yが独立に正規分布に従うという仮定から, 誤差に対してもiidの仮定が自然に継承される.