Bioequivalence and Bioavailability Forum

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

posted by d_labes  – Berlin, Germany, 2015-09-29 16:26 (3901 d 06:36 ago) – Posting: # 15490
Views: 33,288

Dear Shuanghe, dear all,

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

# Get the data according to Shuanghe
emadata  <- "https://dl.dropboxusercontent.com/u/7685360/permlink/emabe4p.csv"
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[!is.na(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")
options(contrasts=c("contr.sum","contr.poly"))
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)))
options(oc)

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 .

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

Complete thread:

• SAS and R (?) for variability comparison (FDA NTID Guidance) Shuanghe 2015-09-25 18:48 
  • SAS and R (?) for variability comparison (FDA NTID Guidance) jag009 2015-09-25 22:42
    • SAS and R (?) for variability comparison (FDA NTID Guidance) Shuanghe 2015-09-28 17:05
      • Missings d_labes 2015-09-30 10:52
        • Missings Shuanghe 2015-09-30 15:34
  • LaTeX, MathML, PNGs Helmut 2015-09-27 14:23
    • OT: LaTeX, MathML, PNGs Shuanghe 2015-09-28 17:32
      • OT: LaTeX, MathML, PNGs Helmut 2015-09-28 17:46
  • Quantiles of the F-distribution d_labes 2015-09-28 09:45
    • Quantiles of the F-distribution Shuanghe 2015-09-28 17:39
  • What's wrong with this picture? Shuanghe 2015-09-28 18:57
    • Suggestion Helmut 2015-09-29 02:22
      • Suggestion Shuanghe 2015-09-29 18:44
        • OT: The Hadleyverse d_labes 2015-09-30 08:44
    • Another suggestion for a R-solutiond_labes 2015-09-29 14:26
      • Wow – another gem from the master! Helmut 2015-09-29 16:31
      • Another suggestion for a R-solution Shuanghe 2015-09-29 18:37
        • Another R-solution d_labes 2015-09-29 21:55
      • @Shuanghe: Outdated URL gvk 2019-05-21 11:36
  • SAS and R (?) for variability comparison (FDA NTID Guidance) M.tareq 2017-04-14 15:04
    • SAS and R (?) for variability comparison (FDA NTID Guidance) M.tareq 2017-04-15 01:01