R project код для расчета доверительных интервалов исследований биоэквивалентности (non-replicated)

Для расчета 90% доверительного интервала для разницы T/R в исследовании биоэквивалентности рекомендуется применение GLM (общая линейная модель):

The Estimate statement in SAS PROC GLM, or equivalent statement in other software, should be used to obtain estimates for the adjusted differences between treatment means and the standard error associated with these differences.

Необходимо учесть, что эта рекомендация справедлива только для нерепликативного дизайна.

Наиболее простой путь для выполнения расчетов – использование функции lm().

В работе Reference datasets for 2-treatment, 2-sequence, 2-period bioequivalence studies. (doi: 10.1208/s12248-014-9661-0. Epub 2014 Sep 12.) показано выполнение расчетов именно на основе lm ().

Для начала воспользуемся датасетом предложенным здесь (Название Treatment изменено на Formulation, AUCt -> AUC).

pkdata

Subject Sequence Period Formulation AUC
1 AB 1 A 364.74595
1 AB 2 B 375.426
2 BA 1 B 595.03875
2 BA 2 A 404.9456
3 BA 1 B 471.1644
3 BA 2 A 702.8314
4 AB 1 A 233.2503
4 AB 2 B 190.39375
5 BA 1 B 257.45585
5 BA 2 A 247.41105
6 AB 1 A 178.19155
6 AB 2 B 175.37495
7 BA 1 B 381.824
7 BA 2 A 246.3878
8 AB 1 A 407.99185
8 AB 2 B 360.8275
9 BA 1 B 218.47035
9 BA 2 A 315.4761
10 AB 1 A 140.1254
10 AB 2 B 91.80655
11 AB 1 A 165.365
11 AB 2 B 269.01575
12 BA 1 B 105.5625
12 BA 2 A 87.9882
13 BA 1 B 290.142
13 BA 2 A 182.77325
14 AB 1 A 122.47815
14 AB 2 B 230.48885
15 BA 1 B 143.5487
15 BA 2 A 67.9815
16 AB 1 A 274.5791
16 AB 2 B 344.4828

[свернуть]

A – Тест

B – Референс

Для импортирования необходимо:

В среде R подготовить команду (не выполнять, только напечатать): pkdata <- read.table(file="clipboard", header=TRUE, colClasses=c('factor', 'factor', 'factor', 'factor', 'numeric'))
Выделить данные вместе с заголовками и скопировать их в буфер обмена (Ctrl+C)
Выполнить команду из п. 1
Проверить правильность импорта данных выполнив: print (pkdata)

Также можно использовать: pkdata <- read.delim("clipboard", header=TRUE, sep=" ", colClasses=c('factor', 'factor', 'factor', 'factor', 'numeric')) , указав флаг импорта заголовков и сепаратор.

Для lm() важно отметить переменные содержащие факторы как “факторы” с помощью colClasses=c() или отдельно с помощью прямого указания: pkdata$Subject <- factor(pkdata$Subject), некоторые процедуры могут быть устойчивы к тому, что переменные не указаны как факторы, но на это лучше не расчитывать и указывать факторы явно.

1
2
3
4
5
6
7
8
9
10
11
12
13
14

pkdata <- read.table(file="clipboard", header=TRUE, colClasses=c('factor', 'factor', 'factor', 'factor', 'numeric'))
lmobj <- lm(-log(AUC)~Formulation+Period+Sequence+Subject, data=pkdata)
typeIII <- drop1(lmobj, test="F") #Anova Type III table
lnPE     <- coef(lmobj )["FormulationB"]
lnCI     <- confint(lmobj, c("FormulationB"), level=0.9)
mse      <- summary(lmobj )$sigma^2
CV <- 100*sqrt(exp(mse)-1)
PE <- exp(lnPE)*100
CI <- exp(lnCI)*100
print (typeIII)
print (CV)
print (PE)
colnames (CI)<- c("Lower level","Upper level")
print(CI)

Выполнив этот код вы получите оценку 90% доверительного интервала, точечную оценку, CV-intra и таблицу ANOVA Type III. Теперь подробно.

lmobj <- lm(-log(AUC)~Formulation+Period+Sequence+Subject, data=pkdata)

lmobj – объект модели lm(), состав можно увидеть вызвав: attributes(lmobj), отдельно стоит отметить выражение -log(AUC). Так как R Project адаптирован к векторным вычислением, то с векторами можно обращаться как с числами. Т.е. log(AUC) логарифмирует каждое значение вектора AUC, что избавляет нас от подготовки датасета и отдельного логарифмирования каждого значения (о минусе позже).

lnPE <- coef(lmobj )["FormulationB"] – точечное значение, которое можно найти вызвав print(lmobj). Значение "FormulationB" зависит от вашего датасета и от того, как названы в нем факторы. К примеру, в датасете указаной выше статьи название переменной-фактора Trt c двумя уровнями (T, R), следовательно, значение в кавычках будет TrtT. Т.е. это фактор, который будет сравниваться – тестируемый фактор (опорный, референсный фактор – первый (R)). Если уровни фактора названы иначе (к примеру: A, B – как в текущем примере), тогда, в случае если фактор FormulationB – референсный, то используются отрицательные (умноженные на -1) значения: -log(AUC), если тестовый, то log(AUC) – не измененные (в данном случае B референс – следовательно, вектор будет “обратным”).

lnCI <- confint(lmobj, c("FormulationB"), level=0.9) – вычисление интервалов с указанным уровнем (так как интервал 90% , то указывается 0.9).

Полный код для R Project статьи (Helmut Schütz, Detlew Labes, Anders Fuglsang) доступен здесь.

Код умышленно применен на другом датасете – по изменениям можно понять, как будет меняться код в зависимости от особенностей данных.

P.S. Не забывайте валидировать код. После того, как вы собрали и подготовили код – анализируйте им как минимум сбалансированный датасет и датасет с дисбалансом по последовательностям. Желательно проводить полуную валидаци по всем датасетам: