Another suggestion for a R-solution [General Statistics]
Dear Shuanghe, dear all,
here my two cents (only base R, code mostly comments
):
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 (
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)
).—
Regards,
Detlew
Regards,
Detlew
Complete thread:
- SAS and R (?) for variability comparison (FDA NTID Guidance) Shuanghe 2015-09-25 18:48 [General Statistics]
- 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
- 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
- OT: LaTeX, MathML, PNGs Shuanghe 2015-09-28 17:32
- 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
- Suggestion Shuanghe 2015-09-29 18: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
- Suggestion Helmut 2015-09-29 02:22
- 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
- SAS and R (?) for variability comparison (FDA NTID Guidance) jag009 2015-09-25 22:42