クロスオーバー試験(混合効果モデル)

混合効果モデル

クロスオーバー試験は、生物学的同等性試験や薬物動態試験で利用されています。ウォッシュアウト期間の設定がとても大切で、治療効果が無くなるのを十分に待つことになります。このデザインはリハビリテーション研究でも利用されていますが、持ち越し効果時期効果が十分に吟味されず、間違った使い方をしているものも見られます。これらの効果は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} $ は実測値を示す

データの準備と要約

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

R
library(lmerTest)
library(exactRankTests)

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

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

Y11, Y21, Y12, Y22のデータの要約

R
#IDを名義尺度に変換
dat$ID <- factor(dat$ID)

参考文献に合わせて基準を変更します

基準を変更する場合にはfactor変数にしなければなりません

Rは数字であれば若い数、アルファベットはAに近い文字が基準となります。差の検定を行うと、II-IやB-Aの解を出力します。折笠(2016)ではI-IIやA-Bの解を出力しています。factor変数に変換して、levels関数で調べたら分かります。

もちろんChat GPTでも教えてくれますよ

R
#factor変数に変更
dat$period <- factor(dat$period)
dat$treatment <- factor(dat$treatment)
#基準を調べてみます(最初に書いてある文字が基準です)
levels(dat$period)
levels(dat$treatment)

基準を変更します

R
dat$period <- relevel(dat$period, "II")
dat$treatment <- relevel(dat$treatment, "B")
#確認
levels(dat$period)
levels(dat$treatment)
R
#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)

箱ひげ図

R
name <- c("Y11", "Y21", "Y12", "Y22")
boxplot(Y11, Y21, Y12, Y22, xaxt="n")
axis(side=1, at=1:4, labels=name) 

参考文献(折笠 2016)の図3の描き方

R
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
)

コメント欄 『間違い』や『分かりにくい部分』などのご意見もお寄せください

タイトルとURLをコピーしました