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

混合効果モデル、反復測定分散分析

ランダム効果と固定効果

通常の線形回帰モデル

通常の線形回帰モデルは、目的変数 $y_i$を説明変数 $x_i$ と誤差項 $e_i$ を用いて説明します。モデルは次のように表されます

$y_i = \beta_0 + \beta_1x_i + e_i$

ここで、$\beta_0$ は切片、$\beta_1$ は傾き（$x_i$ の効果）を表し、これらは固定効果です。誤差項 $e_i$ は、通常、平均 0 と一定の分散を持つ正規分布に従うと仮定されます。このモデルでは、$\beta_0$ と $\beta_1$ の固定効果をデータから推定します。

線形混合効果モデル

このモデルでは、切片と傾きにランダム効果（変量効果）を加えることで、データ内のグループ間の変動を考慮します。モデルは次のように拡張されます：

$y_i = (\beta_0+\alpha_0) + (\beta_1+\alpha_1)x_i + e_i$

$y_i = \beta_0+\beta_1x_i+(\alpha_0+\alpha_1x_i)+e_i$

ここで、$\alpha_0$ はランダム切片、$\alpha_1$ はランダム傾きを表し、これらはそれぞれのグループや個体が持つ固有の変動を捉えるためのものです。$\alpha_0+\alpha_1x_i$ はランダム効果（変量効果）であり、これらも通常、正規分布に従うと仮定されます。

線形混合効果モデルでは、固定効果 $(\beta_0、\beta_1)$ とランダム効果 $(\alpha_0, \alpha_1)$ の両方を同時に推定します。こうして、データ内の固定された傾向とランダムな変動の両方をモデル化できます。そのため、これを「混合効果」と呼びます。

データの準備と要約

例題）20名の脳卒中患者に理学療法を6週間実施して、実施前（ベースライン）、3週後、6週後の各時点でのFIM得点に差があるかどうか、有意水準は5%で検証せよ。

僕が作った架空のデータで練習しましょう。例題の内容や結果については色々とご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。この例題は反復測定の解析になります。各観測時点のデータが独立していないので〇元配置分散分析で解析するのは間違いです。

データセット glmm をdatに格納します（ファイルの読み込み）

glmm サンプルのダウンロード

dat <- read.csv("glmm.csv", header=T, fileEncoding = "UTF-8")

$\blacksquare$ stage：ブルンストロームステージ
$\blacksquare$ period: 観測時点（bl=ベースライン、w3=3週後、w6=6週後）
$\blacksquare$ fim: FIM運動項目の得点

View(dat)

使用するパッケージ（パッケージのインストール）

#混合効果モデルのパッケージ
library(lmerTest)
library(nlme)

#混合効果モデルの多重比較
library(multcomp)

#tidyrとggplot2を使用 
library(tidyverse)

#多重比較のためのパッケージ
library(multcomp)

変数の変更と基準の設定

ファクターに変更

#id, period, stageをファクターに変更
dat$id <- factor(dat$id)
dat$period <- factor(dat$period)
dat$stage <- factor(dat$stage)

基準の変更

# periodoの基準を bl に設定
dat$period <- relevel(dat$period, ref="bl") 

# stageの基準を IV に設定
dat$stage <- relevel(dat$stage, ref="IV")

データの要約（縦のデータセット）

head(dat)
summary(dat)
dim(dat)

＊w3とw6に欠損があります

横データの作成

yoko <- tidyr::spread(dat, key=period, value=fim)
View(yoko)

id12 の 3w, 6w と id16 の 6w の欠損が確認できます

グラフでデータを確認

# 基本的なグラフの設定
g1 <- ggplot2::ggplot(dat, aes(x = period, y = fim)) +
        scale_x_discrete(
            "観測時点", 
            labels = c(
                "bl" = "ベースライン", 
                "w3" = "3週後",
                "w6" = "6週後"
            )
        )

# ジッターの追加とテーマ性の設定
g2 <- g1 +
        geom_jitter(height = 0, width = 0.1, size = 2) +
        theme(
            axis.title.x = element_text(size = 12),
            axis.title.y = element_text(size = 12),
            axis.text.x = element_text(size = 12),
            axis.text.y = element_text(size = 12)
        )

# グラフの出力
print(g2)

ネストしているデータ構造（マルチレベル）

上のグラフに stage と id の情報を追加します

g3 <- g1 + 
    geom_jitter(height = 0, width = 0.1, size = 3, aes(colour = stage)) +
    geom_text(mapping = aes(label = id), vjust = 1.5) 

print(g3)

見づらい図ですが、おおよそこのようなイメージです（ageは記載してません）。これが「ネストしている」と言われるデータの構造です。図中の数字はidです（各観測時点で対応していることが理解できます）。

各stageごとに色分けして各時点を線で結びます。ただしデータを重ねて、IDは省略します。

# 基本的なグラフの設定とステージごとの色分け
g3 <- g1 + 
  geom_jitter(aes(colour = stage), height = 0, width = 0, size = 3) +
  geom_text(mapping = aes(label = ""), vjust = 1.5)

# ステージの色に合わせて線を追加
g4 <- g3 + 
  geom_line(aes(group = id, colour = stage), linewidth = 1)

print(g4)