①まず,応答変数が従う確率分布を指定し,モデルを特定しよう
②自分が興味のある問題(説明変数Xは応答変数Yに効果をもつかなど)に基づき,帰無仮説と対立仮説を立て, それぞれの仮説のもとでのモデルを構築しよう
③モデルのパラメータを推定しよう(ML,LSM)
④推定値と実測値の差:残差から導いた統計量を特定の確率分布に近似し, 「帰無仮説におけるモデル」のあてはまりの程度に対する「対立仮説におけるモデル」のあてはまりの程度の差が統計的に有意な差であるかを検定しよう
(ここでの仮説検定は逸脱度(対数尤度比統計量)を用いたいわゆる尤度比検定ではないと思います。 また,上記の内容のうち仮説検定は第5章で詳細に説明されています)
データ:オーストラリアの農村部と都市部にすむ女性の慢性病状の数(年齢層,社会的地位,一般開業医への来院回数は統一)
sick<- read.csv("2.2_sick.csv")
summary(sick)
## num_sick group
## Min. :0.000 rural:23
## 1st Qu.:0.000 urban:26
## Median :1.000
## Mean :1.184
## 3rd Qu.:2.000
## Max. :4.000
head(sick)
## num_sick group
## 1 0 urban
## 2 1 urban
## 3 1 urban
## 4 0 urban
## 5 2 urban
## 6 3 urban
疑問:都市群と農村群の間で診察の必要性(慢性病状数により評価)は等しいか否か?
問題の定式化: j群のk番目の女性の病状数を\(Y_{jk}\)とする
ただし
j=1(都市群), 2(農村群)
k=1,…,\(K_{j}\)であり, \(K_{1}\)=26, \(K_{2}\)=23
\(Y_{jk}\)はポアソン分布(群内標本平均=標本分散,計数データが従う分布)に従うと仮定
帰無仮説\(H^0\):\(\theta_{1}\)=\(\theta_{2}\)=\(\theta\)が真だと仮定したモデル:
\[E(Y_{jk})=\theta; Y_{jk} ~ Poisson(\theta)\]GLM的表記なら sick~1
GLM的表記なら sick~group
モデルのあてはまりを比較
\(H^0\)が真のとき,\(Y_{jk}\)の対数尤度関数は:
\[l_{0}=l(\theta;y)=\sum_{j=1}^{J}\sum_{k=1}^{K_{j}}(y_{jk}log\theta-\theta-logy_{jk}!)\]
∴\(最尤推定値:\hat\theta\)=1.184
∴\(対数尤度関数最大値:\hat{l_{0}}\)=-68.3868
\(H^1\)が真のとき,\(Y_{jk}\)の対数尤度関数は: \[l_{1}=l(\theta_{1}, \theta_{2}, y)=\sum_{k=1}^{K_{1}}(y_{1k}log\theta_{1}-\theta_{1}-logy_{1k}!)+\sum_{k=1}^{K_{2}}(y_{2k}log\theta_{2}-\theta_{2}-logy_{2k}!)\]
最尤推定値\(\hat{\theta}\)は:∴\(\hat\theta_{1}\)=1.423
∴\(\hat\theta_{1}\)=0.913
∴\(\hat{l_{1}}\)=-67.0230
\(\hat{l_{0}}\)<\(\hat{l_{1}}\)となっているが,これは\(l_{1}\)が\(l_{0}\)よりもパラメータを1つ多くもつため.その差が有意かを知るには対数尤度関数の標本分布を知る必要あり(第45章)
残差による比較
\(Y\)~\(Poisson(\theta)\)のとき, E(Y)=var(Y)=\(\theta\)となり,
E(Y)の推定値\(\hat\theta\)を当てはめ値(fitted value),
Y-\(\hat\theta\)を残差(residual)とよぶ
ポアソン分布における近似的な標準化残差は:
\[\gamma=\frac{Y-\hat\theta}{\sqrt{\hat\theta}}\]
*通常の適合度カイ2乗統計量と同じ
\(H^0\)が真と仮定したモデルの下では:
\(\sum r_i^2\)=46.759
\(H^1\)が真と仮定したモデルの下では:
\(\sum r_i^2\)=43.659
いずれの統計量も,各モデル下の自由度における\(\chi^2\)分布のもとで生じうる確率が極めて大きい
∴\(H^1\)が真と仮定したモデルは,\(H^0\)が真と仮定したモデルと比べてそれほどよくデータを表現するわけではない, すなわちデータは\(H^1\):\(\theta_{1}\)≠\(\theta_{2}\)を指示する根拠を与える
fit0 <- glm(num_sick~1,data=sick,family="poisson")
fit1 <- glm(num_sick~group,data=sick,family="poisson")
anova(fit0,fit1,test="Chisq")
## Analysis of Deviance Table
##
## Model 1: num_sick ~ 1
## Model 2: num_sick ~ group
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 48 56.033
## 2 47 53.306 1 2.7277 0.09862 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ちなみにRで「カイ2乗近似による仮説検定(2つのモデルの比較)」を実行すると, p=0.098となり有意ではない
ただしDobson本の例では残差(標準化残差)を, カイ2乗近似の尤度比検定では逸脱度(deviance)を用いている
データ:ある病院で生まれた12人ずつの男児と女児の出生時体重と妊娠期間(週)
d<-read.csv("2.2_weight_duration.csv")
summary(d)
## duration weight sex
## Min. :35.00 Min. :2412 female:12
## 1st Qu.:37.00 1st Qu.:2785 male :12
## Median :38.50 Median :2952
## Mean :38.54 Mean :2968
## 3rd Qu.:40.00 3rd Qu.:3184
## Max. :42.00 Max. :3473
head(d)
## duration weight sex
## 1 40 2968 male
## 2 38 2795 male
## 3 40 3163 male
## 4 35 2925 male
## 5 36 2625 male
## 6 37 2847 male
疑問:妊娠期間の増加に伴う出生時体重の増加は男児と女児の間で同じか否か?
問題の定式化:
第j群のk番目の用事の出生時体重を表す確率変数を\(Y_{jk}\)とする。
j=1(男児), 2(女児), k=1,…,12
\(Y_{jk}\)は正規分布(平均\(\mu_{jk}=E(Y_{jk})\), 分散\(\sigma^2\))に従うと仮定
\(\alpha\):切片 \(\beta\):傾き
帰無仮説\(H^0\):\(\beta_{1}=\beta_{2}=\beta\)
対立仮説\(H^1\):\(\beta_{1}≠\beta_{2}\)
\(H^{0}\)が真であると仮定:
\[E(Y_{jk})=\mu_{jk}=\alpha_{j}+\beta x_{jk}; Y_{jk}~N(\mu_{jk}),\sigma^2\]
\(H^{1}\)が真であると仮定: \[E(Y_{jk})=\mu_{jk}=\alpha_{j}+\beta_{j}x_{jk}; Y_{jk}~N(\mu_{jk}),\sigma^2\]
\(Y_{jk}\)の確率密度関数:
\[f(y_{jk}; \mu_{jk})=\frac{1}{\sqrt{2\pi\sigma^2}}exp[-\frac{1}{2\sigma^2}(y_{jk}-\mu_{jk})^2]\]
\(H^{1}\)のもとでモデルを当てはめると:
\[l_{1}(\alpha_{1}, \alpha_{2}, \beta_{1}, \beta_{2})=\sum_{j=1}^{J}\sum_{k=1}^{K}[-\frac{1}{2}log(2\pi\sigma^2)-\frac{1}{2\sigma^2}(y_{jk}-\mu_{jk}^2)]=-\frac{1}{2}JKlog(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_{j=1}^{J}\sum_{k=1}^{K}(y_{jk}-\alpha_{j}-\beta_{j}x_{jk})^2\]
ここでJ=2, K=12であり,パラメータ\(\sigma^2\)は既知または局外パラメータ(nuisance parameter)として扱い,推定しない
そのうえで\(\alpha_{1}\),\(\alpha_{2}\),\(\beta_{1}\),\(\beta_{2}\)を推定する
最尤推定値の解は以下の連立方程式の解として与えられる: \[\frac{\delta l_{1}}{\delta\alpha_{j}}=\frac{1}{\sigma^2}\sum_{k}(y_{jk}-\alpha_{j}-\beta_{j}x_{jk})=0\] \[\frac{\delta l_{1}}{\delta\beta_{j}}=\frac{1}{\sigma^2}\sum_{k}x_{jk}(y_{jk}-\alpha_{j}-\beta_{j}x_{jk})=0\]
最尤推定に代替する方法として最小二乗法(MLS)を用いる
\(H^{1}\)のもとでのモデルにおいて: \[S_1=\sum_{j=1}^{J}\sum_{k=1}^{K}(y_{jk}-\mu_{jk})^2=\sum_{j=1}^{J}\sum_{k=1}^{K}(y_{jk}-\alpha_j-\beta_jx_{jk})^2\] これを最小にすることを考えるので, 最小二乗推定値は以下の式を解くことで求められる:
\[\frac{\delta S_{1}}{\delta\alpha_{j}}=-2\sum_{k=1}^{K}(y_{jk}-\alpha_{j}-\beta_{j}x_{jk})=0\]
\[\frac{\delta S_{1}}{\delta\beta_{j}}=-2\sum_{k=1}^{K}x_{jk}(y_{jk}-\alpha_{j}-\beta_{j}x_{jk})=0\]
上記の推定方程式を以下のように簡素化:
\[\sum_{k=1}^{K}y_{jk}-K\alpha_j-\beta_j\sum_{k=1}^{K}x_{jk}=0\] \[\sum_{k=1}^{K}x_{jk}y_{jk}-K\alpha_j\sum_{k=1}^{K}-\beta_j\sum_{k=1}^{K}x^2_{jk}=0\]
*j=1, 2
これらの式を正規方程式(normal equations)とよび, 解は
\[b_j=\frac{K\sum_k x_{jk}y_{jk}-(\sum_k x_{jk})(\sum_k y_{jk})}{K\sum_k x^2_{jk}-(\sum_k x_{jk})^2}\]
\[\alpha_j=\bar{y_{j}}-b_j\bar{x_{j}}\]
最小化すべき最小2乗基準は:
\[S_0=\sum_{j=1}^{J}\sum_{k=1}^{K}(y_{jk}-\alpha_j-\beta x_{jk})^2\]
上式より最小2乗推定値は以下の連立方程式を解くことで得られる:
\[\frac{\delta S_{0}}{\delta\alpha_{j}}=-2\sum_{k=1}^{K}(y_{jk}-\alpha_{j}-\beta x_{jk})=0\]
\[\frac{\delta S_{0}}{\delta\beta}=-2\sum_{j=1}^{J}\sum_{k=1}^{K}x_{jk}(y_{jk}-\alpha_{j}-\beta x_{jk})=0\]
ここでj=1,2 であり,解は以下のようになる
\[b=\frac{K\sum_j\sum_k x_{jk}y_{jk}-\sum_j(\sum_k x_{jk}\sum_k y_{jk})}{K\sum_j\sum_k x^2_{jk}-\sum_j(\sum_k x_{jk})^2}\]
\[a_j=\bar{y_{j}}-b\bar{x_{j}}\]
\(H_0を真とした場合のモデル\)
m0 <- lm(weight~duration+sex,data=d)
par(mfrow=c(2,2))
plot(m0)
\(H_1を真とした場合のモデル\)
m1 <-lm(weight~duration*sex,data=d)
par(mfrow=c(2,2))
plot(m1)
これらの図から,以下のことがいえる:
1. 標準化残差~当てはめ値, 標準化残差~説明変数(妊娠期間)のいずれの散布図にも系統的パターンは認められない
2. Normal Q-Q plotより, 標準化残差は近似的に正規分布に従っているといえる
3. 2つのモデルの間にはほとんど差がない
「2つのモデルの間にはほとんど差がない」ことを\(H_1\)に対する\(H_0\)の検定によって調べる
予想:\(H_1\)が正しければ,\(\hat{S}_0\)は\(\hat{S}_1\)に比べて明らかに大きくなるはず
\(\hat{S}_1\)は\(\hat{S}_0\)の相対的大きさを評価するために以下の2つの確率変数の標本分布を利用する必要あり
\[\hat{S}_1=\sum_{j=1}^{J}\sum_{k=1}^{K}(Y_{jk}-\alpha_j-b_jx_{jk})^2\]
\[\hat{S}_0=\sum_{j=1}^{J}\sum_{k=1}^{K}(Y_{jk}-\alpha_j-b x_{jk})^2\]
このとき,
\[\hat{S}_1=\sum_{j=1}^{J}\sum_{k=1}^{K}[Y_{jk}-(\alpha_j+b_jx_{jk})]^2-K\sum_{j=1}^{J}(\bar{Y_j}-\alpha_j-\beta_j\bar{x_j})^2-\sum_{j=1}^{J}(b_j-\beta_j)^2(\sum_{k=1}^{K}x^2_{jk}-K\bar{x^2_j})\]
となり, 確率変数\(Y_{jk}\), \(\bar{Y_j}\), \(b_j\)は次の分布に従う
\[Y_{jk}~N(\alpha_j+\beta_j x_{jk}, \sigma^2)\]
\[\bar{Y_j}~N(\alpha_j+\beta_j\bar{x_j}, \sigma^2/K)\]
\[b_j~N(\beta_j, \sigma^2/(\sum_{k=1}^{K}x^2_{jk}-K_bar{x^2_j}))\]
このとき, \(\hat{S_1}/\sigma^2\)は, 標準正規分布に従う確率変数の2乗和の差の形になっている
カイ2乗分布の性質から,
\[\hat{S_1}/\sigma^2~\chi^2(JK-2J)\]が導かれる
また\(H^0\)のもとでは,
\[\hat{S_0}/\sigma^2~\chi^2[JK-(J+1)]\]となる
ここで,統計量\(\hat{S}_0\)-\(\hat{S}_1\)は,\(H^0\)の下でのモデルに対する\(H^1\)のもとでのモデルの当てはまりの改良の大きさを表す
\(H_0\)が真のとき,
\[\frac{1}{\sigma^2}(\hat{S}_0-\hat{S}_1)~\chi^2(J-1)\]
\(H_1\)が真のとき,
\[\frac{1}{\sigma^2}(\hat{S}_0-\hat{S}_1)~非心カイ2乗分布\]
このとき\(\sigma^2\)は未知なので, 比をとることで\(\sigma^2\)を取り除く
\[F=\frac{(\hat{S_0}-\hat{S_1})/\sigma^2}{(J-1)}/\frac{\hat{S_1}/\sigma^2}{(JK-2J)}=\frac{(\hat{S_0}-\hat{S_1})/(J-1)}{\hat{S_1}/(JK-2J)}\]
\(H_0\)が真のとき, Fは中心F分布 \(F(J-1, JK-2J)\)に従う
\(H_1\)が真のとき, Fは非心F分布に従い, 計算されたF値は中心F分布から期待されるものより大きくなると予想される
anova(m0,m1,test="F")
## Analysis of Variance Table
##
## Model 1: weight ~ duration + sex
## Model 2: weight ~ duration * sex
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 658771
## 2 20 652425 1 6346.2 0.1945 0.6639
実際に計算してみると, F値は0.19となり, F(1, 20)分布において有意ではない(p=0.66)であるとわかる
すなわちこのデータは\(H_0\)を否定する根拠を与えない