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

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

C
C

もちろん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
)
クロスオーバー試験で想定している結果

対象者を各順序群に無作為に割り付けることで、両群の治療Aと治療Bのが同じくらいになることを想定しています

効果の検証

持ち越し効果(carry-over effect)がある場合

washout期間は十分に確保できていたか?

*後発医薬品の生物学的同等性試験ガイドラインでは、washout期間 (休薬期間) が設定されています

持ち越し効果の検定

各sequenceにおける治療A+治療Bの平均値の差の検定(対応のないt検定)により検証

検出力が低いので有意水準は10%が良い(折笠, 2016)

検定の結果 p<0.1 になった場合は、クロスオーバー試験は適用できません

$\blacklozenge$ サンプルの持ち越し効果の検定

$Y_{11k}+Y_{12k}$ と $Y_{21k}+Y_{22k}$ の平均の差の検定

参考文献(折笠 2016)の図4の結果

#t検定
t1 <- t.test(c(Y11 + Y12), c(Y21 + Y22) , var.equal=T)

#t検定の結果
print(t1)

#AB群平均-BA群平均
mean(Y11 + Y12) - mean(Y21 + Y22)

#AB群平均-BA群平均の標準誤差
t1$stderr

p=0.1253(>0.1)のため持ち越し効果はないと判定します

治療効果(treatment effect)の検証 <t検定と混合効果モデル>

対応のあるt検定

治療Aと治療Bの比較

参考文献(折笠 2016)の図8の結果

#対応のあるt検定(治療A-治療B)
t2 <- t.test(c(Y11, Y22) - c(Y12, Y21) ) 

#t検定の結果
print(t2)

#AB群平均-BA群平均の標準誤差
t2$stderr

p=0.042なので治療Bが有意に高いことになります

混合効果モデルでも対応のあるt検定と同じ解を得ることができます

t3 <- lmerTest::lmer(
    point ~  treatment + (1|ID),
    data = dat
)
summary(t3)

A-B=-0.2647, t値=-2.21, p値=0.0421

ウィルコクソンの符号順位検定では以下のような結果となります (パッケージexactRankTestsを使用)

exactRankTests::wilcox.exact(
    c(Y11, Y22) , c(Y12, Y21), 
    paired=T, conf.int=T
)

検定統計量の結果が参考文献(折笠 2016)の図8の結果と異なります。詳細が分かれば追加で記載します。

混合効果モデル(対象者ID=変量効果)

参考文献(折笠 2016)の図5の結果

period effectを考慮した治療Aと治療Bの差の検定

対応のある検定になるので対象者(ID)を変量効果とした混合効果モデルを適用

fit1 <- lmerTest::lmer(
    point ~ treatment + period + (1|ID),
    data = dat
)
summary(fit1)

period効果を考慮した場合の治療A-治療B平均値の推定値は-0.2565

p値=0.0472なので治療Bが有意に高いことが分かります

時期効果(period effect)の検証

参考文献(折笠 2016)の図6の結果

period と treatment の位置を変えてますが、上述した混合効果モデルと同じです

fit2 <- lmer(
    point ~ period + treatment + (1|ID),
    data = dat
)
summary(fit2)

treatmentで調整したI期—II期の平均値の推定値は0.139

p=0.2595なので有意差は認められません

反復測定分散分析:repeated measures ANOVA (混合効果モデル)

ID(被験者)を変量効果とした混合効果モデル

参考文献(折笠 2016)の図10の結果

$\blacklozenge$ 固定効果の検定

fit3 <- lmer(
    point ~ sequence + treatment + period + (1|ID), 
    data = dat
)
anova(fit3)

平方和は分散分析表から求めることができます

$\blacklozenge$ 分散分析表 (Between Subject) の平方和

Between SubjectですのでIDを説明変数に投入した線形回帰モデル

fit4 <- lm(
    point ~ sequence + ID,
    data = dat
)
anova(fit4)

$\blacklozenge$ 分散分析表 (Within-subject) の平方和

Within-subjectですのでIDを変量効果とした混合効果モデル

fit5 <- lmer(
    point ~ treatment + period + (1|ID),
    data = dat
)
anova(fit5)

事後検定

参考文献(折笠 2016)の図11の結果

分散分析を適用します

$\blacklozenge$ 効果の検定

fit6 <- lm(
    point ~ period + treatment + ID,
    data = dat
)
anova(fit6)

periodの結果が参考文献(折笠 2016)の図11の結果と異なります。原因が分かり次第修正します。

最小二乗平均

lmerTestのls_means関数を使用

ls_means(fit5)
$\blacklozenge$ I期、II期、治療A、治療Bの平均
mean(dat[dat$period=="I", "point"])
mean(dat[dat$period=="II", "point"])
mean(dat[dat$treatment=="A", "point"])
mean(dat[dat$treatment=="B", "point"])
最小二乗平均差

lmerTestのdifflsmeans関数

#参考文献に合わせます

dat$period <- relevel(dat$period, "I")
dat$treatment <- relevel(dat$treatment, "A")
#確認
levels(dat$period)
levels(dat$treatment)

fit5_2 <- lmer(
    point ~ treatment + period + (1|ID),
    data = dat
)

difflsmeans(fit5_2)

求め方・・・

$\blacklozenge$ 各群の最小二乗法推定量

繰り返しありと仮定して二元配置分散分析から求めます

fit7 <- lm(
    point ~ period + treatment,
    data=dat
)
summary(fit7)$coef

Y11(IのA):1.9611029 + 0.1390278 – 0.2565278 = 1.843603

Y12(IIのB):1.9611029

Y21(IのB):1.9611029 + 0.1390278 = 2.100131

Y22(IIのA):1.9611029 – 0.2565278 = 1.704575

I期:(1.8436029+2.100131)/2=1.971867

II期:(1.961103+1.704575)/2=1.832839

治療A:(1.8436029+1.704575)/2=1.774089

治療B:(1.961103+2.100131)/2=2.030617

$\blacklozenge$ 差の平均の最小二乗法推定量

summary(fit6)$coef[c("periodI", "treatmentA"), ]

I期-II期:0.1390278

治療A-治療B:-0.2565278

$\blacklozenge$ 差の平均の最小二乗法推定量95%信頼区間

confint ( fit6 , level = 0.95 )[c("periodI", "treatmentA"), ]

参考文献

折笠秀樹. クロスオーバー試験の計画および解析. 薬理と治療, 2016, 44.9: 1261-1276.

角間辰之, 服部聡; 臨床試験のデザインと解析: 薬剤開発のためのバイオ統計. 近代科学社. 2012

リハビリテーション研究に必要な統計学について、R(Windows, ChatGPT)を使って紹介してます。サンプルは全て架空のデータで作成しています。したがって解析結果は事実とは異なりますのでご了承ください。間違いなどのご指摘はコメント欄にご記入いただければ助かります。

統計学備忘録をフォローする

ダメ出し 間違い、分かりにくい部分などのご意見をお待ちします

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