2023-06-12

【統計】割合の信頼区間（WaldとClopper-pearson）

割合の信頼区間について

【目次】

計算式
計算例
コード
被覆確率 (Coverage Probability)
参考

計算式

Wald (正規近似を用いた方法)
$\quad \widehat{p}\pm z\sqrt{\dfrac{\widehat{p}\left( 1-\widehat{p}\right) }{n}}$

二項分布 B (n, p) は n が十分大きいとき、平均 np , 分散 np(1-p) の正規分布に近づく（ラプラスの定理）。一般的に、np>5 かつ n(1-p)>5 の場合を n が十分に大きいと考えるようである。

n が十分に大きくない場合は近似が悪い。また、０や１の両端も近似が悪いので不均衡データ（Imbalanced Data）に注意が必要となる。

Clopper-pearson (F分布を用いた方法)
$\quad Lower = \Biggl( 1+ \dfrac{n-x+1}{(x)F\bigl( 1-\frac{\alpha}{2}; \ 2(x) , \ 2(n-x+1) \bigr) } \Biggr)^{-1}$
$\quad Upper = \Biggl( 1+ \dfrac{n-x}{(x+1)F\bigl( \frac{\alpha}{2}; \ 2(x+1) , \ 2(n-x) \bigr) } \Biggr)^{-1}$

二項分布（離散確率分布）をベータ分布（連続値型確率分布）に拡張し、F分布で表したものを利用した信頼区間。

計算例

Wald
=D4-NORM.S.INV(1-G4/2)*SQRT(D4*(1-D4)/C4)
=D4+NORM.S.INV(1-G4/2)*SQRT(D4*(1-D4)/C4)

Clopper-pearson
=(1+(C4-B4+1)/(B4*F.INV.RT(1-G5*0.5, 2*B4, 2*(C4-B4+1))))^-1
=(1+(C4-B4)/((B4+1)*F.INV.RT(G5*0.5, 2*(B4+1), 2*(C4-B4))))^-1

コード

library(DescTools)
DescTools::BinomCI(x, n, conf.level=.95,method=c("wald"))
DescTools::BinomCI(x, n, conf.level=.95,method=c("clopper-pearson"))

proc freq data=ADS;
    tables YN / binomial;
    weight CNT / zeros;
run;

被覆確率 (Coverage Probability)

被覆確率（信頼区間が真の値を含む確率）を調べてみる。
p = ０～１を 0.01 刻みで、二項分布 B(n, p) に従う乱数を100個ずつ生成（PCのスペックがないのでこのくらいで）。
n = 10, 100, 1000 の３種類を試してみる。

被覆確率 (Coverage Probability) をグラフ化し、95％のところに横線を引いている。
n が小さいとき、Waldは被覆確率95％を大きく下回る。両端の不均衡データにあたる部分でその傾向が大きく見られる。
Clopper-pearsonは n が小さくても95%から大きく下回ることはなさそう。

参考

https://www.lexjansen.com/phuse/2013/sp/SP05.pdf
https://www.sas.com/content/dam/SAS/ja_jp/doc/event/sas-user-groups/usergroups14-d-05.pdf
統計学入門−第3章

2023-06-10

【統計】割合の差の信頼区間（リスク差の信頼区間）

R、R Studio SAS 統計

SASとRでの使い方メモ。

	R DescTools::BinomDiffCI method = TYPE	SAS PROC FREQ RISKDIFF (CL = TYPE)
Agresti-Caffo	"ac"	AC
Brown, Li's Jeffreys	"blj"	-
Farrington and Manning	-	FM
Haldane	"hal"	-
Hauck-Anderson	"ha"	HA
Jeffreys-Perks	"jp"	-
Mee	"mee"	MN (CORRECT= NO)
Miettinen and Nurminen	"mn"	MN
Newcombe (Corrected)	"scorecc"	NEWCOMBE (CORRECT) \| SCORE (CORRECT) \| WILSON (CORRECT)
Newcombe	"score"	NEWCOMBE \| SCORE \| WILSON
Wald (Corrected)	"waldcc"	WALD (CORRECT)
Wald	"wald"	WALD

Rでの実行

c00 <- 50; c01 <- 75; N0 <- c00+c01;
c10 <- 40; c11 <- 90; N1 <- c10+c11;
 
library(DescTools)
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="ac")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="ha")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="mee")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="mn")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="score")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="scorecc")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="wald")

SASでの実行

data ADS;
    length FACTOR YN CNT 8.;
    input FACTOR YN CNT @@;
    cards;
0 0 50  0 1 75  1 0 40  1 1 90
;
run;
  
proc freq data=ADS;
    table FACTOR * YN / riskdiff(cl=(wald ha fm ac mn newcombe));
    table FACTOR * YN / riskdiff(cl=(wald(correct) mn(correct=no) newcombe(correct)));
    weight CNT / zeros;
run;

参考

https://cran.r-project.org/web/packages/DescTools/DescTools.pdf
https://www.sas.com/content/dam/SAS/ja_jp/doc/event/sas-user-groups/usergroups14-d-08.pdf
https://www.sas.com/content/dam/SAS/ja_jp/doc/event/sas-user-groups/usergroups2015-b-01.pdf
https://www.lexjansen.com/wuss/2016/127_Final_Paper_PDF.pdf
https://www.seikei.ac.jp/university/rikou/pdf/V0410202.pdf
https://www.ms.uky.edu/~mai/sta635/FagerlandLydersenLaake2011---RecommendedCIsForTwoIndependent....pdf
https://onlinelibrary.wiley.com/doi/epdf/10.1002/pst.2162

2023-06-05

【統計】ロジスティック回帰のFisher Scoringによる最適化

R、R Studio SAS 統計機械学習

【目次】

ロジスティック回帰の特徴
ロジスティック回帰の尤度・誤差関数
ロジスティック回帰のFisher Scoringによる最適化
手計算
参考

以下でロジスティック回帰に触れていたがFisher Scoringにて再び。
【統計】ロジスティック回帰分析 - こちにぃるの日記
 【統計】ロジスティック回帰分析（ホスマー・レメショウ検定、ROC、AUC等） - こちにぃるの日記
 【統計】一般化線形モデル（GLM） - こちにぃるの日記
 【Python】Kerasでロジスティック回帰 - こちにぃるの日記
 【Kaggle】タイタニックの分析をやってみる - こちにぃるの日記

ロジスティック回帰の特徴

ロジスティック回帰は一般化線形モデルの一つ。

目的変数	二値の離散値
分布族	二項分布
リンク関数	ロジットリンク関数
逆リンク関数	ロジスティック関数（シグモイド関数）

・誤差構造が二項分布に従うと仮定。
・ロジット関数で実測値 (０～１) をロジット値 $(-∞～∞)$ に変換し、線形モデルを作成。

$\qquad logit(θ_{i})=ln (\dfrac{θ_{i}}{1-θ_{i}})=X_{i}' \beta$

・ロジスティック関数で線形予測子から予測値（０～１）を得る。

$\qquad θ_{i} = logistic( X_{i}' \beta )=\dfrac{1}{1+exp( -X_{i}' \beta )}=\dfrac{exp( X_{i}' \beta )}{1+exp( X_{i}' \beta )}$

・二値データに対する通常の回帰モデルでは、予測値が０～１の範囲を超えてしまい解釈が難しいが、ロジスティック回帰の予測値は０～１の範囲に収まり、解釈しやすい。

（青線：ロジスティック回帰、緑線：単回帰分析）

ロジスティック回帰の尤度・誤差関数

尤度

$\quad L \left( \beta \right) =\Pi \{ θ_{i}^{\ y_{i}} \times \left( 1-θ_{i} \right)^{ \left( 1 - y_{i} \right) } \}$

対数尤度

$\quad ln\ L\left( \beta\right)=\sum \{ y_{i} \cdot ln\left(θ_{i}\right) + \left(1-y_{i} \right) \cdot ln \left( 1-θ_{i} \right) \}$

交差エントロピー誤差（分類問題の損失関数）

$\quad E\left(\beta\right)=-\sum \{ y_{i} \cdot ln\left(θ_{i}\right) + \left(1-y_{i} \right) \cdot ln \left( 1-θ_{i} \right) \} =0$

下図の通り、同一式で Y=１ or ０の実測値と予測値の差（誤差）が求まる。

ロジスティック回帰のFisher Scoringによる最適化

勾配

$\quad \nabla E( \beta )=\dfrac{\partial E( \beta) }{\partial θ}\dfrac{\partial θ}{\partial \widehat{y}}\dfrac{\partial \widehat{y}}{\partial \beta }$

$\qquad \qquad =-\sum \{ y_{i} \dfrac{1}{θ_{i}} θ_{i}( 1-θ_{i}) x_{i} + (1-y_{i}) \dfrac{1}{1-θ_{i}} ( -θ_{i}( 1-θ_{i} ) ) x_{i} \}$

$\qquad \qquad =-\sum \{ y_{i} ( 1-θ_{i}) + (1-y_{i}) ( -θ_{i} )\} x_{i}$

$\qquad \qquad =-\sum \{ y_{i} - y_{i} θ_{i} - θ_{i} + y_{i} θ_{i}\} x_{i}$

$\qquad \qquad =-\sum \{ y_{i} - θ_{i} \} x_{i}$

$\qquad \qquad =\sum ( θ_{i} - y_{i}) x_{i}=0$

$\quad \nabla E( \beta )= X' (θ-Y)$

Fisher情報行列

$\quad I=-E(\nabla^{2} E( \beta ))$

$\quad \nabla^{2} E( \beta ) =-\sum ^{n}_{i=1}\{x_{ik}' \ \theta _{i}(1-\theta _{i}) \ x_{ik})\}$

$\quad \nabla^{2} E( \beta ) =-X'WX$

更新式（Fisher Scoring）
Hessian行列を期待値に置き換えたFisher Scoring法。

$\quad \beta^{(t+1)}=\beta^{(t)}-(I^{(t)})^{-1}\nabla E( \beta )^{(t)}$

$\qquad\quad\ = \beta^{(t)} - (X'WX)^{-1} X'(θ-Y)$

手計算

Kaggleチュートリアルのタイタニックデータを用いる。
簡単に以下のモデルとする。

目的変数	Survived（生存有無）	０= 死亡, １= 生存
説明変数	Pclass（チケットのクラス）	１= 1st, ２= 2nd, ３= 3rd ※今回は簡単に０= 1st or 2nd, １= 3rd
	Sex（性別）	０= Female, １= Male

Rでの実行結果

Fisher Scoring iterations は 4回。

titanic$Event   <- 1 - titanic$Survived
titanic$Pclass3 <- if_else(titanic$Pclass=="3",1,0)
titanic$Sex2    <- if_else(titanic$Sex=="male",1,0)
  
Call:
glm(formula = Event ~ Pclass3 + Sex, family = binomial, data = titanic)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.1603  -0.5310   0.4516   0.8910   2.0143  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.8878     0.1776 -10.629   <2e-16 ***
Pclass3       1.5125     0.1783   8.483   <2e-16 ***
Sexmale       2.6067     0.1818  14.337   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1186.66  on 890  degrees of freedom
Residual deviance:  838.86  on 888  degrees of freedom
AIC: 844.86

Number of Fisher Scoring iterations: 4

手計算
初期値は経験ロジットによる回帰パラメータを与える（LINEST関数で推定）。

【更新１回目】

【更新２回目】

【更新３回目】

【更新４回目】

更新４回目で回帰パラメータの前後に差がなくなったため終了する。
推定パラメーターはRと同じ。

	β	exp(β)
X0	-1.88780	0.15140
X1	1.51251	4.53810
X2	2.60674	13.55479

チケットが1st/2ndクラスに比べて3rdの方が4.5倍オッズが高い。
女性に比べて男性の方が13.6倍オッズが高い。

参考

https://ocw.tsukuba.ac.jp/wp/wp-content/uploads/2019/10/5b43d44b7c8f6475622e63173d226184.pdf
newton
Estimating Logistic Regression Coefficents From Scratch (R version) - The Pleasure of Finding Things Out
Estimating Logistic Regression Coefficents From Scratch (Python version) - The Pleasure of Finding Things Out
http://homepage.stat.uiowa.edu/~luke/classes/STAT7400-2021/_book/optimization.html

2023-04-10

【統計】ウィルコクソンの符号付き順位検定（Exact）

R、R Studio SAS 統計

Wilcoxonの符号付き順位検定のExact法（正確確率検定）のメモ。

【目次】

計算例
プログラムコード
参考

正規近似とT近似での検定は以下記事。
cochineal19.hatenablog.com

計算例

簡単なデータを準備。
下図は正規近似・T近似での符号付き順位検定の記事で使った表。
ここにSinged Rank列（符号付き順位、前後差がマイナスなら順位にマイナスを付ける）とSinged Rank列の合計を追加（紫色のセル）。

観測値（前後差の絶対値）の順位（今回は１、2.5、2.5、４、５位）を固定した場合に、符号付き順位がとり得る全パターンとその順位和を計算する（ゼロ処理はWilcoxon法：前後差ゼロは含めない、タイは平均順位）。下図の2⁵=32パターンが考えられる。

上図の全パターンを順位和別に頻度集計し、生起確率を計算する。
今回、観測値の符号付順位和は７なので、順位和=７とそれより珍しいパターンの確率を足し合わせる。
順位和=７～15が該当し、この確率の和が片側確率、その２倍が両側確率となる。

Rで実行。
coinパッケージの wilcoxsign_test を使用。
zero.method="Wilcoxon" でWilcoxon法を指定。
distribution="exact" で正確確率を指定。

> ads0 <- data.frame(
+   before=c(10,13,13,13,20,21,21)
+   ,after=c(19,18,18,15,14,21,21)
+ )
> coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="exact")

    Exact Wilcoxon Signed-Rank Test

data:  y by x (pos, neg) 
     stratified by block
Z = -0.9482, p-value = 0.4375
alternative hypothesis: true mu is not equal to 0

プログラムコード

library(coin)
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="exact")

proc univariate data=ads;
    var diff;
    output out=out1 probs=probs;
    ods output TestsForLocation=out2;
run;

SASでは、n<20はExact（正確確率）、n>=20ならWilcoxon法・t近似のp値を返す。

参考

SAS Help Center
符号検定とWilcoxonの符号付き順位検定

2023-04-09

【統計】Wilcoxonの順位和検定（Exact）

SAS R、R Studio 統計

Wilcoxonの順位和検定のExact法（正確確率検定）のメモ。

【目次】

計算例
プログラムコード
参考

正規近似とT近似での検定は以下記事。
cochineal19.hatenablog.com

計算例

簡単なデータを準備。
A群が３例で順位和=７、B群が４例で順位和=21。

観測値の順位（今回は１～７位）を固定した場合に、A群が配置されるパターンとその順位和を計算する。
A群３例の配置は下図の35パターンが考えられる（A群=１、B群=０）。

上図の全パターンを順位和別に頻度集計し、生起確率を計算する。
今回、A群の観測値の順位和は７なので、順位和=７とそれより珍しいパターンの確率を足し合わせる。
順位和=６と７が該当し、この確率の和が片側確率、その２倍が両側確率となる。

なお、正規近似での検定は、上記ヒストグラムに正規分布を当てはめるイメージ。

Rで実行。
coinパッケージの wilcox_test を使用。distribution="exact" で正確確率を指定。

> ads <- data.frame(
+   Group <- factor(c("A","A","B","A","B","B","B"))
+   ,AVAL <- seq(10:16))
> wilcox_test(AVAL ~ Group, data=ads, method="exact", distribution="exact")

    Exact Wilcoxon-Mann-Whitney Test

data:  AVAL by Group (A, B)
Z = -1.7678, p-value = 0.1143
alternative hypothesis: true mu is not equal to 0

SASではexactステートメントで指定する。

data ads;
  do AVAL = 10,11,13; Group=1; output; end;
  do AVAL = 12,14,15,16; Group=2; output; end;
run;
proc npar1way data=ads wilcoxon correct=no /*連続修正なし*/;
  class Group;
  var   AVAL;
  exact wilcoxon;
run;

プログラムコード

library(coin)
coin::wilcox_test(AVAL ~ Group, data=ads, method="exact", distribution="exact")

proc npar1way data=ads wilcoxon correct=no ;
  class TRT;
  var   AVAL;
  exact wilcoxon;
run;

参考

Wilcoxon Rank Sum Exact Test | Real Statistics Using Excel

2023-04-09

【統計】ウィルコクソンの符号付き順位検定（Wilcoxon Signed Rank Test）

SAS R、R Studio

ウィルコクソンの符号付き順位検定（Wilcoxon Signed Rank Test）のメモ。
ノンパラメトリックな対応のある差の検定（２時点）。

【目次】

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

正規近似、T近似について。exactは記載なし。

帰無仮説、対立仮設

・帰無仮説 $H_{0}:$ ２時点に差がない

・対立仮設 $H_{1}:$ ２時点に差がある

計算式等

■ 計算の流れ
① 前後差（変化量）を算出する。
② 前後差の絶対値に平均順位を付ける。前後差ゼロの扱いは手法により異なる。
　 Wilcoxon法：ゼロを除いて平均順位を付ける。
　 Pratt法　　：ゼロを含めて平均順位を付ける。
③ 前後差が正の平均順位和 $T^{+}$ を用いて検定統計量 $S$ を算出し、p値を得る。

■ 期待値・分散
　< Wilcoxon法 >

$\quad E(T)=\dfrac{n(n+1)}{4}$

$\quad V(T)=\dfrac{n(n+1)(2n+1) - 0.5 \sum t_{i}(t_{i}+1)(t_{i}-1)}{24}$

　< Pratt法 >

$\quad E(T)=\dfrac{n(n+1) -t_{0}( t_{0}+1)}{4}$

$\quad V(T)=\dfrac{n(n+1)(2n+1) - t_{0}(t_{0}+1)(2t_{0}+1) - 0.5 \sum t_{i}(t_{i}+1)(t_{i}-1)}{24}$

　※ n=被験者数、 $t_{i}$ =前後差の絶対値のタイ数、 $t_{0}$ =前後差ゼロの数

■ 統計量

$\quad S=T^{+} - E(T)$

　※ $T^{+}$ =前後差が正の平均順位和（Rank Sum）

■ t値、z値

$\quad t=\dfrac{S}{\sqrt{\dfrac{n \cdot V(T)-S^{2}}{n-1}}}$

$\quad z=\dfrac{|S| - Correct}{\sqrt{V(T)}}$

　※ Correct=連続修正。有：0.5、無：0。

計算例

架空データで計算。

① A列とB列が観測値。C列が前後差。

② 前後差の絶対値に平均順位を付ける。Wilcoxon法はD～I列、Pratt法はJ～O列。
　・abs(diff)：前後差の絶対値
　　　D列=IF(C5<>0,ABS(C5),"")
　　　J列=ABS(C5)
　・RANK.AVG：平均順位
　　　E列=IF(C5<>"",RANK.AVG(D5,D:D,1),"")
　　　K列=IF(C5<>"",RANK.AVG(J5,J:J,1),"")
　・Positive RANK.AVG：前後差が正の平均順位
　　　F列=IF(C5>0,RANK.AVG(D5,D:D,1),"")
　　　L列=IF(C5>0,RANK.AVG(J5,J:J,1),"")
　・Negative RANK.AVG：前後差が負の平均順位
　　　G列=IF(C5<0,RANK.AVG(D5,D:D,1),"")
　　　M列=IF(C5<0,RANK.AVG(J5,J:J,1),"")
　・t：タイの数
　　　H列=IF(D5<>"",IF(COUNTIF($E$4:E5, E5)=1, COUNTIF(E:E,E5),""),"")
　　　N列=IF(COUNTIF($K$4:K5, K5)=1, COUNTIF(K:K,K5),"")
　・(t³-t)：タイの補正項用
　　　I列=IF(H5<>"", H5^3-H5, "")
　　　O列=IF(N5<>"", N5^3-N5, "")

③ 前後差が正の平均順位和 $T^{+}$ を用いて検定統計量 $S$ を算出し、p値を得る。
　今回は前後差ゼロを３つ含めたので各手法で結果が異なる。

Rで検算。
coinパッケージのデフォルトは、Pratt法・Z近似のp値を返す。連続修正はなし。
zero.method="Wilcoxon" でWilcoxon法。distribution="exact"でexact。

> library(coin)
> ads0 <- data.frame(
+   before=c(10,13,12,12,20,19,18,13,16,20,10,14,20,17,20,18,20,15,13,20,13,18,17,23,16)
+   ,after=c(17,18,16,15,13,13,18,20,10,19,10,20,16,13,13,18,14,13,17,11,11,17,11,12,12)
+ )
> coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="asymptotic")

    Asymptotic Wilcoxon Signed-Rank Test

data:  y by x (pos, neg) 
     stratified by block
Z = 1.4171, p-value = 0.1564
alternative hypothesis: true mu is not equal to 0

> coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Pratt", distribution="asymptotic") #--default

    Asymptotic Wilcoxon-Pratt Signed-Rank Test

data:  y by x (pos, neg) 
     stratified by block
Z = 1.4988, p-value = 0.1339
alternative hypothesis: true mu is not equal to 0

プログラムコード

library(coin)

#-- Z近似
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="asymptotic")
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Pratt", distribution="asymptotic") #--default
  
#-- Exact
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="exact")
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Pratt", distribution="exact")

proc univariate data=ads;
    var diff;
    output out=out1 probs=probs;
    ods output TestsForLocation=out2;
run;

SASでは、n<20はexact, n>=20ならWilcoxon法・t近似のp値を返す。
（ただし、JMPはPratt法のよう。）

参考

SAS Help Center
JMP Help
https://support.sas.com/resources/papers/proceedings-archive/SUGI94/Sugi-94-172%20Tian.pdf
Solved: Wilcoxon signed rank test : Test statistic S - JMP User Community
Wilcoxon Test • Simply explained - DATAtab
https://www.sciencedirect.com/topics/mathematics/wilcoxon-signed-rank-test

2023-03-21

【統計】一般化線形混合モデル (GLMM)

R、R Studio SAS 統計

一般化線形混合モデル (Generalized Linear Mixed Model, GLMM) のメモ。
マルチレベル分析とも言う。

【目次】

計算式
プログラムコード
G行列の構造
R行列の構造（反復測定の誤差構造）
参考

・線形回帰モデルでは、誤差が互いに独立に正規分布に従うことを前提としている。
・クラスター構造を持つデータでは、データの属する集団内での相関が考えられるため、誤差が互いに独立とは考えにくい。
・GLMMではクラスター毎の相関構造をモデリングできる。
・クラスター構造を持つデータは次のようなものが考えられる。A),B)はG行列のランダム効果、C)は反復効果。
　A) 組織をクラスター単位としたデータ（組織毎の性格・性質によるズレ）
　B) 評価者をクラスター単位としたデータ（評価者毎の採点の偏り。甘辛採点）
　C) 個人をクラスター単位とした反復測定データ/パネルデータ

計算式

固定効果モデルと混合効果モデルをそれぞれ記載。
β=固定効果, γ=G行列のランダム効果, G=回帰係数γの分散共分散行列, R=誤差εの分散共分散行列。

固定効果モデルではβとσを推定する。
混合効果モデルではβとσに加え、γ, G, Rを推定する。

【固定効果モデル】
　■ LM (Linear Model)

$\qquad Y=X' \beta+ \varepsilon , \quad \varepsilon \sim N(0,\sigma^{2})$

$\qquad E(Y) = X' \beta$

$\qquad V(Y) = V(\varepsilon)$

　■ GLM (Generalized Linear Model)

$\qquad E(Y) = g^{-1}(X' \beta)$

　【統計】一般化線形モデル（GLM） - こちにぃるの日記

【混合効果モデル】
　■ LMM (Linear Mixed Model)

$\qquad Y_{i} = X_{i}' \beta + Z_{i}' \gamma_{i} + \varepsilon_{i} , \quad \gamma\sim N(0,G),\quad \varepsilon\sim N(0,R_{i})$

$\qquad E(Y_{i}) = X_{i}' \beta$

$\qquad V(Y_{i}) = V ( Z_{i}' \gamma_{i} ) + V (\varepsilon_{i}) = Z_{i} \ G \ Z_{i}' + R_{i}$

　　クラスター単位（i番目の水準）の期待値・分散

$\qquad E(Y_{i}|\gamma_{i}) = X_{i}' \beta$

$\qquad E(Y_{i}|\gamma_{i}) = V (\varepsilon_{i}) = R_{i}$

　■ GLMM (Generalized Linear Mixed Model)

$\qquad E(Y_{i}|\gamma_{i}) = g^{-1}(X_{i}' \beta + Z_{i}' \gamma_{i})$

プログラムコード

*-- ランダム効果を含めたモデル例 ;
proc mixed data=ADS method=reml;
    class Factor Cluster ;
    model Y = Factor / s ddfm=kr ; /* 固定効果 */
    random Intercept Cluster / type=&type. g gcorr ; /* ランダム効果(ランダム切片+ランダム傾き) */
run;
  
*-- 反復効果を含めたモデル例 ;
proc mixed data=ADS method=reml;
    class Group Time SubjectId ;
    model Y = Group Time Group*Time Covariate / s ddfm=kr ; /* 固定効果 */
    repeated Time / subject=SubjectId type=&type. r rcorr ; /* 反復効果 */
    lsmeans Group*Time / pdiff cl ;
run;

SASのMIXEDプロシジャで混合効果モデルが作成できる。
・Modelステートメントで「固定効果」を指定する。
・Randomステートメントで「ランダム効果」を指定する。TYPE=でG行列の構造を指定する。オプションにgとgcorrを入れるとG行列とその相関行列も出力できる。
・Repeatedステートメントで「反復効果」を指定する。TYPE=でR行列の構造を指定する。ARやARMAなど前時点の影響度合いを考慮した構造も指定可能。オプションにrとrcorrを入れるとG行列とその相関行列も出力できる。
・推定方法のデフォルトは制限付き最尤法 (REstricted Maximum Likelihood, REML) で、最尤法 (Maximum Likelihood, ML) への切り替えもできる。

正規分布以外の分布はGENMODプロシジャやGLIMMIXプロシジャで扱えるようである。

#-- nlme
library(nlme)
fit1 <- nlme::lme(
  data=ads
  ,fixed=Y~Group*Time
  ,random=~1|SubjectId
  ,corr=nlme::corCompSymm(form=~1|SubjectId)
  ,method="REML"
)
summary(fit1)
anova(fit1)  #-- ANOVAテーブル
nlme::VarCorr(fit1)  #-- 共分散パラメータの推定
  
#-- lmer
library(lme4)
library(lmerTest)
fit1 <- lmer(Y~ Factor + (1 | Subject), ADS)) #-- ランダム切片
fit2 <- lmer(Y~ Factor + (1 + Factor | Subject), ADS)) #-- ランダム切片＋ランダム傾き
summary(fit1)
anova(fit1, ddf="Kenward-Roger")  #-- ANOVAテーブル
lme4::ranef(fit1)  #-- random effect
lme4::fixef(fit1)  #-- fixed effect