Nothing
####################################################
# Calculate CV according to the EMA's Q&A-document #
####################################################
CV.calc <- function(alpha = 0.05, path.in, path.out, file, set = "",
ext, na, sep = ",", dec = ".", logtrans = TRUE,
regulator, ola = FALSE, details = FALSE,
adjust = FALSE, print, verbose = FALSE, ask = FALSE,
theta1 = theta1, theta2 = theta2, plot.bxp = FALSE,
fence = 2, data) {
if (missing(path.in)) path.in <- NULL
if (missing(data)) data <- NULL
called.from <- as.character(sys.call(-1))[1]
ret <- get.data(path.in = path.in, path.out = path.out, file = file,
set = set, ext = ext, na = na, sep = sep, dec = dec,
logtrans = logtrans, print = print,
plot.bxp = plot.bxp, data = data)
logtrans <- ret$logtrans
ow <- options("digits")
options(digits=12) # more digits for anova
on.exit(ow) # ensure that options are reset if an error occurs
if (logtrans) { # use raw data and log-transform internally
modCVR <- lm(log(PK) ~ sequence + subject%in%sequence + period,
data = ret$ref)
} else { # use already log-transformed data
modCVR <- lm(logPK ~ sequence + subject%in%sequence + period,
data = ret$ref)
}
aovCVR <- anova(modCVR)
msewR <- aovCVR["Residuals", "Mean Sq"]
DfCVR <- aovCVR["Residuals", "Df"]
CVwR <- mse2CV(msewR)
if (ret$design == "full") { # not for the EMA but the WHO
if (logtrans) {
modCVT <- lm(log(PK) ~ sequence + subject%in%sequence + period,
data = ret$test)
} else {
modCVT <- lm(logPK ~ sequence + subject%in%sequence + period,
data = ret$test)
}
aovCVT <- anova(modCVT)
msewT <- aovCVT["Residuals", "Mean Sq"]
DfCVT <- aovCVT["Residuals", "Df"]
CVwT <- mse2CV(msewT)
sw.ratio <- sqrt(msewT)/sqrt(msewR)
sw.ratio.CI <- c(sw.ratio/sqrt(qf(0.1/2, df1 = DfCVT, df2 = DfCVR,
lower.tail = FALSE)),
sw.ratio/sqrt(qf(1-0.1/2, df1 = DfCVT, df2 = DfCVR,
lower.tail = FALSE)))
names(sw.ratio.CI) <- c("lower", "upper")
}
outlier <- FALSE
BE.rec <- rep(NA, 2)
if (ola & !called.from == "ABE") { # check for outliers
stud.res <- rstudent(modCVR) # studentized (SAS)
stud.res <- stud.res[!is.na(stud.res)] # get rid of NAs and zeros
stud.res <- stud.res[which(stud.res != 0)]
stud.res <- stud.res[c(TRUE, FALSE)] # need only the 1st occasions
bp1 <- boxplot(stud.res, range=fence, plot=FALSE)
names.ol1 <- names(bp1$out)
stand.res <- rstandard(modCVR) # standardized (SAS, PHX/WNL)
stand.res <- stand.res[!is.na(stand.res)]
stand.res <- stand.res[which(stand.res != 0)]
stand.res <- stand.res[c(TRUE, FALSE)] # need only the 1st occasions
bp2 <- boxplot(stand.res, range=fence, plot=FALSE)
names.ol2 <- names(bp2$out)
if (length(names.ol1) > 0) { # >=1 detected
outlier <- TRUE
ol.value1 <- as.numeric(bp1$out)
ol.seq1 <- as.character(ret$ref[names(bp1$out), "sequence"])
ol.subj1 <- as.character(ret$ref[names(bp1$out), "subject"])
ol.value2 <- as.numeric(bp2$out)
ol.seq2 <- as.character(ret$ref[names(bp2$out), "sequence"])
ol.subj2 <- as.character(ret$ref[names(bp2$out), "subject"])
pars <- list(boxwex = 0.5, boxfill = "lightblue", medcol = "blue",
outpch = 21, outcex = 1.35, outcol = "red", outbg = "#FFCCCC")
overwrite <- TRUE
if (as.logical(capabilities("png"))) {
if (plot.bxp) { # save in PNG format to path.out
if (ask & file.exists(ret$png.path)) {
answer <- tolower(readline("Boxplot already exists. Overwrite the PNG [y|n]? "))
if(answer != "y") overwrite <- FALSE
}
if (overwrite) { # either the file does not exist or should be overwritten
png(ret$png.path, width = 720, height = 720, pointsize = 18)
}
}
} else {
message("png-device is not available; changed to plot.bxp = FALSE")
}
bxp(bp1, xlim = c(0, 3), ylim = c(-1, 1)*max(abs(c(ol.value1, ol.value2))),
las = 1, ylab = "residual", pars = pars, main = "")
title(main = expression(paste("EMA\u2019s model for ", CV[wR], ":")),
line = 3, cex.main = 1.1)
title(main = expression(paste("log(response) ~ sequence + subject(sequence) + period; data = R")),
line = 2, cex.main = 1.05)
title(main = bquote(paste("Outlier fence ", .(fence), "\u00D7IQR")),
line = 1, cex.main = 1.05)
if (length(names.ol1) == 0) {
lab.txt <- "no outlier"
} else {
ifelse (length(ol.value1) == 1,
lab.txt <- "1 outlier",
lab.txt <- paste(length(ol.value1), "outliers"))
} # Note: not /exactly/ equal. SAS uses 'type=2'
mtext(paste0("studentized\n(R, ~SAS)\n\n", lab.txt), 1, line = 4, at = 1)
text(rep(1.1, 2), bp1$stats[c(1, 5)], adj = c(0, 0.25), cex = 0.8,
sprintf("%+.3f", bp1$stats[c(1, 5)]))
if (!identical(ol.value1, numeric(0))) { # only if stud. outlier
text(rep(0.9, length(ol.value1)), ol.value1, adj = c(1, 0.25), cex = 0.8,
paste0("# ", ol.subj1, " (", ol.seq1, ")"))
text(rep(1.1, length(ol.value1)), ol.value1, adj = c(0, 0.25), cex =0.8,
sprintf("%+.3f", ol.value1))
}
bxp(bp2, axes = FALSE, at = 2, add = TRUE, pars = pars)
if (length(names.ol2) == 0) {
lab.txt <- "no outlier"
} else {
ifelse (length(ol.value2) == 1,
lab.txt <- "1 outlier",
lab.txt <- paste(length(ol.value2), "outliers"))
} # Note: not /exactly/ equal. SAS uses 'type=2' and PHX/WNL 'type=6'
mtext(paste0("standardized\n(R, ~SAS,\n~Phoenix WinNonlin)\n", lab.txt), 1,
line = 4, at = 2)
text(rep(2.1, 2), bp2$stats[c(1, 5)], adj = c(0, 0.25), cex = 0.8,
sprintf("%+.3f", bp2$stats[c(1, 5)]))
if (!identical(ol.value2, numeric(0))) { # only if stand. outlier
text(rep(1.9, length(ol.value2)), ol.value2, adj = c(1, 0.25), cex = 0.8,
paste0("# ", ol.subj2, " (", ol.seq2, ")"))
text(rep(2.1, length(ol.value2)), ol.value2, adj = c(0, 0.25), cex = 0.8,
sprintf("%+.3f", ol.value2))
}
abline(h = 0, lty = "dotted")
if (plot.bxp && file.exists(ret$png.path)) {
if (!is.null(dev.list()["png"])) {
invisible(dev.off(dev.list()["png"]))
}
}
# recalculate CVwR for data without outlier(s)
ol <- ret$ref[names(bp1$out), "subject"]
excl <- ret$ref[!ret$ref$subject %in% ol, ]
if (logtrans) { # use the raw data and log-transform internally
modCVR.rec <- lm(log(PK) ~ sequence + subject%in%sequence + period,
data = excl)
} else {
modCVR.rec <- lm(logPK ~ sequence + subject%in%sequence + period,
data = excl)
}
aovCVR.rec <- anova(modCVR.rec)
msewR.rec <- aovCVR.rec["Residuals", "Mean Sq"]
DfCVR.rec <- aovCVR.rec["Residuals", "Df"]
CVwR.rec <- mse2CV(msewR.rec)
if (ret$design == "full") {
sw.ratio.rec <- sqrt(msewT)/sqrt(msewR.rec)
sw.ratio.rec.CI <- c(sw.ratio.rec/sqrt(qf(0.1/2, df1 = DfCVT,
df2 = DfCVR.rec,
lower.tail = FALSE)),
sw.ratio.rec/sqrt(qf(1-0.1/2, df1 = DfCVT,
df2 = DfCVR.rec,
lower.tail = FALSE)))
names(sw.ratio.rec.CI) <- c("lower", "upper")
}
if (verbose) {
stud.res.whiskers <- signif(range(bp1$stats[, 1]), 7)
stud.res.outliers <- data.frame(ol.subj1, ol.seq1, signif(ol.value1, 7))
names(stud.res.outliers) <- c("subject", "sequence", "stud.res")
stand.res.whiskers <- signif(range(bp2$stats[, 1]), 7)
stand.res.outliers <- data.frame(ol.subj2, ol.seq2, signif(ol.value2, 7))
names(stand.res.outliers) <- c("subject", "sequence", "stand.res")
cat(paste0("\nOutlier analysis\n (externally) studentized residuals",
"\n Limits (", fence, "\u00D7IQR whiskers): ",
stud.res.whiskers[1], ", ", stud.res.whiskers[2],
"\n Outliers:\n")); print(stud.res.outliers, row.names = FALSE)
cat(paste0("\n standarized (internally studentized) residuals\n Limits (",
fence, "\u00D7IQR whiskers): ", stand.res.whiskers[1], ", ",
stand.res.whiskers[2], "\n Outliers:\n"))
# since standardized residuals are more liberal,
# we have to deal with such a special case
if (nrow(stand.res.outliers) == 0) {
cat(" none detected\n")
} else {
print(stand.res.outliers, row.names = FALSE)
}
} # EO verbose
} # EO >= 1 outlier
} # EO outlier analysis (only if called from method.A()/method.B() & ola=TRUE)
if (!called.from == "ABE") {
if (regulator == "EMA") {
reg_set <- reg_const("EMA")
BE <- as.numeric(scABEL(CV = CVwR, regulator = "EMA"))
}
if (regulator == "GCC") {
reg_set <- reg_const("GCC")
BE <- as.numeric(scABEL(CV = CVwR, regulator = "GCC"))
}
if (regulator == "HC") {
reg_set <- reg_const("HC")
if (!alpha == 0.5) {
BE <- as.numeric(scABEL(CV = CVwR, regulator = "HC"))
} else { # Cmax
BE <- c(0.80, 1.25)
}
}
} else {
BE <- c(theta1, theta2)
}
txt <- ret$txt
if (called.from != "ABE") { # only for scaling
if (regulator == "EMA") {
txt <- paste0(ret$txt,
"\nRegulator : EMA",
"\nSwitching CV : ",
sprintf("%6.2f%%", 100*reg_set$CVswitch),
"\nScaling cap : ",
sprintf("%6.2f%%", 100*reg_set$CVcap),
"\nRegul. constant (k): ",
sprintf("%.3f", reg_set$r_const))
}
if (regulator == "GCC") {
txt <- paste0(ret$txt,
"\nRegulator : GCC",
"\nSwitching CV : ",
sprintf("%6.2f%%", 100*reg_set$CVswitch))
}
if (regulator == "HC") {
if (!alpha == 0.5) {
txt <- paste0(ret$txt,
"\nRegulator : HC",
"\nSwitching CV : ",
sprintf("%6.2f%%", 100*reg_set$CVswitch),
"\nScaling cap : ",
sprintf("%6.2f%%", 100*reg_set$CVcap),
"\nRegul. constant (k): ",
sprintf("%.3f", reg_set$r_const))
} else { # Cmax
txt <- paste0(ret$txt,
"\nRegulator : HC")
}
}
if (CVwR > 0.3 & !regulator == "HC" & !alpha == 0.5) {
txt <- paste0(txt, "\nGMR restriction : 80.00% ... 125.00%")
}
}
if (ret$design == "full") {
# sw.ratio <- CV2se(CVwT)/CV2se(CVwR) # we should already have it, right?
txt <- paste0(txt, "\nCVwT : ",
sprintf("%6.2f%%", 100*CVwT))
if (called.from != "ABE") { # not needed for ABE
txt <- paste0(txt, "\nswT : ",
sprintf("%.5f", CV2se(CVwT)))
}
}
txt <- paste0(txt, "\nCVwR : ", sprintf("%6.2f%%", 100*CVwR))
if (called.from != "ABE") { # only for scaling
txt <- paste0(txt, " (reference-scaling ")
if (CVwR <= 0.3) txt <- paste0(txt, "not ")
txt <- paste0(txt, "applicable)")
if (CVwR <= 0.3) {
txt <- paste0(txt, "\nBE-limits : 80.00% ... 125.00%")
} else {
txt <- paste0(txt, "\nswR : ",
sprintf("%.5f", CV2se(CVwR)))
if (regulator == "EMA") {
txt <- paste0(txt, "\nExpanded limits : ",
sprintf("%6.2f%% ... %.2f%%",
100*BE[1], 100*BE[2]), " [100exp(\u00B1",
sprintf("%.3f", reg_set$r_const), "\u00B7swR)]")
}
if (regulator == "GCC") {
txt <- paste0(txt, "\nWidened limits :",
sprintf("%6.2f%% ... %.2f%%",
100*BE[1], 100*BE[2]))
}
}
if (ret$design == "full") {
txt <- paste0(txt, "\nswT / swR : ",
sprintf("%.4f", sw.ratio))
if (sw.ratio >= 2/3 & sw.ratio <= 3/2) { # like in PBE/IBE
txt <- paste0(txt, " (similar variabilities of T and R)")
} else {
ifelse (sw.ratio < 2/3,
txt <- paste0(txt, " (T lower variability than R)"),
txt <- paste0(txt, " (T higher variability than R)"))
}
txt <- paste0(txt, "\nsw-ratio (upper CL): ",
sprintf("%.4f", sw.ratio.CI[["upper"]]))
if (sw.ratio.CI[["upper"]] <= 2.5) { # like in the FDA's warfarin guidance
txt <- paste0(txt, " (comparable variabilities of T and R)")
} else {
txt <- paste0(txt, " (T higher variability than R)")
}
}
if (ola) {
if (outlier) {
if (ret$design == "full") sw.ratio.rec <- sqrt(msewT)/sqrt(msewR.rec)
BE.rec <- as.numeric(scABEL(CV=CVwR.rec, regulator="EMA"))
txt <- paste0(txt, "\n\nOutlier fence : ", fence,
"\u00D7IQR of studentized residuals.")
txt1 <- paste0("\nRecalculation due to presence of ",
length(ol))
ifelse (length(ol) == 1,
txt1 <- paste0(txt1, " outlier (subj. "),
txt1 <- paste0(txt1, " outliers (subj. "))
txt1 <- paste0(txt1, paste0(ol, collapse="|"), ")")
txt <- paste0(txt, txt1, "\n",
paste0(rep("\u2500", nchar(txt1)-1), collapse=""))
txt <- paste0(txt,
"\nCVwR (outl. excl.) : ", sprintf("%6.2f%%", 100*CVwR.rec),
" (reference-scaling ")
if (CVwR.rec > 0.3) {
txt <- paste0(txt, "applicable)")
txt <- paste0(txt, "\nswR (recalculated) : ",
sprintf("%.5f", CV2se(CVwR.rec)))
txt <- paste0(txt, "\nExpanded limits : ",
sprintf("%6.2f%% ... %.2f%%",
100*BE.rec[1], 100*BE.rec[2]), " [100exp(\u00B1",
sprintf("%.3f", reg_set$r_const), "\u00B7swR)]")
} else {
txt <- paste0(txt, "not applicable)")
txt <- paste0(txt, "\nUnscaled BE-limits : 80.00% ... 125.00%")
}
if (ret$design == "full") {
txt <- paste0(txt, "\nswT / swR (recalc.): ",
sprintf("%.4f", sw.ratio.rec))
if (sw.ratio.rec >= 2/3 & sw.ratio.rec <= 3/2) { # like in PBE/IBE
txt <- paste0(txt, " (similar variabilities of T and R)")
} else {
ifelse (sw.ratio.rec < 2/3,
txt <- paste0(txt, " (T lower variability than R)"),
txt <- paste0(txt, " (T higher variability than R)"))
}
txt <- paste0(txt, "\nsw-ratio (upper CL): ",
sprintf("%.4f", sw.ratio.rec.CI[["upper"]]))
if (sw.ratio.rec.CI[["upper"]] <= 2.5) { # like in the FDA's warfarin guidance
txt <- paste0(txt, " (comparable variabilities of T and R)")
} else {
txt <- paste0(txt, " (T higher variability than R)")
}
}
} else {
txt <- paste0(txt, "\n\nOutlier fence : ", fence,
"\u00D7IQR of studentized residuals.")
txt <- paste0(txt, "\nNo outlier detected.",
"\n", paste0(rep("\u2500", 49), collapse=""))
}
} # EO ola
} else { # called from ABE()
txt <- paste0(txt, "\nBE-limits : ",
sprintf("%6.2f%% ... %.2f%%", 100*BE[1], 100*BE[2]))
}
ret$txt <- txt
ret <- c(ret, CVswitch=ifelse(called.from != "ABE", reg_set$CVswitch, NA),
CVcap=ifelse(called.from != "ABE", reg_set$CVcap, NA),
r_const=ifelse(called.from != "ABE", reg_set$r_const, NA), BE=BE,
CVwT=ifelse(ret$design == "full", CVwT, NA), CVwR=CVwR,
swT=ifelse(ret$design == "full", sqrt(msewT), NA), swR=sqrt(msewR),
sw.ratio=ifelse(called.from != "ABE" & ret$design == "full",
sw.ratio, NA),
sw.ratio.upper=ifelse(called.from != "ABE" & ret$design == "full",
sw.ratio.CI[["upper"]], NA),
ol=ifelse(called.from != "ABE" & outlier, list(ol.subj1), NA),
CVwR.rec=ifelse(called.from != "ABE" & outlier, CVwR.rec, NA),
swR.rec=ifelse(called.from != "ABE" & outlier, sqrt(msewR.rec), NA),
sw.ratio.rec=ifelse(called.from != "ABE" & outlier &
ret$design == "full", sw.ratio.rec, NA),
sw.ratio.rec.upper=ifelse(called.from != "ABE" & outlier &
ret$design == "full",
sw.ratio.rec.CI[["upper"]], NA),
BE.rec=BE.rec)
return(ret)
} # end of function CV.calc()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.