ここからはやや複雑になりますので、実際に解析されるときは、統計の専門家にご相談ください。参考までにメモを残しておきます。
複合対称型（Compound symmetry）
各観測時点における分散が等しく（等分散）、異なる観測時点間の共分散が等しいことを仮定したモデル
$\sigma^2=\sigma_{id}^2+\sigma_e^2$
$\rho=\dfrac{\sigma_{id}^2}{\sigma_{id}^2+\sigma_e^2}=ICC$（級内相関係数）
この仮定は、例えば反復測定における測定間隔が一定である場合や、クラスター内の観測値が同じように相関すると期待される場合などに適用されます。
Compound symmetryは、モデルの解釈を容易にするためにしばしば使用されますが、解析対象のデータの共分散構造に適合しているか検証しなければなりません。適合していない場合は、モデルの適合度や推定結果に誤りが生じる可能性があります。
$cov(Y_{ij}) =
\begin{pmatrix} \sigma_{id}^2+\sigma_e^2 & \sigma_{id}^2 & \sigma_{id}^2 \\ & \sigma_{id}^2+\sigma_e^2 & \sigma_{id}^2 \\ & & \sigma_{id}^2+\sigma_e^2 \\ \end{pmatrix}$
$= \sigma^2 \begin{pmatrix} 1 & \rho & \rho \\ & 1 & \rho \\ & & 1 \ \end{pmatrix}$
lme関数による Compound symmetry を指定したモデル。欠損値を除去したdat2を使用します。
model3 <- nlme::lme(
fim ~ period,
random = list(id = pdCompSymm(~ period)),
data = dat2
)
summary(model3)
list(id = pdIdent(~ period))
~ period はランダム切片とランダム傾きを考慮します。~ period – 1 は -1 によりランダム切片を除外し、ランダム傾きのみを持つモデルを指定します。
AR1型
$\blacksquare$ 一次自己回帰モデル
$X_t=\phi X_{t-1}+\epsilon_t$
$X_t: $時点tにおける観測値\
$\phi: $自己回帰パラメータ
$\epsilon_t: $時点tにおける誤差項
AR1型の共分散構造は、一次自己回帰モデル(AR1)に由来しており、隣接する観測値間の相関が減衰するという構造を持っています。
$Corr(X_t,X_{t−1})=\phi$
$\phi*Corr(X_t,X_{t−2})=\phi^2$
一次自己回帰モデル（AR1モデル）は、時系列データのモデリングで使用され、現在の観測値が前の時点の観測値に依存するという構造を持っています。AR1型の共分散構造は、一次自己回帰モデル(AR1)に由来しており、隣接する観測値間の相関が減衰するという構造を持っています。AR1型の共分散構造は、混合効果モデルや一般化推定方程式などで使用され、観測値間の相関をモデル化する際に、隣接する観測値間の相関が最も強く、次第に離れた時点間の相関が減少する構造になっています。
経時データの解析で使用される場合が多い共分散構造
最も隣接する時点（”bl” と “w3″）の相関が最も強く、次第に離れた時点（”bl” と “w6″）の相関が減少する共分散構造。
$cov(Y_i) = \sigma^2 \begin{pmatrix} 1 & \phi & \phi^2 \\ & 1 & \phi \\ & & 1 \\ \end{pmatrix} $
例えば、AR1型の共分散構造において隣接する観測値間の相関が0.5の場合、1時点離れた観測値間の相関は0.5^2=0.25、2時点離れた観測値間の相関は0.5^3=0.125となる。
ランダム効果の共分散構造は、対象者ごとの違いやグループ内の変動を捉えるために使用され、 random引数 で指定されます。観測誤差の相関構造は、隣接する観測誤差（モデルの予測と実際の観測値との差）の相関を考慮するために使用され、これは correlation引数 で指定されます。
観測誤差（モデルの予測と実際の観測値との差）の相関とは、例えば「blの予測値と実際の観測値の差」と「w3の予測値と実際の観測値の差」の相関ということです！
model4 <- lme(
fim ~ period,
random = ~1 | id,
correlation = corAR1(form = ~1 | id),# (form = ~1 | id) は、相関構造が個体idごとに存在することを示しています
data = dat2
)
summary(model4)Random effects:
Formula: ~1 | id
(Intercept) Residual
StdDev: 8.753775 1.268974
Correlation Structure: AR(1)
Formula: ~1 | id
Parameter estimate(s):
Phi
0
Fixed effects: fim ~ period
Value Std.Error DF t-value p-value
(Intercept) 61.58824 2.1452941 16 28.708527 0.0000
periodw6 1.70588 0.4352546 16 3.919274 0.0012Correlation Structureは、 AR1 の相関構造を示し、そのパラメータの推定値Phi = 0 ということは、時点間の相関がゼロ であることを意味します。AR(1) 構造を適用する必要がないことを示しています。
無相関型（indepedence）
切片のランダム効果を除外して、bl, 3w, 6wの各期間の分散が等分散であり、また各期間には関連性がない（無相関）ことを仮定したモデル。したがって、今回のような経時データの解析には適用できない場合が多いモデルです。
$cov(Y_{ij}) = \begin{pmatrix} \sigma^2 & 0 & 0 \\ & \sigma^2 & 0 \\ & & \sigma^2 \\ \end{pmatrix} $
この共分散構造は、異なる期間同士のランダム効果が互いに無相関であるという仮定に基づいています。（~ period-1）の -1 で、各idのランダム効果の切片を除外して、periodに対するidのランダム効果のみを設定。
# nlmeパッケージを使用して線形混合効果モデルを作成
model5 <- nlme::lme(
fim ~ period,
random = list(id = pdIdent(~ period - 1)),
data = dat2
)
# モデルの要約を表示
summary(model5)pdIdent関数は等分散の共分散構造を定義するための関数です。ランダム効果のすべての成分が同じ分散を持ち、他の成分との共分散が0であるという構造を持つため、実質的には対角行列となります。
Random effects:
Formula: ~period - 1 | id
Structure: Multiple of an Identity
periodw3 periodw6 Residual
StdDev: 8.282093 8.282093 3.105785
Fixed effects: fim ~ period
Value Std.Error DF t-value p-value
(Intercept) 61.58824 2.145295 16 28.708510 0.0000
periodw6 1.70588 3.033906 16 0.562273 0.5817分散
var2 <- VarCorr(model5)
print(var2)> print(var2)
id = pdIdent(period - 1)
Variance StdDev
periodw3 68.59307 8.282093
periodw6 68.59307 8.282093
Residual 9.64590 3.105785無構造型 (unstructured)
$cov(Y_{ij}) = \begin{pmatrix} \sigma_{bl} & \sigma_{bl, w3} & \sigma_{bl, w6} \\ & \sigma_{w3} & \sigma_{w3, w6} \\ & & \sigma_{w6} \\ \end{pmatrix} $
model6 <- nlme::lme(
fim ~ period,
random = list(id = pdSymm(~ period - 1)),
data = dat2
)
# モデルの要約を表示
summary(model6)Random effects:
Formula: ~period - 1 | id
Structure: General positive-definite
StdDev Corr
periodw3 8.862784 perdw3
periodw6 8.689775 0.995
Residual 1.099260
Fixed effects: fim ~ period
Value Std.Error DF t-value p-value
(Intercept) 61.58824 2.1660118 16 28.433934 0.0000
periodw6 1.70588 0.4352545 16 3.919276 0.0012以下、ChatGPTの答えです
共変量で調整したモデル
$\blacksquare$ 混合効果モデルは重回帰分析と同様に多くの変数を投入して変数選択していくことも可能です。lme関数を使用して線形混合効果モデルを作成します。期間、ステージ、年齢を固定効果として、個体差を複合対称構造でランダム効果として考慮します。
multi1 <- lme(
fim ~ period + stage,
random = list(id = pdCompSymm(~ period)),
data = dat2
)
summary(multi1)Random effects:
Formula: ~period | id
Structure: Compound Symmetry
StdDev Corr
(Intercept) 5.9081534326
periodw6 5.9081534326 -0.045
Residual 0.0002571355
Fixed effects: fim ~ period + stage
Value Std.Error DF t-value p-value
(Intercept) 56.68754 2.410391 16 23.517984 0.0000
periodw6 1.70588 1.432938 16 1.190479 0.2512
stageV 5.87779 3.283618 14 1.790036 0.0951
stageVI 10.54180 3.809781 14 2.767037 0.0151
Correlation:
(Intr) perdw6 stageV
periodw6 -0.027
stageV -0.734 0.000
stageVI -0.632 0.000 0.464$\blacksquare$ 交互作用の効果も検証可能
multi2 <- lme(
fim ~ period * stage,
random = ~ 1 | id,
data = dat2
)
summary(multi2)Random effects:
Formula: ~1 | id
(Intercept) Residual
StdDev: 8.326358 1.326887
Fixed effects: fim ~ period * stage
Value Std.Error DF t-value p-value
(Intercept) 56.66667 3.442113 14 16.462755 0.0000
periodw6 2.16667 0.766079 14 2.828256 0.0134
stageV 5.90476 4.690810 14 1.258794 0.2287
stageVI 10.58333 5.442459 14 1.944587 0.0722
periodw6:stageV -0.59524 1.043989 14 -0.570157 0.5776
periodw6:stageVI -0.91667 1.211277 14 -0.756777 0.4617
Correlation:
(Intr) perdw6 stageV stagVI prd6:V
periodw6 -0.111
stageV -0.734 0.082
stageVI -0.632 0.070 0.464
periodw6:stageV 0.082 -0.734 -0.111 -0.052
periodw6:stageVI 0.070 -0.632 -0.052 -0.111 0.464参考文献
丹後 俊郎(著); 医学統計学シリーズ 10, 経時的繰り返し測定デザイン, ―治療効果を評価する混合効果モデルとその周辺―
丹後, 俊郎／松井, 茂之（編）; 医学統計学ハンドブック 新版
柳川 堯 (著); 医療・臨床データチュートリアル: 医療・臨床データの解析事例集
コメント欄 『間違い』や『分かりにくい部分』などのご意見もお寄せください