Another suggestion for a R-solution [General Sta­tis­tics]

posted by d_labes  – Berlin, Germany, 2015-09-29 14:26  – Posting: # 15490
Views: 19,491

Dear Shuanghe, dear all,

here my two cents (only base R, code mostly comments :cool:):

# Get the data according to Shuanghe
emadata  <- ""
ema <- read.csv(file = emadata, stringsAsFactors = FALSE)

# sort by subject, treatment and period to get replicate number
# this works in contrast to the SAS code and Shuanghe's R solution also for
# other sequences than TRTR|RTRT

ema <- ema[order(ema$subj, ema$treat, ema$per),]
ema$repl <- replicate(nrow(ema),1)
for (i in c(2:nrow(ema))){
  if (ema$subj[i] != ema$subj[i-1] | ema$treat[i] != ema$treat[i-1]) ema$repl[i] <- 1 else
        ema$repl[i] <- ema$repl[i-1]+1
# make a column with treatment and repl, i.e T1,T2, R1,R2
ema$code2 <- paste0(ema$treat, ema$repl)
ema <- ema[order(ema$subj, ema$per),]
#now reshape to wide with cols subj, seq, pk.R1, pk.R2, pk.T1, pk.T2
ema_r <- reshape(data=ema, direction="wide", v.names="logpk", idvar=c("subj","seq"),
                 timevar="code2", drop=c("pk", "per", "treat", "repl"))
# change logpk.T1, T2 ... to T1, T2 ... to avoid keystrokes
names(ema_r) <- sub("logpk.", "", names(ema_r) )

# now T-R analysis, ilat in the SAS code
# next will give TR = NA if any T1, T2, R1, R2 is missing
# must adapted if 3-period TRT|RTR is used (not recommended by the FDA

ema_r$TR <- 0.5 * (ema_r$T1 + ema_r$T2 - ema_r$R1 - ema_r$R2)
# now get rid of subjects not having all 4 periods which have TR = NA
ema_r <- ema_r[!$TR),] # ilat analysis, ANOVA with seq as effect
ema_r$seq <- as.factor(ema_r$seq)
# with standard contr.treatment we get not the wanted intercept!
# with that the intercept is intercept+seq1

oc <- options("contrasts")
m1 <- lm(TR ~ seq, data=ema_r)
# intercept is estimate of µt-µR
est <- coef(m1)[1]
pe  <- exp(est)
# lower, upper 90% CI
alpha <- 0.05
CI    <- confint(m1, 1, level=1-2*alpha)
dfTR  <- m1$df.residual
# stderr, "the unknown x" for unbiased estimate of (µT-µR)^2
stderr <- summary(m1)$coefficients["(Intercept)","Std. Error"]
# linearized scABE criterion component x
x      <- est^2-stderr^2
boundx <- max(abs(CI))^2
# now the equivalent of SAS code for R-R and T-T, dlat analysis
ema_r$RR <- ema_r$R1 - ema_r$R2
ema_r$TT <- ema_r$T1 - ema_r$T2
m2   <- lm(RR ~ seq, data=ema_r)
dfRR <- m2$df.residual
s2wR <- summary(m2)$sigma^2/2
# regulatory setting for NTID's, Warfarin guidance
theta <- ((log(1.11111))/0.1)^2
# progesterone guidance for HVD's
#theta <- ((log(1.25))/0.25)^2
# linearized scABE criterion component y
y         <- -theta*s2wR
boundy    <- y*dfRR/qchisq(0.95,dfRR);
swR       <- sqrt(s2wR)
# linearized scABE criterion, 95% upper bound
crit      <- x+y
critbound <- (x+y)+sqrt(((boundx-x)^2)+((boundy-y)^2))
# now TT variability
m3   <- lm(TT ~ seq, data=ema_r)
dfTT <- m3$df.residual
s2wT <- summary(m3)$sigma^2/2
swT  <- sqrt(s2wT)
# ratio of swT/swR and 90% CI
alpha2  <- 0.1
sRatio  <- swT/swR
sRatioCI <- c(swT/swR/sqrt(qf(alpha2/2,df1=dfTT, df2=dfRR, lower.tail=FALSE)),
              swT/swR/sqrt(qf(1-alpha2/2,df1=dfTT, df2=dfRR, lower.tail=FALSE)))

Encapsulating this in a function, adding the BE decision (critbound <0 and sRatioCI[2]<=2.5 in case of NTID's) and output (printing) of details for use in a statistical report is left as your homework :-D.

BTW: EMAI is not NTID with low variability but HVD (CVwR = se2CV(swR)).



Complete thread:

 Admin contact
20,471 posts in 4,300 threads, 1,415 registered users;
online 10 (2 registered, 8 guests [including 3 identified bots]).
Forum time (Europe/Vienna): 12:35 CEST

We should not speak so that it is possible
for the audience to understand us,
but so that it is impossible
for them to misunderstand us.    Quintilian

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz