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,983

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,646 posts in 4,328 threads, 1,436 registered users;
online 29 (0 registered, 29 guests [including 16 identified bots]).
Forum time: 09:58 CEST (Europe/Vienna)

A little Learning is a dang’rous Thing;
Drink deep, or taste not the Pierian Spring:
There shallow Draughts intoxicate the Brain,
And drinking largely sobers us again.    Alexander Pope

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