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

対応のあるｔ検定
混合効果モデル
多重比較

対応のあるｔ検定

3週目と6週目の差の検定を実行してみましょう（6週目－3週目）

これは対応のある2群の差の検定ですので、t検定で大丈夫です。ただし欠損しているデータが含まれているので、除外して検定しなければなりません。

id12とid16に欠損があるので、id12とid16のデータを削除して検定します

欠損のある対象者（id12, id16）とperiod=bl を全て除外したデータセット dat2 を作成

dat2 <- subset(
    dat,
    id != 12 & id != 16 & period != "bl"
)
summary(dat2)

グラフ作成

#dat2を使用してグラフを作成
g4 <- ggplot(dat2, aes(x = period, y = fim, group = id, colour = id)) +
        geom_line() +
        scale_x_discrete(
            "観測時点", 
            labels = c(
                "w3" = "3週後", 
                "w6" = "6週後"
            )
        )

print(g4)

欠損値のある対象者（id12, id16）を削除したので、全て関連性があります

$\blacksquare$ 対応のあるｔ検定

横のデータセット yoko から欠損のある対象者（ id12 、 id16）を除外したセット yoko2 を作成（n=18）

# データyokoの行（12行目と16行目）を除外して新しいデータセットを作成
yoko2 <- yoko[c(-12, -16), ]
View(yoko2)

観測時点6週後と3週後のデータ間のt検定を実施

t.test(yoko2$w6 - yoko2$w3)

p<0.05で６週後のFIMが有意に高い結果となりました

混合効果モデル

反復測定分散分析（Repeated Measures ANOVA）と混合効果モデル
ANOVAにおいて、切片や傾きを被験者ごとに異なるランダム変数としてモデル化する場合、その分析は混合効果モデル（Mixed Effects Model）の一形態とみなされます。観測時点の効果を固定効果、各対象者の効果をランダム効果（変量効果）とした混合効果モデルを適用することで、各観測時点のデータの関連性を配慮した解析が可能となります。

例題を簡単なモデルにすると以下のようになります（簡単のためstageのランダム効果は除く）

$Y_{ij}=\mu+\pi_i+\gamma_j+b_i+\varepsilon_{ij}$

$i=1,…,n_j \quad j=bl, w3, w6$

$n_{bl}=20,\, n_{w3}=19, \, n_{w6}=18$

$\mu:$ 母平均

$\pi_i: $ 対象者 $i$ が属する群の効果（この例では stage や age の固定効果）

$\gamma_j:$ 時期の固定効果 $\quad$ j=bl (ベースライン),$\,$ j=w3 (３週間後),$\,$ j=w6 (6週間後)

$b_i: $ 対象者 $i$ の変量効果（各対象者 $i$ のばらつき） $\quad \sim N(0,\, \sigma_{id}^2)$

$\varepsilon_{ij}:$ 誤差（各対象者 $i$ の各観測時点 $j$ における効果） $\quad \sim N(0,\, \sigma_e^2)$

通常の線形回帰モデルは、固定効果のパラメータ（$\mu, \pi_i, \gamma_j, \varepsilon_{ij}$）のみを推定するために使用されます。このモデルでは、対象者の効果を表すランダム効果 $b_i$（または変量効果）が正規分布を仮定した確率分布としてモデル化されます。固定効果とランダム効果の両パラメータを推定するために、混合効果モデルと称されます。

混合効果モデルによる対応のある2群の比較

欠損値をカットしたデータセット dat2 を使用

対象者のデータが反復して測定されているので、対象者をランダム効果（変量効果）としてモデルに投入します

# lmer関数を使用した線形混合効果モデルによる解析
fit1 <- lmer(
    fim ~ period + (1 | id), # (1 | id) これで対象者をランダム効果として指定
    data = dat2
)

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

固定効果の結果に着目

periodw6 の Estimate は 6w と 3w 差の平均値を示します

ｔ値、p値ともに 対応のあるt検定 と同じ値を示しています

欠損のある対象者も含めて解析

混合効果モデルを利用することで、欠損のある対象者も含めた解析が可能になります

データセット dat の period=3w, 6w を抜き出したデータセット（dat3）を使用

dat3

# datからbl（ベースライン）を除いたデータセットdat3を作成
dat3 <- subset(dat, dat$period != "bl")
summary(dat3)

混合効果モデル

# lmer関数 を使用して混合効果モデルで解析
fit2 <- lmer(
    fim ~ period + (1 | id),
    data = dat3
)

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

欠損値のあるデータを除外した場合とはやや異なる結果となります

経時データの相関

横のデータ（yoko）を使用

# yokoデータセットから"bl", "w3", "w6"のカラムを選択し、
# 上部パネルを削除して散布図行列を作成
pairs(
    yoko[, c("bl", "w3", "w6")],
    upper.panel = NULL
)

欠損しているデータを除いた yoko2 を使用

#分散・共分散行列
cov(yoko2[,c("bl","w3","w6")])

#相関行列
cor(yoko2[,c("bl","w3","w6")])

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

$\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”