Код IBM SPSS/SAS/R Project для расчета доверительных интервалов для исследования биоэквивалентности проведенному по адаптивному (adaptive/sequential) дизайну

Код для определения доверительных интервалов исследования биоэквивалентности для T-R контраста. С учетом структуры данных использована модель (ALL FIXED):
Y = Subject(Sequence*Stage) + Period(Stage) + Formulation + Sequence + Sequence(Stage) + Stage
Фактор: Formulation(Stage) – не включен.
Альфа = 0.05 (должна быть переопределена в зависимости от использованного алгоритма).

SPSS GLM:

1
2
3
4
5
6
7
8
9

GLM
LnVar BY Period Sequence Formulation Subject Stage
/CONTRAST(Formulation)=Simple(1)
/METHOD=SSTYPE(3)
/INTERCEPT=INCLUDE
/CRITERIA=ALPHA(.1)
/PRINT ETASQ PARAMETER
/EMMEANS=TABLES(Formulation) COMPARE ADJ(LSD)
/DESIGN= Subject(Sequence*Stage) Period(Stage) Formulation Sequence Sequence(Stage) Stage.

SPSS UNIANOVA:

1
2
3
4
5
6
7
8
9

UNIANOVA
LnVar BY Period Formulation Sequence Stage Subject
/CONTRAST(Formulation)=Simple
/METHOD=SSTYPE(3)
/INTERCEPT=INCLUDE
/CRITERIA=ALPHA(0.1)
/PRINT ETASQ PARAMETER
/EMMEANS=TABLES(Formulation) COMPARE ADJ(LSD)
/DESIGN=Subject (Sequence*Stage) Period(Stage) Formulation Sequence Sequence(Stage) Stage.

SAS:

1
2
3
4
5
6
7
8

PROC GLM data=pkdata;
CLASS Period Sequence Formulation Subject Stage;
MODEL logCMAX = Sequence Subject(Sequence*Stage) Period(Stage) Sequence(Stage) Stage Formulation;
*RANDOM SUBJECT(SEQUENCE*STAGE) /TEST;
LSMEANS Formulation / PDIFF=control("R") CL ALPHA=0.1;
ESTIMATE "Test - Reference" Formulation -1 1;
ods output LSMeanDiffCL=LSMeanDiffCL;
RUN;

R Project (см. комментарий):

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

lmobj <- lm(log(Param)~Formulation+Period%in%Stage+Sequence+Sequence%in%Stage+Stage+Subject%in%Sequence:Stage, data=pkdata)
lnCI <- confint(lmobj, c("FormulationT"), level=0.9)
lnPE <- coef(lmobj )["FormulationT"]
typeIII <- drop1(lmobj, test="F")
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)

Оценка ДИ соответствует значениям, полученным в Phoenix WinNonlin.

Комментарий по R Projet:
Если смотреть на описаие модели, то указанный код R Project должен ей соответствовать. Но print (typeIII) не дает полных результатов. Вызов car::Anova(lmobj, type = 3) отдает “there are aliased coefficients in the model”. На результаты оценки доверительного интервала некоторые особенности формулировки модели не влияют, если заменить на:

1	lmobj <- lm(log(Param)~Formulation+Period:Stage+Sequence+Sequence:Stage+Stage+Subject%in%Sequence:Stage, data=pkdata)

Или:

1	lmobj <- lm(log(Param)~Formulation+Period:Stage+Sequence+Sequence:Stage+Stage+Subject, data=pkdata)

То оценка доверительного интервала не изменится. Разумеется изменится оценка значимости в ANOVA.
Fixed VS Random: На мой взгляд позиция EMA (all fixed) более взвешена. Но опять же – на оценке ДИ это не сказывается.