混合効果モデル
クロスオーバー試験は、生物学的同等性試験や薬物動態試験で利用されています。ウォッシュアウト期間の設定がとても大切で、治療効果が無くなるのを十分に待つことになります。このデザインはリハビリテーション研究でも利用されていますが、持ち越し効果や時期効果が十分に吟味されず、間違った使い方をしているものも見られます。これらの効果はP値に直接関連してきますので、十分に注意しましょう。
ここでは、折笠(2016)のなかにあるサンプルデータを使用して、同じ解析結果をRで算出してみます。FEVの値をpointとして、有意水準は5%で検定します。
クロスオーバー試験の研究デザイン
2群2期のクロスオーバーデザイン、2×2 デザイン、AB/BA デザイン と呼ばれています
ここでは $ Y_{ij} $ という標記を使います(2×2の行列と考えてください)
$\blacklozenge \quad i $ は順序群(行)を示す ➡ 1=Sequence1 (AB群)、2=Sequence2 (BA群)
$\blacklozenge \quad j $ は時期(列)を示す ➡ 1=I期、2=II期
$\blacklozenge \quad Y_{ijk} $ は実測値を示す
データの準備と要約
使用するパッケージ(パッケージのインストール)
library(lmerTest)
library(exactRankTests)
データセット crossover をdatに格納します(ファイルの読み込み)
dat <- read.csv("crossover.csv", header=T, fileEncoding = "UTF-8")
Y11, Y21, Y12, Y22のデータの要約
#IDを名義尺度に変換
dat$ID <- factor(dat$ID)
参考文献に合わせて基準を変更します
基準を変更する場合にはfactor変数にしなければなりません
Rは数字であれば若い数、アルファベットはAに近い文字が基準となります。差の検定を行うと、II-IやB-Aの解を出力します。折笠(2016)ではI-IIやA-Bの解を出力しています。factor変数に変換して、levels関数で調べたら分かります。
もちろんChat GPTでも教えてくれますよ
#factor変数に変更
dat$period <- factor(dat$period)
dat$treatment <- factor(dat$treatment)
#基準を調べてみます(最初に書いてある文字が基準です)
levels(dat$period)
levels(dat$treatment)
基準を変更します
dat$period <- relevel(dat$period, "II")
dat$treatment <- relevel(dat$treatment, "B")
#確認
levels(dat$period)
levels(dat$treatment)
#I期
#sequence1の治療A
Y11 <- dat[dat$sequence=="s1" & dat$treatment=="A", "point"]
#sequence2の治療B
Y21 <- dat[dat$sequence=="s2" & dat$treatment=="B", "point"]
#II期
#sequence1の治療B
Y12 <- dat[dat$sequence=="s1" & dat$treatment=="B", "point"]
#sequence2の治療A
Y22 <- dat[dat$sequence=="s2" & dat$treatment=="A", "point"]
#データの要約
summary(Y11)
summary(Y21)
summary(Y12)
summary(Y22)
箱ひげ図
name <- c("Y11", "Y21", "Y12", "Y22")
boxplot(Y11, Y21, Y12, Y22, xaxt="n")
axis(side=1, at=1:4, labels=name)
参考文献(折笠 2016)の図3の描き方
fig <- function(){
plot(
c(1, 1.2, 1.8, 2),
c(mean(Y11), mean(Y21), mean(Y12), mean(Y22)),
xlab = "", ylab = "",
ylim = c(0,5), xlim = c(0.5, 2.5),
xaxt="n")
name <- c("I期", "II期")
axis(side=1, at=c(1.1, 1.9), labels=name)
}
fig()
#標準誤差の関数を作成
stm <- function(x) sd(x)/sqrt(length(x))
arrows(
c(1, 1.2, 1.8, 2),
c(
mean(Y11) - stm(Y11),
mean(Y21) - stm(Y21),
mean(Y12) - stm(Y12),
mean(Y22) - stm(Y22)
),
c(1, 1.2, 1.8, 2),
c(
mean(Y11) + stm(Y11),
mean(Y21) + stm(Y21),
mean(Y12) + stm(Y12),
mean(Y22) + stm(Y22)
),
code=3, angle=90
)
コメント欄 『間違い』や『分かりにくい部分』などのご意見もお寄せください