経時データの解析（混合効果モデル、反復測定分散分析）

相関を考慮しなければならない理由

$\blacksquare$ 検定力低下

相関を考慮しない場合、標準誤差（SE）の推定が不正確になる可能性がある（過大または過小な推定）

結果として、検定力が低下する可能性がある

例えば、対応のある２群（A1, A2）の差の不偏標本分散は・・・

$V[A2-A1]=V[A2]+V[A1]-2cov(A2, A1)$

対応ある２群の相関を無視してしまうと不偏標準偏差が間違った値となることが理解できます。証明については「確率変数の差の分散」などをキーワードで検索してみてください

$\blacksquare$ 偏りの問題

ランダム効果を無視すると、固定効果の推定値が偏る可能性がある

固定効果やランダム効果が過剰な値に推定される可能性があり、非常に小さい効果でも統計的に有意とされる可能性がある

$\blacksquare$ Type I エラーの増加

複数の比較（例えば、idごと、periodごと、stageごとなど）が行われる場合、それぞれの比較でType Iエラーが生じる確率が独立に高まる可能性が生じる

$\blacksquare$ 信頼区間の不正確性

上述した問題から信頼区間が不正確になる場合がある

切片モデル（ランダム切片モデル）

対象者ごとのランダム効果の切片の変動を取り入れて、そのばらつき（分散）を定量化するモデル

# lmerTestパッケージを使用して線形混合効果モデルを作成
model1 <- lmerTest::lmer(
    fim ~ period + (1 | id), # 固定効果とランダム効果の定義
    data = dat               # 使用データセット
)

# モデルの要約を表示
summary(model1)

3週目とblには差があり、また6週目とblには差があることが理解できます

3Wと6Wの比較は基準を変更するか、多重比較で検定可能です

model1から分散分析表の算出も可能です

anova(model1)

Random effects の Variance

ランダム効果の切片の分散：$\hat{\sigma}_{id}^2=61.04$（母分散を推定した分散）

Correlation of Fixed Effects (固定効果パラメータの相関行列)

固定効果のパラメータが互いにどの程度関連しているかを示す

# 固定効果パラメータの分散共分散行列を取得
cov_mod <- vcov(model1)
# 標準偏差を計算
std_mod <- sqrt(diag(cov_mod))
# 相関行列を計算
cor_mod <- cov_mod / (std_mod %o% std_mod)
print(cor_mod)

混合効果モデルの分散や相関行列の導出、またパラメータの推定や標準偏差の導出（Kenward-Roger、Satterthwaite）については、またの機会にアップしたいと思います。以下、ChatGPTの答えです。

切片・傾きモデル（ランダム切片、ランダム傾きモデル）

注意）
傾きを評価するためにperiodを連続変数として扱います（連続変数period_nを追加）
ベースラインを0、３週目を3、６週目を6と設定します

library(dplyr)
dat <- dat %>%
  mutate(period_n = case_when(
    period == "bl" ~ 0,
    period == "w3" ~ 3,
    period == "w6" ~ 6,
    TRUE ~ NA_real_  # それ以外の場合はNAを代入
  ))

ランダム切片、ランダム傾きモデル

model2 <- lmerTest::lmer(
    fim ~ period_n + (1 + period_n| id), 
    data = dat             
)
summary(model2)

多重比較

periodに有意な効果を認めたのでmodel1の多重比較

mc <- multcomp::glht(
    model1,
    linfct = mcp(period = "Tukey")
)

summary(mc)

調整（bonferroni）

summary(mc, test = adjusted("bonferroni"))

以下の調整が可能

“single-step”, “Shaffer”, “Westfall”, “free”, “holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”

複合対称型（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 を指定したモデル

model3 <- nlme::lme(
    fim ~ period,
    random = list(id = pdCompSymm(~ period)),
    data = dat
)

summary(model3)

list(id = pdIdent(~ period)) 

(~ period) と書くと、切片も傾きもランダム効果となります。(~ period-1) と書くと、 -1 がランダム切片をモデルから除外することを意味することになり、対象者ごとのランダム傾きのみを持つモデルを指定します。

periodとstageの両方にCompound symmetry構造を想定した場合

model3_2 <- lme(
    fim ~ period, 
    random = list(id = pdCompSymm(~ period), id = pdCompSymm(~ stage)), 
    data = dat
)

summary(model3_2)