【統計】時系列分析（ARIMAモデル）

統計 R、R Studio

本記事は、ARIMAモデル等の特徴について。
【目次】

計算式等
ノート
参考

計算式等

＜前提＞

■ 弱定常性（Weakly stationary）
　時系列の期待値、分散および自己共分散が常に一定である状態（時点 t および k に依存しない状態）。単に定常性と言ったらこちらを指す。

$\quad E[y_{t}]=\mu_{t}=\mu_{t-k}=\mu$
$\quad Var[y_{t}]=\sigma^{2} _{t}=\sigma^{2} _{t-k}=\sigma^{2}$
$\quad Cov \left(y_{t},\ y_{t-s} \right) = Cov \left(y_{t-k},\ y_{t-s-k} \right)$

■ 強定常性（Strongly stationary）
　任意の時点の集合 {t1, t2, ..., tm} と k に対する時系列の同時密度関数 f が一定の状態（同時分布が同一となる場合）。

$\quad f \left(y_{t_{1}} \ _{,}\ y_{t_{2}} \ _{,}\ ...\ _{,}\ y_{t_{m}} \right) = f \left( y_{t_{1-k}}\ _{,}\ y_{t_{2-k}}\ _{,}\ ...\ _{,}\ y_{t_{m-k}} \right)$

■ ホワイトノイズ（White noise、白色雑音）
　時系列の誤差項。正規分布であり、他の時点と自己相関がない。

$\quad E[\varepsilon _{t}]=0, \quad Var[\varepsilon _{t}]=\sigma^{2}, \quad Cov[\varepsilon _{t},\ \varepsilon _{t-s}]=0$

■ ランダムウォーク（Random walk）
　ホワイトノイズの累積和により得られる時系列。

$\quad y_{t} = \sum \varepsilon _{t}$

＜各モデル＞

■ AR(p)モデル（AutoRegressive model、自己回帰モデル）
　ラグpまでの過去時点（AR(p) 成分）を説明変数に投入して現時点を表すモデル。pは自己回帰の次数、Φは自己回帰係数。

$\quad y_{t} = c + \sum ^{p}_{s=1} \phi _{s} y_{t-s} + \varepsilon _{t}, \quad \varepsilon _{t}\sim W.N.\left( \sigma ^{2}\right)\quad\verb|※W.N.=ホワイトノイズ|$

　※AR(1) の時：
$\qquad y_{t} = c + \phi _{1} y_{t-1} + \varepsilon _{t}$

■ MA(q)モデル（Moving Average model、移動平均モデル）
　ラグqまでの過去のホワイトノイズ（MA(q) 成分）を説明変数に投入して現時点を表すモデル。qは移動平均の次数、θは移動平均係数。

$\quad y_{t} = c + \varepsilon _{t} + \sum ^{q}_{s=1}\theta _{s}\varepsilon_{t-s},\quad \varepsilon _{t}\sim W.N.\left( \sigma ^{2}\right)$

　※MA(1) の時：
$\qquad y_{t} = c + \varepsilon _{t} + \theta _{1}\varepsilon_{t-1}$

■ ARMA(p, q)モデル（AutoRegressive Moving Average model、自己回帰移動平均モデル）
　AR(p) 成分とMA(q) 成分の両方を説明変数に投入して現時点を表すモデル。

$\quad y_{t} = c + \sum ^{p}_{s=1}\phi _{s}y_{t-s} + \varepsilon _{t} + \sum ^{q}_{s=1}\theta _{s}\varepsilon_{t-s},\quad \varepsilon _{t}\sim W.N.\left( \sigma ^{2}\right)$

　※ARMA(1,1) の時：
$\qquad y_{t} = c + \phi _{1}y_{t-1} + \varepsilon _{t} + \theta _{1}\varepsilon_{t-1}$

■ ARIMA(p, d, q)モデル（AutoRegressive Integrated Moving Average model、自己回帰和分移動平均モデル）
　観測値の時系列に対して、近似的に定常な時系列が得られるまで d階の階差（和分プロセス）をとった後に、階差の時系列に対して、ARMA(p,q) モデルを当てはめるモデル。

和分プロセス：
　[ ラグ i の１階差 ]
$\quad \Delta _{i} y_{t} = y_{t} - y_{t-i}$

　[ ラグ i, j の２階差 ]
$\quad \Delta _{i} \Delta _{j} y_{t} = \Delta _{j} y_{t} - \Delta _{j} y_{t-i} = \left(y_{t} - y_{t-j}\right) - \left(y_{t-i} - y_{t-i-j}\right)$

　※ラグ１の１階差の場合：
$\qquad \Delta_{1} y_{t} = y_{t} - y_{t-1}$

ノート

上に示した過程はARIMAモデルで表すことができる。

・ARIMA(0,0,0)：ホワイトノイズ
・ARIMA(0,1,0)：ランダムウォーク
・ARIMA(p,0,0)：AR(p)モデル
・ARIMA(0,0,q)：MA(q)モデル
・ARIMA(p,0,q)：ARMA(p,q)モデル
・ARIMA(p,d,q)：ARIMA(p,d,q)モデル

また、MA(q) モデルは時点 t に依存しないため常に定常だが、
AR(p) モデルは、次に示す特性方程式の解（複素解含む）の絶対値がすべて１未満であることが、定常性の必要十分条件と言われる。

$\quad 1 = \sum \phi_{s} z^{s} = \phi_{1} z^{1} + \phi_{2} z^{2} + ... + \phi_{s} z^{s}$

AR(1) モデルなら 1 ＜ |Φ1| の場合。つまり、Φ1 が１を超えると非定常になってしまう。

ここからは、自己回帰係数Φと平均移動係数θの違いによる振る舞いを見てみる。

次のように適当な定数項 c とホワイトノイズ（ =NORM.INV(RAND(),0,0.5) ）を準備し、n=150 の観測値の時系列を作成。
※Excelでの例示は雑なのでイメージ程度にとどめてほしい。また、ARモデルの初期値は０にしている。
f:id:cochineal19:20210307162430p:plain

この観測値の時系列に対して、AR(1) とMA(1) モデルを当てはめて、自己回帰係数Φと平均移動係数θによる変化を見てみる。

f:id:cochineal19:20210307171843p:plain
f:id:cochineal19:20210307171921p:plain
f:id:cochineal19:20210307171951p:plain
f:id:cochineal19:20210307172016p:plain
f:id:cochineal19:20210307172120p:plain
f:id:cochineal19:20210307172145p:plain
f:id:cochineal19:20210307172215p:plain
f:id:cochineal19:20210307172238p:plain

先ほどの AR(1) モデルの定常性条件（1 ＜ |Φ1| ）を頭に入れて見てみると、自己回帰係数Φの絶対値が１以上の Φ=1.1 と Φ=－1.1 で非定常になっていることが良く分かる（予測モデルとして使えない）。
また、MA(1) モデルは定数項 c 付近でバラつく定常な時系列であることが分かる。

次に、先ほどの観測値の時系列にトレンド（時点毎に 0.5 増える）を加えてみる。このトレンド付きの時系列について、ラグ１の１階差（d=1）をとったときの振る舞いを見てみる。
※こちらも例示は雑なのでイメージ程度にとどめてほしい。また、初期値は０にしている。
f:id:cochineal19:20210307165011p:plain

１階差をとった後、階差の時系列が定常になったことが分かる。

今回、AR(1) とMA(1) モデルでいくつかと１階差を１つ例示したが、ARIMA(p,d,q) でどのようなパラメタ設定が良いか（ベストモデルは何か）は、AICやRMSE等の評価指標を用いて判断する。
r では forecastパッケージでここらができるので次回記事にしたい。
SAS にも proc ARIMA があるので余裕があれば。

なお、ARIMA(p,d,q) に従う観測値の時系列は、r の stats::arima.sim関数で比較的簡単に作れるので備忘録がてらメモする。

stats::arima.sim(model=list(order=c(p,d,q), ar=c(x,...), ma=c(x,...), n=x))

AR(p)、MA(q) に従う観測値の時系列：

library(tidyverse)
set.seed(1)
ads1 <- data.frame()
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(1,0,0), ar=c(0.3)),         n=200), TYPE="01. AR(1)  Φ1 = 0.3"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(2,0,0), ar=c(0.3,0.5)),     n=200), TYPE="02. AR(2)  Φ1 = 0.3, Φ2 = 0.5"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(3,0,0), ar=c(0.3,0.5,0.1)), n=200), TYPE="03. AR(3)  Φ1 = 0.3, Φ2 = 0.5, Φ3 = 0.1"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(0,0,1), ma=c(0.3)),         n=200), TYPE="04. MA(1)  θ1 = 0.3"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(0,0,2), ma=c(0.3,0.5)),     n=200), TYPE="05. MA(2)  θ1 = 0.3, θ2 = 0.5"))
ads1 <- union_all(ads1, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(0,0,3), ma=c(0.3,0.5,0.1)), n=200), TYPE="06. MA(3)  θ1 = 0.3, θ2 = 0.5, θ3 = 0.1"))

ggplot(ads1, aes(x=x, y=AVAL)) +
  geom_line() + facet_wrap(. ~ TYPE, ncol=3)

f:id:cochineal19:20210307180214p:plain

ARMA(p,q) に従う観測値の時系列：

set.seed(1)
ads2 <- data.frame()
ads2 <- union_all(ads2, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(1,0,1), ar=c(0.3),         ma=c(0.3)),         n=200), TYPE="01. ARIMA(1,0,1) \n Φ1 = 0.3 \n θ1 = 0.3"))
ads2 <- union_all(ads2, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(2,0,2), ar=c(0.3,0.5),     ma=c(0.3,0.5)),     n=200), TYPE="02. ARIMA(2,0,2) \n Φ1 = 0.3, Φ2 = 0.5 \n θ1 = 0.3, θ2 = 0.5"))
ads2 <- union_all(ads2, data.frame(x=1:200, AVAL=arima.sim(model=list(order=c(3,0,3), ar=c(0.3,0.5,0.1), ma=c(0.3,0.5,0.1)), n=200), TYPE="03. ARIMA(3,0,3) \n Φ1 = 0.3, Φ2 = 0.5, Φ3 = 0.1 \n θ1 = 0.3, θ2 = 0.5, θ3 = 0.1"))
ggplot(ads2, aes(x=x, y=AVAL)) +
  geom_line() + facet_wrap(. ~ TYPE, ncol=3)

f:id:cochineal19:20210307181223p:plain
ARIMA(p,d,q) に従う観測値の時系列：

set.seed(1)
ARIMA212 <- data.frame(y=arima.sim(model=list(order=c(2,1,2), ar=c(0.3,0.5), ma=c(0.3,0.5)), n=200), id=1:201)
ggplot(data=ARIMA212, mapping=aes(x=id, y=y)) +
  geom_line() +
  labs(title="ARIMA(2,1,2)", subtitle="（AR: Φ1 = 0.3, Φ2 = 0.5,  MA: θ1 = 0.3, θ2 = 0.5）"
       ,x="時点", y="", caption="🄫Cochineal19") +
  theme(plot.title=element_text(size=30)
        ,plot.subtitle=element_text(size=25)
        ,axis.title=element_text(size=25)
        ,axis.text=element_text(size=20))
 
set.seed(1)
ARIMA222 <- data.frame(y=arima.sim(model=list(order=c(2,2,2), ar=c(0.3,0.5), ma=c(0.3,0.5)), n=200), id=1:202)
ggplot(data=ARIMA222, mapping=aes(x=id, y=y)) +
  geom_line() +
  labs(title="ARIMA(2,2,2)", subtitle="（AR: Φ1 = 0.3, Φ2 = 0.5,  MA: θ1 = 0.3, θ2 = 0.5）"
       ,x="時点", y="", caption="🄫Cochineal19") +
  theme(plot.title=element_text(size=30)
        ,plot.subtitle=element_text(size=25)
        ,axis.title=element_text(size=25)
        ,axis.text=element_text(size=20))

f:id:cochineal19:20210307182908p:plain:w400
f:id:cochineal19:20210307182949p:plain:w400

本日は時系列の特徴について。ここまで。

参考

統計学（出版：東京図書）, 日本統計学会編
カルマンフィルタ Rを使った時系列予測と状態空間モデル（出版：共立出版）, 野村俊一
http://elsur.jpn.org/202004MMM/chap5.pdf
http://yuchimatsuoka.github.io/seminar/timeseries1.pdf

★サイトマップ

cochineal19.hatenablog.com

2021-02-28

【統計】時系列分析（自己相関、自己共分散）

統計 R、R Studio

自己相関、自己共分散について、視覚的に分かるようにメモ。

【目次】

計算式
ノート
参考

自己相関、「自己」とつくので、ちょっと分かりづらい。
自分自身に対する相関なのだなと何となく分かるが、それで？となる。

計算式

■ 自己共分散
$\quad \gamma\left( s \right)=Cov \left(x_{t},\ x_{t-s} \right) =\dfrac{1}{N} \sum \{ \left(x_{t} - \overline{x} \right) \left(x_{t-s} - \overline{x} \right) \}$

■ 自己相関
$\quad \rho\left( s \right) = \dfrac{\gamma\left( s \right)}{\gamma\left( 0 \right)} = \dfrac{Cov \left(x_{t},\ x_{t-s} \right)}{Cov \left(x_{t},\ x_{t} \right)} = \dfrac{Cov \left(x_{t},\ x_{t-s} \right)}{Var \left(x_{t} \right) }$
$\qquad\quad = \dfrac{\dfrac{1}{N}\sum \{ \left(x_{t} - \overline{x} \right) \left(x_{t-s} - \overline{x} \right) \}}{\dfrac{1}{N}\sum \left(x_{t} - \overline{x} \right)^{2}} = \dfrac{\sum \{ \left(x_{t} - \overline{x} \right) \left(x_{t-s} - \overline{x} \right) \}}{\sum \left(x_{t} - \overline{x} \right)^{2}}$

　※t=時期、s=ラグ

ノート

サンプルデータで動きを見てみる。
Rで準備したデータ。３期で１サイクルの何らかのデータ。特に値に意味はない。

ads <- data.frame(
  AVAL=c(10,8,6,10,7,5,9,8,5,10,8,5)
)

f:id:cochineal19:20210227222714p:plain:w300

自己相関と自己共分散は acf 関数 で計算できる。
plot=FALSEにするとコンソールに自己相関などの値を出してくれ、
plot=TRUEにするとコレログラムを描画してくれる。

acf(ads, plot=F, type="correlation")  #-- 自己相関
acf(ads, plot=F, type="covariance")  #-- 自己共分散
#-- 自己偏相関もある（type="partial"）

plot=FASEで出した結果。

> acf(ads, plot=F, type="correlation")
Autocorrelations of series ‘ads’, by lag
     0      1      2      3      4      5      6      7      8      9     10 
 1.000 -0.322 -0.466  0.695 -0.208 -0.436  0.465 -0.032 -0.283  0.235 -0.002 

> acf(ads, plot=F, type="covariance")
Autocovariances of series ‘ads’, by lag
       0        1        2        3        4        5        6        7        8        9       10 
 3.57639 -1.15336 -1.66782  2.48438 -0.74537 -1.55845  1.66319 -0.11516 -1.01157  0.84201 -0.00579

plot=TRUEで出したコレログラム（自己相関のみ掲載）。
f:id:cochineal19:20210227185549p:plain:w300

これらは時期 t とそこから s 期分ずらした時期 t-s（LAG s）の AVAL について相関や共分散を出している。
つまり、過去のどの時点と関係性があるかを示すことができる。
今回の架空データは３期で１サイクルのため、特に LAG 3 に中程度以上の相関がみられる。
LAG 0 は現時点と現時点を見比べるため、当然ながら自己相関＝１になる。
これを見ることで周期性を探れる。

Excelで計算方法を見てみる。
f:id:cochineal19:20210228004053p:plain 緑列は基本情報。
timeは時間単位、AVALは本記事の最初に示したもの。
偏差と偏差平方は、例えば、time=1 では AVAL=10、AVALの平均≒7.58 のため、
偏差=(10 - 7.58)=2.42、偏差平方=2.42²≒5.84 になる。

黄色列は偏差積。
例えば、LAG 1 の time=2 は１期前の time=1 との偏差積をとり、0.42×2.42≒1.007 になる。
LAG 2 の time=3 は、２時期ずらして time=1 との偏差積をとり、(-1.6)×2.42≒-3.826 となる。

最後に青列は、自己共分散と自己相関。
黄色列の LAG s をそれぞれ総和（偏差積和）して、N=12 で割ったものが自己共分散。
LAG0 は分散 Var[AVAL]≒3.5764 になる。
そして、LAG s の偏差積和 ÷ AVALの偏差平方和、もしくは LAG s の共分散 ÷ 分散 Var[AVAL] とすれば、LAG s の自己相関になる。

今回はここまで。

参考

統計学（出版：東京図書）, 日本統計学会編

2021-02-27

【R、機械学習メモ】Rでのデータ分割

機械学習 R、R Studio

Rでのデータ分割の方法のメモ。

データ準備
rsample で分割
dplyr で分割
参考

データ準備

1500行の何の意味もないデータ。

ads <- data.frame(id=seq(1:1500))

rsample で分割

library(rsample)
ads.sp1  <- rsample::initial_split(ads, prop=3/4) 
adtrain1 <- rsample::training(ads.sp1)  #-- 訓練データ
adtest1  <- rsample::testing(ads.sp1)  #-- テストデータ

initial_split関数で分割する割合を決める。
training関数とtesting関数で訓練データとテストデータを作る。

dplyr で分割

library(tidyverse)
adtrain2 <- dplyr::sample_frac(ads, size=3/4, replace=FALSE)  #-- 訓練データ
adtest2  <- dplyr::anti_join(ads, adtrain2, by=c("id"))  #-- テストデータ。byはユニークキー

sample_frac関数で訓練データ分を作る。
anti_join関数でテストデータ（訓練データ以外のデータセット）を作る。

参考

Simple Training/Test Set Splitting — initial_split • rsample

2021-02-23

【統計】共分散分析（ANCOVA）

R、R Studio SAS Python 統計

本記事では、共分散分析について記載します。

【目次】

帰無仮説、対立仮設
計算式等
プログラムコード
雑談
参考

共分散分析は、分散分析と重回帰分析を組み合わせた手法と言われます。
分散分析では説明変数に集団などカテゴリー変数（要因）のみを投入しますが、共分散分析ではカテゴリー変数（要因）と連続変数（共変量）の両方を投入するのが特徴です。

分散分析を ANOVA (ANalysis Of VAriance) と言うのに対し、共分散分析は ANCOVA (ANalysis of COVAriance) と呼称されます。

帰無仮説、対立仮設

帰無仮説 $H_{0}:$ 要因の水準間の主効果が同じ（平均に差がない）
対立仮設 $H_{1}:$ 要因の水準間の主効果に差がある（平均に差がある）

計算式等

ANOVAでは、水準間の変動が残差の変動（ランダムな変動）より有意に大きいかどうかで、水準間の差を評価します。
ANCOVAでは、目的変数と相関（線形関係）を持つ共変量を投入し、その影響を調整（共通の傾きを持つ回帰直線を作成）した上で、水準間の変動が残差の変動より有意に大きいかを評価します。

以下、回帰式を記載しますが、重回帰分析の記事と同じです。
※カテゴリー変数１つ（２水準）、連続変数（共変量）１つの例としています。

■ 水準毎の回帰直線（水準毎の条件付き期待値）
$\quad \widehat{ys}_{ki} = \beta_{0_{k}} + \beta_{1_{k}}c_{ki}$
　※ k：要因の水準、c：共変量
　※ 例えば、要因の水準がA薬・B薬の２つなら薬剤毎に計２つの回帰直線を作成。

■ 共通の回帰直線（共通の条件付き期待値）
$\quad \widehat{yc}_{i} = \gamma_{0} + \gamma_{1}c_{i} + \gamma_{2} x_{i}$
　※ c：共変量、x：要因
　※ 共変量 c の傾きを一定として、要因 x による違いを表現した直線。
　※ 例えば、要因の水準がA薬・B薬の２つなら、x=0、x=1の2値で $\gamma_{2}$ 分が期待値の差になる。

■ 共通の回帰直線の重心（総平均）
$\quad \delta_{0} = \sum \{ \left( \gamma_{0} + \gamma_{2} x_{k} \right) N_{k} \} / N$
$\quad G = \delta_{0} + \gamma_{1} \overline{c}$
　※ $\gamma_{2} x_{k}$ ：水準kの回帰係数（xは０or１のダミー変数）、 $N_{k}$ =水準kのN数

■ 分散分析表

変動要因		平方和	自由度	平均平方	F値
B ：共通傾きでの回帰直線(水準毎) vs. 共通の回帰直線の重心		$S_{B}=\sum\left( \widehat{yc}_{i} - G \right)^{2}$	$df_{B}=a-1$	$V_{B}=S_{B}/df_{B}$	$F_{B}=V_{B}/V_{W}$
W：観測値 vs. 共通傾きでの回帰直線(水準毎)		$S_{W}=\sum\left( y_{i} - \widehat{yc}_{i} \right)^{2}$	$df_{W}=df_{T}-df_{B}$	$V_{W}=S_{W}/df_{W}$	$F_{B}=V_{B}/V_{W}$
	R：回帰直線(水準毎) vs. 共通傾きでの回帰直線(水準毎)	$S_{R}=\sum\sum\left( \widehat{ys}_{ik} - \widehat{yc}_{i} \right)^{2}$	$df_{R}=a-1$	$V_{R}=S_{R}/df_{R}$	$F_{R}=V_{R}/V_{e}$
	e：観測値 vs. 回帰直線(水準毎)	$S_{e}=\sum\sum\left( y_{i} - \widehat{ys}_{ik} \right)^{2}$	$df_{e}=df_{T}-df_{R}$	$V_{e}=S_{e}/df_{e}$	$F_{R}=V_{R}/V_{e}$
T ：観測値 vs. 共通の回帰直線の重心		$S_{T}=\sum\left( y_{i} - G \right)^{2}$	$df_{T}=n-1$	$V_{T}=S_{T}/df_{T}$	---

　※k=水準k、a=要因の水準数＋共変量の数。

（各平方和の図解）

f:id:cochineal19:20210222182649p:plain

■ 調整平均（LS mean：Least Square mean）

水準毎の調整平均　※要因 x を水準AとBで考えた場合。
$\quad \overline{yc}_{A} = \gamma_{0} + \gamma_{1}\overline{c} + \gamma_{2} x_{A}$
$\quad \overline{yc}_{B} = \gamma_{0} + \gamma_{1}\overline{c} + \gamma_{2} x_{B}$
調整平均の差とその信頼区間
$\quad d_{y_{\verb|adj|}} = \gamma_{2}$
$\quad CL=d_{y_{\verb|adj|}} \pm \verb|T.INV|\left(α,\ n-3\right) \cdot se[\gamma_{2}]$
　※ $se[\gamma_{2}]$ ：共通の回帰直線における、共変量の回帰係数 $\gamma_{2}$ の標準誤差（Standard Error）

（補足）通常の平均の差とその信頼区間
$\quad d_{y}=\overline{y_{B}} - \overline{y_{A}}$
$\quad CL=d_{y} \pm \verb|T.INV|\left(α,\ n-2\right) \cdot \sqrt{\dfrac{\sigma ^{2}}{n_{A}}+\dfrac{\sigma ^{2}}{n_{B}}}$

ノート

今回も架空データで計算手順を確認します。
　BASE = {21, 15, 18, 16, 26, 25, 22, 21, 16, 17, 18}
　TRT01AN = {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1}
　CHG = {-7, -2, -5, -4, -12, -15, -12, -12, -4, -7, -7}
　AVAL = {14, 13, 13, 12, 14, 10, 10, 9, 10, 10, 11}

特定時点における検査Xのベースラインからの変化量 (CHG) を目的変数として、A薬（TRT01AN=0）とB薬（TRT01AN=1）による効果の違いを評価するシナリオとします。この時、共変量として検査Xのベースライン値 (BASE) を投入します。臨床試験で見かけるもの。

■ 分析モデル：
　ベースラインからの変化量＝投与群（要因）＋ベースライン値（共変量）

回帰の同質性の検定

先ず、水準毎の回帰直線（CHG=BASE）を描画してみます。
図のとおり、若干傾きが異なり、A薬とB薬の差をどこで見るかで結果が変わってしまうことが分かります。
$\quad \hat{ys_{Ai}}=10.325 - 0.850 c_{Ai}$
$\quad \hat{ys_{Bi}}=11.011 - 1.051 c_{Bi}$

f:id:cochineal19:20210222231956p:plain

ここに、共通の傾きの回帰直線（CHG=BASE+TRT01AN）を重ねます。
この赤・青の実線であれば、共変量の値によらずA薬とB薬の差を一定に保つことができます。
$\quad \hat{yc_{i}} = 11.972 - 0.936 c_{i} - 3.240 x_{i}$

f:id:cochineal19:20210223115640p:plain

共分散分析では、この共通の傾きの回帰直線（赤・青の実線）を用いますが、そのためには、水準毎の回帰直線（赤・青の点線）が平行であるという前提が必要です。
この前提の検証として、下図のとおり、 $S_{W}$ （共通の傾きの回帰直線からの残差平方和）を $S_{R}$ （水準毎の回帰直線からの残差平方和）と $S_{e}$ （２本の回帰直線の偏差平方和）に分解して F検定を行うことで、水準毎の回帰直線が平行と解釈できるかを評価します（回帰の同質性の検定、帰無仮説：水準毎の回帰直線が平行）。

f:id:cochineal19:20210222233701p:plain

それでは計算に入ります。先ほどの回帰式を個々の値に当てはめて、期待値と偏差平方を出したものが下表です。
※Y共通傾き列：共通の傾きを用いた各水準の回帰直線での期待値、Y単体列：各水準単体の回帰直線での期待値、R列：Y共通傾きとY単体の偏差平方、e列：YCHGとY単体の偏差平方
＜再掲＞
$\quad \hat{ys_{Ai}}=10.325 - 0.850 c_{Ai}$
$\quad \hat{ys_{Bi}}=11.011 - 1.051 c_{Bi}$
$\quad \hat{yc_{i}} = 11.972 - 0.936 c_{i} - 3.240 x_{i}$

C_BASE	X_TRT01AN	Y_CHG	Y_共通傾き	Y_単体	R	e
21	0	-7	-7.68	-7.53	0.02	0.28
15	0	-2	-2.07	-2.43	0.13	0.18
18	0	-5	-4.88	-4.98	0.01	0.00
16	0	-4	-3.00	-3.28	0.08	0.52
26	0	-12	-12.37	-11.78	0.34	0.05
25	1	-15	-14.67	-15.26	0.35	0.07
22	1	-12	-11.86	-12.11	0.06	0.01
21	1	-12	-10.93	-11.06	0.02	0.88
16	1	-6	-6.25	-5.80	0.19	0.04
17	1	-7	-7.18	-6.86	0.11	0.02
18	1	-7	-8.12	-7.91	0.04	0.82

R列、e列をそれぞれ足し合わせ平方和を算出し、F値、p値を求めます。

変動要因	平方和	自由度	平均平方	F値	p値
R：回帰直線(水準毎) vs. 共通傾きでの回帰直線(水準毎)	1.357	2	0.679	1.4139	0.3140
e：観測値 vs. 回帰直線(水準毎)	2.880	6	0.480	1.4139	0.3140

p > 0.05 で非有意であれば、水準毎の回帰直線は平行であると解釈して、以降、共通の傾きでの回帰直線を用いて共分散分析を行います。
今回の架空データでは p=0.3140で非有意のため、A薬・B薬の回帰直線は平行と解釈し、共分散分析に進みます。

（※ 水準毎の回帰直線が平行であることの評価方法として、交互作用項を含めたモデルを作り、交互作用項が非有意なら平行と解釈する方法もあります。雑談に回します）

共分散分析

先ず、共通の回帰直線における重心（総平均）を考えます。
※今回、A薬はN=5, B薬はN=6の全体N=11。A薬を x=0、B薬を x=1 としています。
$\quad \delta_{0} = \sum \{ \left( \gamma_{0} + \gamma_{2} x_{k} \right) N_{k} \} / N=\{ 11.972 \times 5 + \left(11.972 - 3.240 \right) \times 6 \} / 11 ≒ 10.205$
$\quad G = \delta_{0} + \gamma_{1} \overline{c} = 10.205 - 0.936 \times 19.546 ≒ -8.091$

重心が算出できたら同質性の検定時と同じ要領で偏差平方を求めます。
※T列：YCHGと重心との偏差平方、B列：Y単体と重心との偏差平方、W列：YCHGとY共通傾きの偏差平方

f:id:cochineal19:20210223002545p:plain

C_BASE	X_TRT	Y_CHG	AVAL	Y_共通傾き	Y_単体	T	B	W
21	0	-7	14	-7.68	-7.53	1.19	0.16	0.47
15	0	-2	13	-2.07	-2.43	37.10	36.27	0.00
18	0	-5	13	-4.88	-4.98	9.55	10.33	0.02
16	0	-4	12	-3.00	-3.28	16.74	25.87	0.99
26	0	-12	14	-12.37	-11.78	15.28	18.27	0.13
25	1	-15	10	-14.67	-15.26	47.74	43.28	0.11
22	1	-12	10	-11.86	-12.11	15.28	14.22	0.02
21	1	-12	9	-10.93	-11.06	15.28	8.03	1.15
16	1	-6	10	-6.25	-5.80	4.37	3.41	0.06
17	1	-7	10	-7.18	-6.86	1.19	0.83	0.03
18	1	-7	11	-8.12	-7.91	1.19	0.00	1.25

T列、B列、W列をそれぞれ足し合わせ平方和を算出し、F値、p値を求めます。

変動要因		平方和	自由度	平均平方	F値	p値
B ：共通傾きでの回帰直線(水準毎) vs. 共通の回帰直線の重心		160.672	2	80.336	151.6721	0.0000
W：観測値 vs. 共通傾きでの回帰直線(水準毎)		4.237	8	0.530	151.6721	0.0000
	R：回帰直線(水準毎) vs. 共通傾きでの回帰直線(水準毎)	1.357	2	0.679	1.4139	0.3140
	e：観測値 vs. 回帰直線(水準毎)	2.880	6	0.480	1.4139	0.3140
T ：観測値 vs. 共通の回帰直線の重心		164.909	10	16.491	---	---

　※薄黄色は先ほどの同質性の検定の部分です。

この表の $S_{B}$ （水準間の平方和）と $S_{W}$ （共通の傾きの回帰直線からの残差平方和）の平均平方を比較することで、水準間の変動がランダムな変動より有意に大きいかを評価します。今回の架空データでは p < 0.001 で水準間に有意な変動があるようでした。

（追記）
SASの Output の Type II または III を見ると F (1,1)=53.64, p<0.0001 で薬剤（TRT01AN）の主効果が有意だったことが分かります。Type X 平方和は、共分散分析モデルの要因・共変量（TRT01AN、BASE）を分解して、要因別の主効果の有無を評価したもの。

f:id:cochineal19:20210223215534p:plain

※ Type II, III 平方和の計算は省略します。平方和の違いはいつかまとめたい。
※ Type I 平方和のTRT01ANは次のとおり。要否別で備忘録として。
$\quad \left( \overline{y_{\verb|A薬|}} - G \right)^{2} \cdot N_{\verb|A薬|} + \left( \overline{y_{\verb|B薬|}} - G \right)^{2} \cdot N_{\verb|B薬|} = \left( -6.000 - \left(-8.091 \right) \right)^{2} \cdot 5 + \left( -9.833 - \left(-8.091 \right) \right)^{2} \cdot 6 ≒ 40.0758$

調整平均（LS mean：Least Square mean）

共分散分析と一緒に調整平均の差とその信頼区間を示すこともありますので、備忘録がてらメモします。
今回の架空データを ExcelのLINEST関数で実行した結果がこちらです：

f:id:cochineal19:20210223104601p:plain
また、共変量（BASE）の平均は19.545だったため、調整平均は以下となります。
水準毎の調整平均
$\quad \overline{yc} = \gamma_{0} + \gamma_{1}\overline{c} + \gamma_{2} x=11.972-0.936\times19.545 + \begin{pmatrix} -3.240\times 0 \\ -3.240\times 1 \end{pmatrix}=\begin{pmatrix} -6.323 \\ -9.564\end{pmatrix}$
調整平均の差とその信頼区間
$\quad d_{y_{\verb|adj|}} = \gamma_{2}=-3.240$
$\quad CL=d_{y_{\verb|adj|}} \pm \verb|T.INV|\left(0.05, 8\right) \cdot se[\gamma_{2}]=-3.240 \pm 2.3060 \times 0.4424=[ -4.2608,\ -2.2202]$

これを通常の平均と比べると下表のとおりです。

評価項目	A薬	B薬	差 (B-A)	95%信頼区間
Y_CHGの平均	-6.000	-9.833	-3.833	-8.9349	1.2682
Y_CHGの調整平均（LS mean）	-6.323	-9.564	-3.240	-4.2608	-2.2202

今回の架空データでは、通常の平均の差の信頼区間は０を挟むのに対し、調整平均では信頼区間の幅が狭まり、０を挟まなくなったことが分かります（信頼区間下限でもB薬の方が効果を示している）。

Rでの実行：

library(tidyverse)
library(car)

#-- サンプルデータ
ADS <- data.frame(
  TRT01AN=c(0,0,0,0,0,1,1,1,1,1,1)
  ,BASE=c(21,15,18,16,26,25,22,21,16,17,18)
  ,AVAL=c(14,13,13,12,14,10,10,9,10,10,11))
ADS$CHG <- ADS$AVAL - ADS$BASE
ADS$TRT01AF <- relevel(factor(ifelse(ADS$TRT01AN==0,"A薬","B薬")), ref="A薬")

#-- 水準毎の回帰分析
ADS.A <- ADS[ADS$TRT01AF=="A薬",]; ANOVA1.A <- lm(CHG ~ BASE, data=ADS.A); summary(ANOVA1.A);
ADS.B <- ADS[ADS$TRT01AF=="B薬",]; ANOVA1.B <- lm(CHG ~ BASE, data=ADS.B); summary(ANOVA1.B)
g1 <- ggplot() +
  geom_point(data=ADS, mapping=aes(x=BASE, y=CHG, color=TRT01AF), size=5) +
  geom_smooth(data=ADS.A, mapping=aes(x=BASE, y=CHG), method="lm", se=F, color="red", lty=2) +
  geom_smooth(data=ADS.B, mapping=aes(x=BASE, y=CHG), method="lm", se=F, color="blue", lty=2) +
  labs(x="ベースラインの検査値（共変量）", y="ベースラインからの変化量"
       ,title="検査Xのベースラインからの変化量", color="投与群",caption="🄫Cochineal19") +
  theme(plot.title=element_text(size=30)
        ,axis.title=element_text(size=25)
        ,axis.text=element_text(size=20)
        ,legend.title=element_text(size=20)
        ,legend.text=element_text(size=20))
plot(g1)

#-- 交互作用項が非有意かどうかで平行性を評価。
ANCOVA0 <- lm(CHG ~ BASE + TRT01AF + BASE*TRT01AF, data=ADS);
summary(ANCOVA0)
car::Anova(ANCOVA0)

#-- ANCOVA
ANCOVA1 <- lm(CHG ~ BASE + TRT01AF, data=ADS)
(res <- summary(ANCOVA1))
car::Anova(ANCOVA1) #-- Type 2 平方和

#-- 共通の回帰直線の重心（総平均）
(b0.adj <- (res$coefficients[1] * nrow(ADS.A) + (res$coefficients[1]+res$coefficients[3]) * nrow(ADS.B)) / nrow(ADS))
(yint <- b0.adj + mean(ADS$BASE) * res$coefficients[2])

#-- 共通傾きでの回帰直線（水準毎）
ANCOVA1.A <- predict(ANCOVA1, newdata=ADS.A, se=F)
ANCOVA1.B <- predict(ANCOVA1, newdata=ADS.B, se=F)

#-- 共通の回帰直線と共通傾きの回帰直線（水準毎）をオーバーラップ
g2 <- g1 +
  geom_smooth(data=ADS, mapping=aes(x=BASE, y=ANCOVA1$fitted.values), stat="smooth", method="lm", se=F, color="darkgreen") +
  geom_smooth(data=ADS.A, mapping=aes(x=BASE, y=ANCOVA1.A), stat="smooth", method="lm", se=F, color="red") +
  geom_smooth(data=ADS.B, mapping=aes(x=BASE, y=ANCOVA1.B), stat="smooth", method="lm", se=F, color="blue") 
  geom_hline(yintercept=yint, lty=2, color="gray60") +
  geom_vline(xintercept=mean(ADS$BASE), lty=2, color="gray60")
plot(g2)

Rの ANCOVAテーブル

> #-- 交互作用項が非有意かどうかで平行性を評価。
> car::Anova(lm(CHG ~ BASE + TRT01AF + BASE*TRT01AF, data=ADS))
Anova Table (Type II tests)
Response: CHG
              Sum Sq Df F value    Pr(>F)    
BASE         120.596  1 293.113 5.692e-07 ***
TRT01AF       28.413  1  69.058 7.144e-05 ***
BASE:TRT01AF   1.357  1   3.299    0.1122    
Residuals      2.880  7                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> #-- ANCOVA
> car::Anova(ANCOVA1) #-- Type 2 平方和
Anova Table (Type II tests)
Response: CHG
           Sum Sq Df F value    Pr(>F)    
BASE      120.596  1 227.682 3.680e-07 ***
TRT01AF    28.413  1  53.642 8.196e-05 ***
Residuals   4.237  8                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SASでの実行：

data ADS;
input BASE TRT01AN CHG AVAL 8. @@;
cards;
21 0 -7 14
15 0 -2 13 
18 0 -5 13
16 0 -4 12
26 0 -12 14
25 1 -15 10
22 1 -12 10
21 1 -12 9
16 1 -6 10
17 1 -7 10
18 1 -7 11
;run;
proc glm data=ADS;
 class TRT01AN; /* 要因を指定 */
 model CHG = TRT01AN BASE / ss1 ss2 ss3 e solution;
 lsmeans TRT01AN / cl pdiff=control('0');
run;

f:id:cochineal19:20210223215441p:plain
f:id:cochineal19:20210223112602p:plain

プログラムコード

■ Rのコード

#-- 交互作用項が非有意かどうかで平行性を評価。
ANCOVA.0 <- lm(Y ~ X1 + C1 + X1*C1, data=ADS)
summary(ANCOVA.0)
car::Anova(ANCOVA.0)

#-- ANCOVA
ANCOVA.1 <- lm(CHG ~ BASE + TRT01AF, data=ADS)
(res <- summary(ANCOVA.1))
car::Anova(ANCOVA.1) #-- Type 2 平方和

■ SASのコード

proc glm data=ADS;
 class X1; /* 要因を指定 */
 model Y = X1 C1;
 lsmeans X1 / cl pdiff=control('XXX'); /* 調整平均 controlでレファレンスを指定*/
 estimate "X1 XXX vs. YYY" X1 -1 1; /* 対比を用いる場合 */
run;

■ Pythonのコード

整備中

雑談

水準毎の回帰直線が平行であることの評価方法
（交互作用項を含めたモデルを作り、交互作用項が非有意なら平行と解釈する方法）

本記事の架空データでの例：
① CHG=BASE + TRT01AN + BASE*TRT01AN を実行する。
② BASE*TRT01AN が非有意なら、CHG=BASE + TRT01AN のモデルでANCOVAを実行する。

参考

統計学（出版：東京図書）, 日本統計学会編

多変量解析実務講座テキスト, 実務教育研究所

https://www.sas.com/content/dam/SAS/ja_jp/doc/event/sas-user-groups/usergroups09-a-17h.pdf

http://nfunao.web.fc2.com/files/R-intro/R-stat-intro_06.pdf

★サイトマップ

cochineal19.hatenablog.com

2021-02-20

【統計】二元配置分散分析（Two-way ANOVA）

SAS 統計 R、R Studio Python

本記事では二元配置分散分析（Two-way ANOVA）について記載します。

【目次】

帰無仮説、対立仮設
計算式等
ノート
プログラムコード
雑談
参考

前回の一元配置分散分析は１要因に対する分析手法でした。
２要因の場合は、二元配置分散分析を用います。
※本記事の二元配置分散分析は、対応なしの被験者間実験計画（ABSタイプ）です。

帰無仮説、対立仮設

要因Aの主効果

帰無仮説 $H_{0}:$ 要因Aの水準間の主効果が同じ（平均が同じ）
対立仮設 $H_{1}:$ 要因Aの水準間のうち、少なくとも一つの水準で主効果に差がある（平均に差がある）

要因Bの主効果

帰無仮説 $H_{0}:$ 要因Bの水準間の主効果が同じ（平均が同じ）
対立仮設 $H_{1}:$ 要因Bの水準間のうち、少なくとも一つの水準で主効果に差がある（平均に差がある）

要因Aと要因Bの交互作用

帰無仮説 $H_{0}:$ 要因Aと要因Bに交互作用がある
対立仮設 $H_{1}:$ 要因Aと要因Bに交互作用がない

二元配置分散分析では、要因Aの水準間の平均の違い、要因Bの水準間の平均の違い、要因AとBの交互作用の有無を評価できます。

また、要因AとBに交互作用が無いと分かっている場合は交互作用を含めないモデルで分析可能です。交互作用が無いとは、例えば「要因A1の時の要因B1とB2の差」と「要因A2の時の要因B1とB2の差」に違いがないと言うことです。

計算式等

※本計算は各水準のサンプルサイズが等しいことを前提とします。また、対応なしの被験者間実験計画（ABSタイプ）です。

交互作用項を含めた分析
■ 平方和
$\quad S_{A} = \sum _{k=1} ^{a} n_{k} \cdot \left( \overline{y_{k}} - \bar{\bar{y}} \right)^{2}$
$\quad S_{B} = \sum _{m=1} ^{b} n_{m} \cdot \left( \overline{y_{m}} - \bar{\bar{y}} \right)^{2}$
$\quad S_{AB} = \sum _{k=1} ^{a} \sum _{m=1} ^{b} n_{km} \cdot \left( \overline{y_{km}} - \overline{y_{k}} - \overline{y_{m}} + \bar{\bar{y}} \right)^{2}$
$\quad S_{e} = \sum _{k=1} ^{a} \sum _{m=1} ^{b} \sum _{i=1} ^{n_{km}} \left( y_{km_{i}} - \overline{y_{km}} \right)^{2}$
$\quad S_{T} = \sum _{i=1} ^{n} \left( y_{i} - \bar{\bar{y}} \right)^{2}$

■ 分散分析表

変動因	平方和	自由度	平均平方	F値（F比）
水準間A	$S_{A}$	$f_{A}=a-1$	$V_{A}=\dfrac{S_{A}}{f_{A}}$	$F_{A}=\dfrac{V_{A}}{V_{e}}$
水準間B	$S_{B}$	$f_{B}=b-1$	$V_{B}=\dfrac{S_{B}}{f_{B}}$	$F_{B}=\dfrac{V_{B}}{V_{e}}$
交互作用AB	$S_{AB}$	$f_{AB}=a\cdot b-a-b+1$	$V_{AB}=\dfrac{S_{AB}}{f_{AB}}$	$F_{AB}=\dfrac{V_{AB}}{V_{e}}$
残差e	$S_{e}$	$f_{e}=n-a \cdot b$	$V_{e}=\dfrac{S_{e}}{f_{e}}$
全体T	$S_{T}$	$f_{T}=n-1$

※ k=要因Aの特定水準、m=要因Bの特定水準、i=特定の観測値、 $\overline{y_{k}}$ =要因Aの特定水準の平均値、 $\overline{y_{m}}$ =要因Bの特定水準の平均値、 $\bar{\bar{y}}$ =全体の平均値（総平均）、a=要因Aの水準数、b=要因Bの水準数

ここで平方和と自由度には次の関係があります。
$\quad S_{T} = S_{A} + S_{B} + S_{AB} + S_{e}$
$\quad f_{T} = f_{A} + f_{B} + f_{AB} + f_{e}$

交互作用項を含めない分析
■ 平方和　※ $S_{e}$ 以外同じ。交互作用項を残差項にプールしています。
$\quad S_{A} = \sum _{k=1} ^{a} n_{k} \cdot \left( \overline{y_{k}} - \bar{\bar{y}} \right)^{2}$
$\quad S_{B} = \sum _{m=1} ^{b} n_{m} \cdot \left( \overline{y_{m}} - \bar{\bar{y}} \right)^{2}$
$\quad S_{e} = \sum _{k=1} ^{a} \sum _{m=1} ^{b} \sum _{i=1} ^{n_{km}} \left( y_{km_{i}} - \overline{y_{k}} - \overline{y_{m}} + \bar{\bar{y}} \right)^{2}$
$\quad S_{T} = \sum _{i=1} ^{n} \left( y_{i} - \bar{\bar{y}} \right)^{2}$

■ 分散分析表　※ $S_{e}$ 以外同じ。

変動因	平方和	自由度	平均平方	F値（F比）
水準間A	$S_{A}$	$f_{A}=a-1$	$V_{A}=\dfrac{S_{A}}{f_{A}}$	$F_{A}=\dfrac{V_{A}}{V_{e}}$
水準間B	$S_{B}$	$f_{B}=b-1$	$V_{B}=\dfrac{S_{B}}{f_{B}}$	$F_{B}=\dfrac{V_{B}}{V_{e}}$
残差e	$S_{e}$	$f_{e}=n-a-b+1$	$V_{e}=\dfrac{S_{e}}{f_{e}}$
全体T	$S_{T}$	$f_{T}=n-1$

ここで平方和と自由度には次の関係があります。
$\quad S_{T} = S_{A} + S_{B} + S_{e}$
$\quad f_{T} = f_{A} + f_{B} + f_{e}$

また、質的変数をダミー変数（０と１の数量変数）にして重回帰分析を行う（正規方程式を解く）ことで、各水準の条件付き期待値を算出できます（数量化理論１類）。
こちらは一元配置分散分析の記事を参照ください。

ノート

簡単なデータで計算方法を確認します。
今回は薬剤（A薬、B薬、C薬の３水準）と前処置（有、無の２水準）の２要因において、評価指標Xの結果に差があるかを見ていきます（架空データです）。交互作用項を含めた分析を行います。なお、評価指標Xは特定のものではありませんが、低いほど良いとしておきます。

　A薬＋前処置無 {20,21,22,20,21}
　A薬＋前処置有 {18,16,17,16,17}
　B薬＋前処置無 {19,18,17,19,19}
　B薬＋前処置有 {19,18,18,17,17}
　C薬＋前処置無 {20,18,16,19,18}
　C薬＋前処置有 {15,14,17,16,18}

まず、全体・薬剤別・前処置別で平均を求め、総平均からのズレ（偏差平方和）を求めます。

薬剤	前処置	N	平均	要因:薬剤によるズレ（偏差平方和）と要因:前処置によるズレ（偏差平方和）
全体	全体	30	18.0	---
A薬	全体	10	18.8	$\left(18.8-18.0\right)^{2}\cdot 10=6.4$
B薬	全体	10	18.1	$\left(18.1-18.0\right)^{2}\cdot 10=0.1$
C薬	全体	10	17.1	$\left(17.1-18.0\right)^{2}\cdot 10=8.1$
全体	無	15	19.1	$\left(19.1-18.0\right)^{2}\cdot 10≒19.267$
全体	有	15	16.9	$\left(16.9-18.0\right)^{2}\cdot 10≒19.267$

同じく、薬剤×前処置別で平均を求め、総平均・各要因の平均からのズレ（交互作用に関する変動）を求めます。

薬剤	前処置	N	平均	総平均・各要因の平均からのズレ（交互作用に関する変動）
A薬	無	5	20.8	$\left(20.8-18.8-19.1+18.0\right)^{2}\cdot 5≒3.756$
A薬	有	5	16.8	$\left(16.8-18.8-16.9+18.0\right)^{2}\cdot 5≒3.756$
B薬	無	5	18.4	$\left(18.4-18.1-19.1+18.0\right)^{2}\cdot 5≒3.472$
B薬	有	5	17.8	$\left(17.8-18.1-16.9+18.0\right)^{2}\cdot 5≒3.472$
C薬	無	5	18.2	$\left(18.2-17.1-19.1+18.0\right)^{2}\cdot 5≒0.006$
C薬	有	5	16.0	$\left(16.0-17.1-16.9+18.0\right)^{2}\cdot 5≒0.006$

各表の最右列のズレを足し合わせて、各要因の水準間平方和と交互作用の平方和を計算します。
$\quad S_{A}=6.4+0.1+8.1=14.6$
$\quad S_{B}=19.267+19.267≒38.533$
$\quad S_{AB}=3.756+3.756+3.472+3.472+0.006+0.006=14.467$

また、今回は薬剤が３水準（A薬、B薬、C薬）、前処置が２水準（有、無）のため自由度は次のとおりです。
$\quad f_{A}=3-1=2$
$\quad f_{B}=2-1=1$
$\quad f_{AB}=3\times2-3-2+1=2$

次に、個々の観測値について、各水準の平均からのズレ（群内の偏差平方）と総平均からのズレ（全体の偏差平方）を求めます。
※ $S_{T} = S_{A} + S_{B} + S_{AB} + S_{e}$ なので $S_{T}$ か $S_{e}$ の１つ分かれば十分ですが、計算の理解のため両方求めます。

薬剤	前処置	評価指標	各水準の平均からのズレ（群内の偏差平方）	総平均からのズレ全体の偏差平方
A薬	無	20	$(20 - 20.8)^{2}=0.64$	$(20-18)^{2}=4$
A薬	無	21	$(21 - 20.8)^{2}=0.04$	$(21-18)^{2}=9$
A薬	無	22	$(22 - 20.8)^{2}=1.44$	$(22-18)^{2}=16$
A薬	無	20	$(20 - 20.8)^{2}=0.64$	$(20-18)^{2}=4$
A薬	無	21	$(21 - 20.8)^{2}=0.04$	$(21-18)^{2}=9$
A薬	有	18	$(18 - 16.8)^{2}=1.44$	$(18-18)^{2}=0$
A薬	有	16	$(16 - 16.8)^{2}=0.64$	$(16-18)^{2}=4$
A薬	有	17	$(17 - 16.8)^{2}=0.04$	$(17-18)^{2}=1$
A薬	有	16	$(16 - 16.8)^{2}=0.64$	$(16-18)^{2}=4$
A薬	有	17	$(17 - 16.8)^{2}=0.04$	$(17-18)^{2}=1$
B薬	無	19	$(19 - 18.4)^{2}=0.36$	$(19-18)^{2}=1$
B薬	無	18	$(18 - 18.4)^{2}=0.16$	$(18-18)^{2}=0$
B薬	無	17	$(17 - 18.4)^{2}=1.96$	$(17-18)^{2}=1$
B薬	無	19	$(19 - 18.4)^{2}=0.36$	$(19-18)^{2}=1$
B薬	無	19	$(19 - 18.4)^{2}=0.36$	$(19-18)^{2}=1$
B薬	有	19	$(19 - 17.8)^{2}=1.44$	$(19-18)^{2}=1$
B薬	有	18	$(18 - 17.8)^{2}=0.04$	$(18-18)^{2}=0$
B薬	有	18	$(18 - 17.8)^{2}=0.04$	$(18-18)^{2}=0$
B薬	有	17	$(17 - 17.8)^{2}=0.64$	$(17-18)^{2}=1$
B薬	有	17	$(17 - 17.8)^{2}=0.64$	$(17-18)^{2}=1$
C薬	無	20	$(20 - 18.2)^{2}=3.24$	$(20-18)^{2}=4$
C薬	無	18	$(18 - 18.2)^{2}=0.04$	$(18-18)^{2}=0$
C薬	無	16	$(16 - 18.2)^{2}=4.84$	$(16-18)^{2}=4$
C薬	無	19	$(19 - 18.2)^{2}=0.64$	$(19-18)^{2}=1$
C薬	無	18	$(18 - 18.2)^{2}=0.04$	$(18-18)^{2}=0$
C薬	有	15	$(15 - 16)^{2}=1$	$(15-18)^{2}=9$
C薬	有	14	$(14 - 16)^{2}=4$	$(14-18)^{2}=16$
C薬	有	17	$(17 - 16)^{2}=1$	$(17-18)^{2}=1$
C薬	有	16	$(16 - 16)^{2}=0$	$(16-18)^{2}=4$
C薬	有	18	$(18 - 16)^{2}=4$	$(18-18)^{2}=0$

このうち、各水準の平均からのズレ（群内の偏差平方）を足し合わせたものが残差平方和で、総平均からのズレ（全体の偏差平方）を足し合わせたものが全体平方和です。
$\quad S_{e}≒0.64 + 0.04 + ... + 0 + 4≒30.4$
$\quad S_{T}≒4+9+...+4+0≒98$

また、全体のN数が30、薬剤が３水準（A薬、B薬、C薬）、前処置が２水準（有、無）のため自由度は次のとおりです。
$\quad f_{e}=n-a\cdot b=30-3\cdot2=24$
$\quad f_{T}=n-1=30-1=29$

最後に、これら平方和・自由度を分散分析表にまとめて F値と p値を算出します。

変動因	平方和	自由度	平均平方	F値（F比）	p値
薬剤	$S_{A}=14.6$	$f_{A}=2$	$V_{A}=7.3$	$F_{A}≒5.7632$	$0.0090$
前処置	$S_{B}≒38.533$	$f_{B}=1$	$V_{B}≒38.533$	$F_{B}≒30.4211$	$0.0000$
交互作用	$S_{AB}≒14.467$	$f_{AB}=2$	$V_{AB}≒7.233$	$F_{AB}≒5.7105$	$0.0094$
残差e	$S_{e}=30.4$	$f_{e}=24$	$V_{e}≒1.267$
全体T	$S_{T}=98$	$f_{T}=29$

分散分析表のとおり、今回の架空データでは薬剤と前処置の交互作用が p<0.01 で有意であり、これら２要因の間には交互作用があるようです（何れかの組み合わせで相乗効果か相殺効果がある）。

なお、薬剤の主効果は p<0.01、前処置の主効果は p<0.001 でともに統計学的に有意でしたが、交互作用があるため参考程度にとどめておきます（独立していないため、薬剤と前処置を切り離して考えない方が良い）。

交互作用プロットを見ると、薬剤毎の前処置有無の平均の差が一定ではないため、交互作用が有意になったと考えられます。特にB薬は前処置に関する効果が薄いように見えます。

f:id:cochineal19:20210220141315p:plain

※変数名：TRT01AN=薬剤、水準：1=A薬、2=B薬、3=C薬
　変数名：PRCFL=前処置有無、水準：0=無、1=有
　変数名：AVAL=評価指標X

二元配置分散分析の結果の解釈と進め方

今回の架空データでは交互作用ありでしたが、二元配置分散分析の結果によって次のアクションに違いがあるため、備忘録がてらメモします。
※多重比較を事後検定と認識するのはどうだ、交互作用があったら単純主効果検定など色々な見解を見かけましたが、こちらのサイト様が分かりやすかったので、参考にさせていただきました。

１．交互作用が認められる
　⇒　薬剤と前処置が互いに影響を与えていると結論。
　　　必要であれば薬剤×前処置の６つの組み合わせで多重比較。
　　　各要因を切り離して考えない。

２．交互作用が認められない＆薬剤の主効果が有意
　⇒　薬剤に有意な変動があったと結論。
　　　薬剤が３水準なので、必要であれば多重比較。

３．交互作用が認められない＆前処置の主効果が有意

　⇒　前処置に有意な変動があったと結論。
　　　２水準なので多重比較する必要なし。

４．交互作用が認められない＆薬剤・前処置ともに主効果が有意

　⇒　薬剤・前処置にともに有意な変動があったと結論。
　　　加えて、薬剤×前処置に交互作用がなく互いに影響を与えていないと判断。
　　　必要であれば、薬剤が３水準なので、薬剤について多重比較。
　　　前処置は２水準なので完結。

Rでの実行：

> library(car) 
> ADS <- data.frame(
+   TRT01AN=c(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2)
+   ,PRCFL=c(0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,1,1,1,1,1)
+   ,AVAL=c(20,21,22,20,21,18,16,17,16,17,19,18,17,19,19,19,18,18,17,17,20,18,16,19,18,15,14,17,16,18))
> ADS$前処置 <- factor(ifelse(ADS$PRCFL==0,"無","有"),levels=c("無","有"))
> ADS$投与群 <- factor(ifelse(ADS$TRT01AN==0,"A薬",ifelse(ADS$TRT01AN==1,"B薬","C薬")),levels=c("A薬","B薬","C薬"))
> fit1 <- lm(AVAL ~ 投与群 + 前処置 + 投与群:前処置, data=ADS)
> car::Anova(fit1, Type="II")
Anova Table (Type II tests)

Response: AVAL
              Sum Sq Df F value    Pr(>F)    
投与群        14.600  2  5.7632  0.009035 ** 
前処置        38.533  1 30.4211 1.134e-05 ***
投与群:前処置 14.467  2  5.7105  0.009363 ** 
Residuals     30.400 24                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SASでの実行：

data ADS;
 do TRT01AN = 1; do PRCFL = 0; do AVAL = 20,21,22,20,21; output; end; end; end;
 do TRT01AN = 1; do PRCFL = 1; do AVAL = 18,16,17,16,17; output; end; end; end;
 do TRT01AN = 2; do PRCFL = 0; do AVAL = 19,18,17,19,19; output; end; end; end;
 do TRT01AN = 2; do PRCFL = 1; do AVAL = 19,18,18,17,17; output; end; end; end;
 do TRT01AN = 3; do PRCFL = 0; do AVAL = 20,18,16,19,18; output; end; end; end;
 do TRT01AN = 3; do PRCFL = 1; do AVAL = 15,14,17,16,18; output; end; end; end;
run;
proc glm data=ADS;
 class TRT01AN PRCFL;
 model AVAL = TRT01AN PRCFL TRT01AN*PRCFL / ss1 ss2 ss3 e1 e2 e3 solution; 
run;

f:id:cochineal19:20210220140739p:plain

f:id:cochineal19:20210220140807p:plain

プログラムコード

Rのコード

library(car)
fit1 <- lm(AVAL ~ FACTOR1 + FACTOR2 + FACTOR1:FACTOR2, data=ADS) #-- FACTOR1:FACTOR2は交互作用項
car::Anova(fit1, Type="II") #-- Type2平方和

SASのコード

proc glm data=ADS1 outstat=OUTDS;
  class FACTOR1; /*class指定必須*/
  model AVAL = FACTOR1 FACTOR2 FACTOR1*FACTOR2 / ss1 ss2 ss3 e1 e2 e3 solution;
run;
/* option: ss1=Type1平方和, ss2=Type2平方和, ss3=Type3平方和 */
/* option: e#= Type#平方和のEstimable Functions（推定可能関数）を図示してくれる */
/* option: solution=回帰モデルの推定値を出力してくれる */

Pythonのコード

整備中

雑談

交互作用有無のイメージです。

f:id:cochineal19:20210220145954p:plain

　⇒　前処置の効果がA～C薬どこでも一定であれば交互作用なし。

f:id:cochineal19:20210220152536p:plain

　⇒　前処置の効果が一定でなければ交互作用あり。
　　　薬剤と前処置の組み合わせでより効果が出たり（相乗効果）、
　　　効果が薄まったり（相殺効果）する。

参考

統計学基礎（出版：東京図書）, 日本統計学会編

統計学（出版：東京図書）, 日本統計学会編

改訂増補版：統計検定を理解せずに使っている人のためにIII

★サイトマップ

cochineal19.hatenablog.com

2021-01-27

【R】ggplot：svg画像の作成、svg画像へのリンク埋め込み

R、R Studio

備忘録がてら、ggplot で svg形式ファイルを作る方法とリンクの埋め込み方法を書きます。

svg形式ファイルとは？
ggplotでの出力方法
svg形式ファイルへのリンク埋め込み
お試し画像

svg形式ファイルとは？

svg形式ファイルは xmlをベースにしたベクター画像です。
画像情報を座標（テキスト）で持っているため、拡大しても荒くならない、容量が軽いなどの特徴があります。
一方、png や jpeg などの形式はラスター画像と言われます。

ggplotでの出力方法

ggplotはsvg形式ファイルに対応していますが、svg形式を出力するために svglite パッケージが必要です。

 library(svglite)

ggplotでの出力方法は次のとおり。

ggsave(filename="./output_file.svg", plot=ggplot_data
       ,width=5, height=5, dpi=300, units="in", device="svg")

ggsave(filename="./output_file.svg", plot=ggplot_data) #-- ファイル名のみでも良い。

svg形式ファイルへのリンク埋め込み

svg形式にすることで画像内にリンクを埋め込むことができます。
annotateで埋め込んだテキストにリンクを貼って別の画像を呼び出すなんてこともできます。

理解のため、試しにやってみます。
irisデータで適当に作った散布図に文字列「テスト」のアノテーションを挿入しました。

library(tidyverse)
library(svglite)
g1 <- ggplot(data=iris, mapping=aes(x=Sepal.Length, y=Sepal.Width)) +
  geom_point() + 
  annotate("text",x=7.5,y=2.5,label="テスト",size=8)
ggsave(filename="./ggplot_test.svg", plot=g1
       ,width=5, height=4, dpi=500, units="in", device='svg')

このsvg形式ファイルをテキストファイルで開いて下にスクロールしていくと、ggplotでアノテーションした文字列「テスト」が書かれた部分があります。

＜置換前のコード（svgファイル内）＞

<g clip-path='url(#cpMzUuNTd8MzU0LjUyfDI1Ni41MXw1LjQ4)'><text x='258.64' y='181.11' style='font-size: 11.04px; font-family: Arial;' textLength='17.79px' lengthAdjust='spacingAndGlyphs'>テスト</text></g>

この文字列「テスト」をaタグに変えればリンク埋め込み済みの svg 画像の完成です。
注意点として、XLinkを使うため、href要素を xlink:href と書きます。
また、target="_blank" で別タブに飛ばすようにします。

＜置換後のコード（svgファイル内）＞

<g clip-path='url(#cpMzUuNTd8MzU0LjUyfDI1Ni41MXw1LjQ4)'><text x='258.64' y='181.11' style='font-size: 11.04px; font-family: Arial;' textLength='17.79px' lengthAdjust='spacingAndGlyphs'><a xlink:href="パス/ファイル名.svg" target="_blank">テスト</a></text></g>

お試し画像

先ほどのサンプルコードで作ったirisデータの散布図です。
画像内の文字列「テスト」にリンクを仕込んでいます。
クリックすると、本ブログの統計サイトマップに飛ぶはず。

以上です。

2021-01-25

【統計】ロジスティック回帰分析（ホスマー・レメショウ検定、ROC、AUC等）

R、R Studio Python SAS 統計

本記事では、ロジスティック回帰の評価について記載します。

【目次】

計算式等
ノート
プログラムコード
雑談
参考

前回記事で理論部分を取り扱いましたが、続きでモデルの評価方法を取り扱います。

cochineal19.hatenablog.com

計算式等

モデル評価には、大きく分けて「Discrimination（判別能力）」と「Calibration（較正）」があります。
「Discrimination（判別能力）」では、モデルが正しくイベントの発生有無を判別できるかの予測性能を評価します。
「Calibration（較正）」では、サブグループ毎に実測値と予測値のズレ具合を評価し、部分的に予測精度の悪さがないかを評価します。

Discrimination（判別能力）

① 真陽性率（感度）、真陰性率（特異度）等

実測と予測の関係を混同行列（Confusion Matrix）に表すと次のとおりです。

		実測
		(+)	(-)
予測	(+)	真陽性 Ture Positive (TP)	偽陽性 False Positive (FP)
予測	(-)	偽陰性 False Negative (FN)	真陰性 True Negative (TN)

この混同行列を参考として、各種評価指標を次のように計算します。

項目	式	類語
真陽性率 (True Positive Rate)	TP / (TP + FN)	感度 (Sensitivity)、検出率 (Recall)
真陰性率 (True Negative Rate)	TN / (FP + TN)	特異度 (Specificity)
偽陽性率 (False Positive Rate)	FP / (FP + TN)	1 - 特異度 (Specificity)
偽陰性率 (False Negative Rate)	FN / (TP + FN)	1 - 感度 (Sensitivity)
陽性適中率 (Positive Predictive Value)	TP / (TP + FP)	適合率 (Precision)
陰性適中率 (Negative Predictive Value)	TN / (FN + TN)	ー
偽発見率 (False Discovery Rate)	FP / (TP + FP)	1 - 適合率 (Precision)
正診率 (Accuracy)	(TP + TN) / N	ー

類語の感度（検出率）、特異度、適合率、検出率の方が聞きなじみがあるかもしれません。どの指標を使うかは目的に応じます。

・感度（検出率）：実際に陽性だったもののうち、どれだけ陽性と言えたか。
・特異度：実際に陰性だったもののうち、どれだけ陰性と言えたか。
・適合率：陽性と予測したもののうち、本当に陽性だったのはどのくらいか。
・正診度：予測全体のうち、本当に陽性または陰性だったのはどのくらいか。

なお、例えば、めったに起きないイベントの有無を予測したい場合、無と答えていれば大体正解できてしまうため予測性能が高いモデルのように見えてしまいます。
このような不均衡データに対応する指標として、マシューズ相関係数（MCC）やF1スコアなどがあります。F1スコアは検出率と適合率の調和平均です。

$\quad MCC=\dfrac{TP\times TN-FP\times FN}{\sqrt{\left( TP+FP\right) \left( TP+FN\right) \left( TN+FP\right) \left( TN+FN\right) }}$

$\quad F1\ Score=2\times\dfrac{Recall \times Precision}{Recall + Precision}=2\times\dfrac{\{TP / \left(TP + FN\right)\}\times \{TP / \left(TP + FP\right)\} }{ \{TP / \left(TP + FN\right)\} + \{TP / \left(TP + FP\right)\} }$

② ROC曲線（Receiver Operating Characteristic Curve）、
　 AUC（Area Under the Curve）

ROC曲線は、縦軸に真陽性率（感度）、横軸に偽陽性率（１－特異度）をプロットした曲線です。横軸を真陰性率（特異度）とすることもあります。
AUCは、ROC曲線下の面積です。０～１の範囲を取り、一つの基準として " 0.7 以上で可、0.8 以上で良、0.9 以上で優" とされます。

f:id:cochineal19:20210130225421p:plain

Calibration（較正）

① ホスマー・レメショウ検定（Hosmer and Lemeshow goodness of fit (GOF) test）

ホスマー・レメショウ検定は、データセットをk個に分割してサブグループ（リスクグループ）を作成します。このサブグループ毎に "実測度数 $O_{k}$ と期待度数 $E_{k}$ のズレ" を算出し、足し合わせた $\chi _{HL}^{2}$ 統計量を用いてモデルの適合度を評価する手法です。
※ 死亡有無を予測するモデルを例にすると、死亡率０～９％、10％台、20％台… とリスクグループを分割することで各リスクグループでの予実の差（ズレ）を評価できます。これにより、特定のリスクグループに予測精度の悪さがないかを評価できます。特に目立ったズレがなければ全体的にモデルが適合していると言えます。

$\qquad \chi _{HL}^{2}=\sum \dfrac{O_{k}-E_{k}}{N_{k}\widehat{p} _{k}\left( 1-\widehat{p}_{k}\right) }=\sum \dfrac{O_{k}-N_{k}\widehat{p} _{k}}{N_{k}\widehat{p} _{k}\left( 1-\widehat{p}_{k}\right) }$

　　※ $\widehat{p}$ は予測値 $\widehat{y}$ の平均値。
　　※ グループ数は明確な基準はない。５～１０くらい？

$\chi _{HL}^{2}$ は自由度 $k-2$ の $\chi ^{2}$ 分布に近似することができ、p>0.05 の時、帰無仮説を棄却せず、モデルが適合していると評価します。
　帰無仮説 $H_{0}: p=\widehat{p}$ （モデルがどのサブグループにも適合している）
　対立仮説 $H_{1}: p≠\widehat{p}$ （モデルがいずれかのサブグループに適合していない）

② Calibration plot

予測値 $\widehat{p}_{k}$ と実測値 $p_{k}$ をプロットしたグラフです。
上のホスマー・レメショウ検定の分割したデータセットを使って作ることができます。

---------
その他にも、McFaddenの疑似決定係数、Brier Score、Log Lossなどの指標があります。

ノート

ロジスティック回帰の理論部分で使ったサンプルデータを用います。
（SASのサンプルデータ作成は前回記事を参照ください）
※機械学習では訓練データ：バリデーションデータ：テストデータなどに分割したりしますが、ここでは取り扱いません（また別のタイミングで）。

ads <- data.frame(
  x1=c(rep(5,20),rep(5,22),rep(10,25),rep(10,30),rep(15,30),rep(15,30))
  ,x2=c(rep(0,20),rep(1,22),rep(0,25),rep(1,30),rep(0,30),rep(1,30))
  ,y=c(rep(0,19),rep(1,1),rep(0,16),rep(1,6),rep(0,16),rep(1,9),rep(0,14),rep(1,16),rep(0,11),rep(1,19),rep(0,2),rep(1,28))
)
fit1 <- glm(y ~ x1 + x2, family=binomial, data=ads)
ads$yhat <- fit1$fitted.values

------------------------------------

Discrimination（判別能力）

ROC曲線を描画し、AUCを算出します。

Rでの実行：

library(tidyverse)
library(pROC)
(roc1 <- pROC::roc(response=fit1$y, predictor=fit1$fitted, ci=T))
roc1.res1 <- (pROC::coords(roc1, x="best", best.method="youden", ret=c("threshold","tpr","fpr","tnr","fnr","npv","ppv"), transpose=F))
roc1.res2 <- (pROC::coords(roc1, x="best", best.method="closest.topleft", ret=c("threshold","tpr","fpr","tnr","fnr","npv","ppv"), transpose=F))
AUC1 <- round(roc1$ci, 4)

ggroc(list(Model1=roc1), legacy.axes=T, size=1.5, color="darkgreen", alpha=.7) +
  geom_abline(intercept=0, slope=1, linetype="dashed") +
  annotate("text", x=1-roc1.res1$tnr, y=roc1.res1$tpr, hjust=0.5, vjust=0.5, size=10, label="×", color="blue") +
  annotate("text", x=1-roc1.res2$tnr, y=roc1.res2$tpr, hjust=0.5, vjust=0.5, size=10, label="〇", color="red") +
  annotate("text", x=0.75, y=0.1, size=6
           ,label=paste0("Area Under the Curve ( 95%CI ) :\n",AUC1[2]," ( ",AUC1[1],", ",AUC1[3]," )")) +
  labs(title="ROC曲線")

f:id:cochineal19:20210130230434p:plain

SASでの実行：

ods output LackFitChiSq=HOSLEM1;
proc logistic data=ads descending;
  model y(event='1') = x1 x2 / outroc=ROC1 lackfit;
  output out=OutDS;
run;

f:id:cochineal19:20210130230629p:plain

Rのグラフでは最適カットオフ値を載せています。
Youden Indexが青色の×印、closest.topleft が赤色の〇印です。
Youden Indexは45度線から一番離れているROC曲線のポイント、closest.topleftは左上の角（AUC=1.0の位置）から一番近いROC曲線のポイントをそれぞれ最適カットオフ値とします。今回のデータでは、両方法とも同じ最適カットオフ値（threshold=0.4497）でした。

今回のデータでは、AUC=0.8085 [ 0.7435, 0.8735] で0.8を超える良好な結果でした。
また、最適カットオフ値での混同行列（Confusion Matrix）と各評価指標は次のとおりです。

		実測
		(+)	(-)
予測	(+)	６３	２７
予測	(-)	１６	５１

・感度（検出率）=63/(63+13)=0.7975
・特異度 =51/(27+51)=0.6538
・適合率 =63/(63+27)=0.7000
・正診率 =(63+51)/(63+27+16+51)=0.7261
・F1スコア=2×0.7975×0.7/(0.7975+0.7)=0.7456
・MCC =(63×51-27×16)/sqrt{(63+27)(63+16)(51+27)(51+16)}=0.4562

------------------------------------

Calibration（較正）

ホスマー・レメショウ検定をRとSASで行います。

Rでの実行：（※グループ数はSASの実行結果に合わせて 6 にしています）

> library(tidyverse)
> library(ResourceSelection)
> (HL1 <- ResourceSelection::hoslem.test(fit1$y, fit1$fitted, g=6))
	Hosmer and Lemeshow goodness of fit (GOF) test
data:  fit1$y, fit1$fitted
X-squared = 2.1597, df = 4, p-value = 0.7064

> cbind(HL1$observed, HL1$expected)
               y0 y1     yhat0     yhat1
[0.0763,0.237] 35  7 35.261211  6.738789
(0.237,0.292]  16  9 17.705136  7.294864
(0.292,0.608]  14 16 11.772441 18.227559
(0.608,0.673]  11 19  9.821233 20.178767
(0.673,0.885]   2 28  3.439979 26.560021

SASでの実行は先ほどのコードから。

f:id:cochineal19:20210130231534p:plain

残念ながら結果は一致せず、どうやら bin の設定の仕方が異なるようです。
ということで、Rのコードを自作します。

> #-----------------------------------------------------------
> #-- 自作 Hosmer and Lemeshow test, and Calibration Plot
> #-----------------------------------------------------------
> grp1 <- 6 #-- グループ数
> sep1 <- unique(c(0,quantile(fit1$fitted.values, probs = seq(0, 1, 1/grp1)))) #-区間
> out1 <- data.frame(y=fit1$y
+                    ,yhat=fit1$fitted.values
+                    ,ycut=cut(fit1$fitted.values, breaks=sep1, include.lowest=T))
> out2 <- out1 %>% 
+   dplyr::group_by(ycut) %>%
+   dplyr::summarise(Act=sum(y)/n()
+                    ,Pred=mean(yhat)
+                    ,N=n()
+                    ,O=sum(y)
+                    ,E=n()*mean(yhat)
+                    ,ChiSq=(sum(y) - n()*mean(yhat))^2 / (n()*mean(yhat)*(1-mean(yhat))))
> out3 <- data.frame(Chi_Square = sum(out2$ChiSq)
+                    ,P_value = 1 - pchisq(sum(out2$ChiSq), grp1 - 2))
> print(out2)
  ycut             Act   Pred     N     O     E ChiSq
                  
1 [0,0.0763]     0.05  0.0763    20     1  1.53 0.197
2 (0.0763,0.237] 0.273 0.237     22     6  5.21 0.156
3 (0.237,0.292]  0.36  0.292     25     9  7.29 0.563
4 (0.292,0.608]  0.533 0.608     30    16 18.2  0.694
5 (0.608,0.673]  0.633 0.673     30    19 20.2  0.210
6 (0.673,0.885]  0.933 0.885     30    28 26.6  0.681
> print(out3$Chi_Square)
[1] 2.500154
> print(out3$P_value)
[1] 0.6446081

自作プログラムではSASの結果と一致しました。
今回は $p≒0.645$ で帰無仮説は却下せず、モデルが適合していると評価します。

Rでキャリブレーションプロットも描いてみます。
先ほどのホスマー・レメショウ検定で使用したグループ毎の予測値と実測値をプロットします。45° 線付近にプロットがあり、大きなズレはなさそうです。

ggplot(data=out2) +
  geom_abline(intercept=0, slope=1, linetype="dashed") +
  scale_x_continuous(limits=c(0,1), breaks=seq(0,1,0.25)) +
  scale_y_continuous(limits=c(0,1), breaks=seq(0,1,0.25)) +
  geom_point(mapping=aes(x=Pred, y=Act)) +
  labs(title="キャリブレーションプロット", x="予測値", y="実測値")

f:id:cochineal19:20210130232214p:plain

プログラムコード

■ Rのコード

#-----------------------------------------------------------------------------------
# Discrimination（判別能力）
# ROC、AUC、感度、特異度
#-----------------------------------------------------------------------------------
library(tidyverse)
library(pROC)
(roc1 <- pROC::roc(response=fit1$y, predictor=fit1$fitted, ci=T)) #-- glmの結果を投入
(roc1.res1 <- (pROC::coords(roc1, x="best", best.method="youden", ret=c("threshold","tpr","fpr","tnr","fnr","npv","ppv"), transpose=F)))
(roc1.res2 <- (pROC::coords(roc1, x="best", best.method="closest.topleft", ret=c("threshold","tpr","fpr","tnr","fnr","npv","ppv"), transpose=F)))
(AUC1 <- round(roc1$ci, 4))

ggroc(list(Model1=roc1), legacy.axes=T, size=1.5) +
  geom_abline(intercept=0, slope=1, linetype="dashed") +
  annotate("text", x=1-roc1.res1$tnr, y=roc1.res1$tpr, hjust=0.5, vjust=0.5, size=10, label="〇", color="blue") +
  annotate("text", x=1-roc1.res2$tnr, y=roc1.res2$tpr, hjust=0.5, vjust=0.5, size=10, label="×", color="blue") +
  annotate("text", x=1-roc1.res1$tnr, y=roc1.res1$tpr-0.05, hjust=0, vjust=0.5, size=5, label="youden",) +
  annotate("text", x=1-roc1.res2$tnr, y=roc1.res2$tpr-0.05, hjust=0, vjust=0.5, size=5, label="closest.topleft") +
  annotate("text", x=0.8, y=0.1, size=4
           ,label=paste0("Area Under the Curve ( 95%CI ) :\n",AUC1[2]," ( ",AUC1[1],", ",AUC1[3]," )")) +
  labs(title="ROC曲線")

#-----------------------------------------------------------------------------------
# Calibration（較正）
# Hosmer and Lemeshow goodness of fit (GOF) test
# Calibration plot
#-----------------------------------------------------------------------------------
grp1 <- 6 #-- グループ数
sep1 <- unique(c(0,quantile(fit1$fitted.values, probs = seq(0, 1, 1/grp1)))) #-- 区間
out1 <- data.frame(y=fit1$y                  #-- glmの結果を投入（実測値）
                   ,yhat=fit1$fitted.values  #-- glmの結果を投入（予測値）
                   ,ycut=cut(fit1$fitted.values, breaks=sep1, include.lowest=T) )
out2 <- out1 %>% 
  dplyr::group_by(ycut) %>%
  dplyr::summarise(Act=sum(y)/n()
                   ,Pred=mean(yhat)
                   ,N=n()
                   ,O=sum(y)
                   ,E=n()*mean(yhat)
                   ,ChiSq=(sum(y) - n()*mean(yhat))^2 / (n()*mean(yhat)*(1-mean(yhat))) )
out3 <- data.frame(Chi_Square = sum(out2$ChiSq)
                   ,P_value = 1 - pchisq(sum(out2$ChiSq), grp1 - 2))
print(out2)
print(out3$Chi_Square)
print(out3$P_value)

ggplot(data=out2) +
  geom_abline(intercept=0, slope=1, linetype="dashed") +
  scale_x_continuous(limits=c(0,1), breaks=seq(0,1,0.25)) +
  scale_y_continuous(limits=c(0,1), breaks=seq(0,1,0.25)) +
  geom_point(mapping=aes(x=Pred, y=Act)) +
  labs(title="キャリブレーションプロット", x="予測値", y="実測値")

■ SASのコード

ods output LackFitChiSq=HOSLEM1;
proc logistic data=ads descending;
  model y(event='1') = x1 x2 / outroc=ROC1 lackfit;
  output out=OutDS;
run;

■ Pythonのコード

整備中

雑談

何かあれば追記します。

参考

https://www.ism.ac.jp/~noma/2018-12-07%20JBS%20Seminar%20Kyoto.pdf

★サイトマップ

cochineal19.hatenablog.com