R/method.B.R
In Helmut01/replicateBE: Average Bioequivalence with Expanding Limits (ABEL)
Defines functions
method.B
Documented in
method.B
#################################################
# EMA's 'Method B' (linear mixed effects) #
# fixed: sequence, period, treatment #
# random: subjects #
# option=1: nlme/lme (Satterthwaite's DF) #
# option=2: lmerTest/lmer (CONTAIN/Residual DF) #
#################################################
method.B <- function(alpha = 0.05, path.in, path.out = tempdir(), file,
set = "", ext, na = ".", sep = ",", dec = ".",
logtrans = TRUE, regulator = "EMA", ola = FALSE,
print = TRUE, details = FALSE, verbose = FALSE,
ask = FALSE, plot.bxp = FALSE, fence = 2,
data = NULL, option = 2) {
exec <- strftime(Sys.time(), usetz=TRUE)
if (!missing(ext)) ext <- tolower(ext) # case-insensitive
ret <- CV.calc(alpha=alpha, path.in=path.in, path.out=path.out,
file=file, set=set, ext=ext, na=na, sep=sep,
dec=dec, logtrans=logtrans, regulator=regulator, ola=ola,
print=print, verbose=verbose, ask=ask,
plot.bxp=plot.bxp, fence=fence, data=data)
logtrans <- ret$logtrans
# Add description of the degrees of freedom to the result file
if (option == 1) DF.suff <- "Satt"
if (option == 2) DF.suff <- "GL"
if (option == 3) DF.suff <- "KR"
results <- paste0(ret$res.file, "_MethodB_DF_", DF.suff, ".txt")
# generate variables based on the attribute
# 2nd condition: Otherwise, the header from a CSV file will be overwritten
if (!is.null(data) & missing(ext)) {
info <- info.data(data)
file <- info$file
set <- info$set
ref <- info$ref
descr <- info$descr
ext <- ""
}
os <- Sys.info()[[1]] # get OS for line-endings in output (Win: CRLF)
od <- options("digits") # save options
options(digits = 12) # increase digits for anova()
on.exit(od) # ensure that options are reset if an error occurs
if (option == 2) { # by nlme (Residual DF)
oc <- options("contrasts") # save options
options(contrasts=c("contr.treatment", "contr.poly")) # assure defaults
on.exit(oc) # ensure that options are reset if an error occurs
if (logtrans) { # use the raw data and log-transform internally
modB <- lme(log(PK) ~ sequence + period + treatment,
random = ~1|subject, na.action=na.omit,
data = ret$data)
} else { # use the already log-transformed data
modB <- lme(logPK ~ sequence + period + treatment,
random = ~1|subject, na.action=na.omit,
data = ret$data)
}
EMA.B <- summary(modB)
PE <- EMA.B$tTable["treatmentT", "Value"]
CI <- exp(as.numeric(intervals(modB, which="fixed",
level=1-2*alpha)[[1]]["treatmentT", c(1, 3)]))
DF <- EMA.B$tTable["treatmentT", "DF"]
if (verbose) {
name <- paste0(file, set)
cat("\nData set", paste0(name, ": Method B (option = 2) by lme()"),
paste0("\n", paste0(rep("\u2500", 41+nchar(name)), collapse="")), "\n")
if (logtrans) cat("Response: log(PK)\n") else cat("\nResponse: logPK\n")
print(anova(modB), digits = 6, signif.stars = FALSE)
cat("\ntreatment T \u2013 R:\n")
print(signif(EMA.B$tTable["treatmentT", c(1:2, 4:5)], 5))
cat(DF, "Degrees of Freedom (equivalent to SAS\u2019 DDFM=CONTAIN)\n\n")
}
} else { # by lmer/lmerTest (Satterthwaite's or Kenward-Roger DF)
if (logtrans) { # use the raw data and log-transform internally
modB <- lmer(log(PK) ~ sequence + period + treatment + (1|subject),
data = ret$data)
} else { # use the already log-transformed data
modB <- lmer(logPK ~ sequence + period + treatment + (1|subject),
data = ret$data)
}
if (option == 1) {
EMA.B <- summary(modB, ddf="Satterthwaite")
} else {
EMA.B <- summary(modB, ddf="Kenward-Roger")
}
PE <- EMA.B$coefficients["treatmentT", "Estimate"]
CI <- exp(PE + c(-1, +1) *
qt(1-alpha, EMA.B$coef["treatmentT", "df"]) *
EMA.B$coef["treatmentT", "Std. Error"])
DF <- EMA.B$coefficients["treatmentT", "df"]
if (verbose) {
if (option == 1) {
cat("\nData set", paste0(file, set, ": Method B (option = 1) by lmer()"),
paste0("\n", paste0(rep("\u2500", 46), collapse="")), "\n")
if (logtrans) cat("Response: log(PK)\n") else cat("\nResponse: logPK\n")
print(anova(modB, ddf="Satterthwaite"), digits=6, signif.stars=FALSE)
cat("\ntreatment T \u2013 R:\n")
print(signif(EMA.B$coefficients["treatmentT", c(1:2, 4:5)], 5))
cat(signif(DF, 6), "Degrees of Freedom (equivalent to SAS\u2019 DDFM=SATTERTHWAITE)\n\n")
} else {
df.txt <- "3; equivalent to Stata\u2019s dfm=Kenward Roger EIM)"
cat("\nData set", paste0(file, set, ": Method B (option = 3) by lmer()"),
paste0("\n", paste0(rep("\u2500", 46), collapse="")), "\n")
if (logtrans) cat("Response: log(PK)\n") else cat("\nResponse: logPK\n")
print(anova(modB, ddf="Kenward-Roger"), digits=6, signif.stars=FALSE)
cat("\ntreatment T \u2013 R:\n")
print(signif(EMA.B$coefficients["treatmentT", c(1:2, 4:5)], 5))
cat(signif(DF, 6), "Degrees of Freedom (equivalent to Stata\u2019s dfm=Kenward Roger EIM)\n\n")
}
}
} # end of evaluation by option=2 (lme) or option=1/3 (lmer)
PE <- exp(PE)
res <- data.frame(ret$type, "A", ret$n, ret$nTT, ret$nRR,
paste0(ret$Sub.Seq, collapse="|"),
paste0(ret$Miss.seq, collapse="|"),
paste0(ret$Miss.per, collapse="|"), alpha,
DF, ret$CVwT, ret$CVwR, ret$swT, ret$swR, ret$sw.ratio,
ret$sw.ratio.upper, ret$BE1, ret$BE2, CI[1], CI[2],
PE, "fail", "fail", "fail", log(CI[2])-log(PE),
paste0(ret$ol, collapse="|"), ret$CVwR.rec,
ret$swR.rec, ret$sw.ratio.rec,
ret$sw.ratio.rec.upper, ret$BE.rec1,
ret$BE.rec2, "fail", "fail", "fail",
stringsAsFactors=FALSE)
names(res)<- c("Design", "Method", "n", "nTT", "nRR", "Sub/seq",
"Miss/seq", "Miss/per", "alpha", "DF", "CVwT(%)",
"CVwR(%)", "swT", "swR", "sw.ratio", "sw.ratio.CL",
"L(%)", "U(%)", "CL.lo(%)", "CL.hi(%)", "PE(%)",
"CI", "GMR", "BE", "log.half-width", "outlier",
"CVwR.rec(%)", "swR.rec", "sw.ratio.rec",
"sw.ratio.rec.CL", "L.rec(%)", "U.rec(%)",
"CI.rec", "GMR.rec", "BE.rec")
if (ret$BE2 == 1.25) { # change column names if not scaling
colnames(res)[which(names(res) == "L(%)")] <- "BE.lo(%)"
colnames(res)[which(names(res) == "U(%)")] <- "BE.hi(%)"
}
if (regulator == "HC" & alpha == 0.5) {
if (option == 2)
stop("Please use option = 1 or option = 3 instead.")
res[which(names(res) == "L(%)")] <- 0.80
res[which(names(res) == "U(%)")] <- 1.25
colnames(res)[which(names(res) == "L(%)")] <- "BE.lo(%)"
colnames(res)[which(names(res) == "U(%)")] <- "BE.hi(%)"
}
if (!is.na(ret$BE.rec2) & ret$BE.rec2 == 1.25) { # change column names if not scaling
colnames(res)[which(names(res) == "L.rec(%)")] <- "BE.rec.lo(%)"
colnames(res)[which(names(res) == "U.rec(%)")] <- "BE.rec.hi(%)"
}
# Convert CVs, limits, PE, and CI (up to here as fractions) to percent
res$"CVwT(%)" <- 100*res$"CVwT(%)"
res$"CVwR(%)" <- 100*res$"CVwR(%)"
if ("BE.lo(%)" %in% names(res)) { # conventional limits
res$"BE.lo(%)" <- 100*res$"BE.lo(%)"
res$"BE.hi(%)" <- 100*res$"BE.hi(%)"
} else { # expanded limits
res$"L(%)" <- 100*res$"L(%)"
res$"U(%)" <- 100*res$"U(%)"
}
if (!is.na(res$"CVwR.rec(%)")) { # only if recalculated CVwR
res$"CVwR.rec(%)" <- 100*res$"CVwR.rec(%)"
if ("BE.rec.lo(%)" %in% names(res)) { # conventional limits
res$"BE.rec.lo(%)" <- 100*res$"BE.rec.lo(%)"
res$"BE.rec.hi(%)" <- 100*res$"BE.rec.hi(%)"
} else { # expanded limits
res$"L.rec(%)" <- 100*res$"L.rec(%)"
res$"U.rec(%)" <- 100*res$"U.rec(%)"
}
}
res$"PE(%)" <- 100*res$"PE(%)"
res$"CL.lo(%)" <- 100*res$"CL.lo(%)"
res$"CL.hi(%)" <- 100*res$"CL.hi(%)"
if (round(res$"CL.lo(%)", 2) >= 100*ret$BE1 &
round(res$"CL.hi(%)", 2) <= 100*ret$BE2)
res$CI <- "pass" # CI within acceptance range
if (!regulator == "HC") {
if (round(res$"PE(%)", 2) >= 80 & round(res[["PE(%)"]], 2) <= 125)
res$GMR <- "pass" # PE within 80.00-125.00%
} else {
if (round(res$"PE(%)", 1) >= 80 & round(res[["PE(%)"]], 1) <= 125)
res$GMR <- "pass" # PE within 80.00-125.00%
}
if (res$CI == "pass" & res$GMR == "pass")
res$BE <- "pass" # if passing both, conclude BE
if (!is.na(res$"CVwR.rec(%)")) {
if (round(res$"CL.lo(%)", 2) >= 100*ret$BE.rec1 &
round(res$"CL.hi(%)", 2) <= 100*ret$BE.rec2)
res$CI.rec <- "pass" # CI within acceptance range
res$GMR.rec <- res$GMR
if (res$CI.rec == "pass" & res$GMR.rec == "pass")
res$BE.rec <- "pass" # if passing both, conclude BE
}
if (details) { # results in full precision
ret <- res
if (as.character(res$outlier) == "NA") {
# remove superfluous columns if ola=FALSE or ola=TRUE
# and no outlier(s) detected
ret <- ret[, !names(ret) %in% c("outlier", "CVwR.rec(%)",
"swR.rec", "sw.ratio.rec",
"L.rec(%)", "U.rec(%)",
"CI.rec", "GMR.rec", "BE.rec",
"sw.ratio.rec.CL")]
}
if (regulator == "HC" & alpha == 0.5) {
# remove columns for Health Canada (PE of Cmax within 80.0-125.0)
ret <- ret[, !names(ret) %in% c("CL.lo(%)", "CL.hi(%)", "CI", "BE",
"log.half-width")]
}
#class(ret) <- "repBE"
return(ret)
}
# Round percents to two decimals according to the GL
res$"CVwT(%)" <- round(res$"CVwT(%)", 2)
res$"CVwR(%)" <- round(res$"CVwR(%)", 2)
if ("BE.lo(%)" %in% names(res)) { # conventional limits
res$"BE.lo(%)" <- round(res$"BE.lo(%)", 2)
res$"BE.hi(%)" <- round(res$"BE.hi(%)", 2)
} else { # expanded limits
res$"L(%)" <- round(res$"L(%)", 2)
res$"U(%)" <- round(res$"U(%)", 2)
}
res$"PE(%)" <- round(res$"PE(%)", 2)
res$"CL.lo(%)" <- round(res$"CL.lo(%)", 2)
res$"CL.hi(%)" <- round(res$"CL.hi(%)", 2)
if (!is.na(res$"CVwR.rec(%)")) { # only if recalculated CVwR
res$"CVwR.rec(%)" <- round(res$"CVwR.rec(%)", 2)
if ("BE.rec.lo(%)" %in% names(res)) { # conventional limits
res$"BE.rec.lo(%)" <- round(res$"BE.rec.lo(%)", 2)
res$"BE.rec.hi(%)" <- round(res$"BE.rec.hi(%)", 2)
} else { # expanded limits
res$"L.rec(%)" <- round(res$"L.rec(%)", 2)
res$"U.rec(%)" <- round(res$"U.rec(%)", 2)
}
}
overwrite <- TRUE # default
if (print) { # to file in UTF-8
if (ask & file.exists(results)) {
answer <- tolower(readline("Results already exists. Overwrite the file [y|n]? "))
if(answer != "y") overwrite <- FALSE
}
if (overwrite) { # either the file does not exist or should be overwritten
# only binary mode supports UTF-8 and different line endings
res.file <- file(description=results, open="wb")
res.str <- info.env(fun="method.B", option=option, path.in, path.out,
file, set, ext, exec, data)
if (os == "Windows") res.str <- gsub("\n", "\r\n", res.str) # CRLF
if (os == "Darwin") res.str <- gsub("\n", "\r", res.str) # CR
writeBin(charToRaw(res.str), res.file)
close(res.file)
}
}
# insert DF of treatment difference and alpha into the text generated by CV.calc()
cut.pos <- unlist(gregexpr(pattern="Regulator", ret$txt))
left.str <- substr(ret$txt, 1, cut.pos-1)
right.str <- substr(ret$txt, cut.pos, nchar(ret$txt))
insert.str <- "Degrees of freedom : "
if (option == 1) insert.str <- paste0(insert.str, sprintf("%7.3f (Satterthwaite)", DF))
if (option == 2) insert.str <- paste0(insert.str, sprintf("%3i", DF))
if (option == 3) insert.str <- paste0(insert.str, sprintf("%7.3f (Kenward-Roger)", DF))
insert.str <- paste0(insert.str, "\nalpha : ", alpha)
if (!alpha == 0.5) { # Health Canada Cmax
insert.str <- paste0(insert.str," (", 100*(1-2*alpha), "% CI)\n")
} else {
insert.str <- paste0(insert.str," (only PE assessed)\n")
}
txt <- paste0(left.str, insert.str, right.str)
if (!is.na(res$"CVwR.rec(%)")) {
txt1 <- paste0("\n\nAssessment based on original CVwR",
sprintf(" %.2f%%", res$"CVwR(%)"))
txt <- paste0(txt, txt1, "\n", paste0(rep("\u2500", nchar(txt1)-2), collapse=""))
}
if (!alpha == 0.5) {
txt <- paste0(txt,
"\nConfidence interval: ", sprintf("%6.2f%% ... %6.2f%%",
res$"CL.lo(%)", res$"CL.hi(%)"),
" ", res$CI,
"\nPoint estimate : ", sprintf("%6.2f%%", res$"PE(%)"),
" ", res$GMR,
"\nMixed (CI & PE) : ", res$BE, "\n")
} else {
txt <- paste0(txt,
"\nPoint estimate : ", sprintf("%6.2f%%", res$"PE(%)"),
" ", res$GMR, "\n")
}
txt <- paste0(txt,
repBE.draw.line(called.from="ABEL", L=ret$BE1, U=ret$BE2,
lo=CI[1], hi=CI[2], PE=PE), "\n")
if (!is.na(res$"CVwR.rec(%)")) {
txt1 <- paste0("\nAssessment based on recalculated CVwR",
sprintf(" %.2f%%", res$"CVwR.rec(%)"))
txt <- paste0(txt, txt1, "\n", paste0(rep("\u2500", nchar(txt1)-1), collapse=""))
txt <- paste0(txt, "\nConfidence interval: ", res$CI.rec,
"\nPoint estimate : ", res$GMR.rec,
"\nMixed (CI & PE) : ", res$BE.rec, "\n")
txt <- paste0(txt,
repBE.draw.line(called.from="ABEL", L=ret$BE.rec1, U=ret$BE.rec2,
lo=CI[1], hi=CI[2], PE=PE), "\n")
}
if (res$Design == "TRR|RTR")
txt <- paste0(txt, "Note: The extra-reference design assumes lacking period effects. ",
"The treatment\ncomparison will be biased in the presence of a ",
"true period effect.\n")
if (print & overwrite) {
res.file <- file(results, open="ab") # line endings
res.str <- txt # LF (UNIXes, Solaris)
if (os == "Windows") res.str <- gsub("\n", "\r\n", res.str) # CRLF (Windows)
if (os == "Darwin") res.str <- gsub("\n", "\r", res.str) # CR (OSX)
writeBin(charToRaw(res.str), res.file)
close(res.file)
}
} # end of function method.B()
Helmut01/replicateBE documentation built on May 9, 2022, 12:15 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions method.B
Documented in method.B
#################################################
# EMA's 'Method B' (linear mixed effects) #
# fixed: sequence, period, treatment #
# random: subjects #
# option=1: nlme/lme (Satterthwaite's DF) #
# option=2: lmerTest/lmer (CONTAIN/Residual DF) #
#################################################
method.B <- function(alpha = 0.05, path.in, path.out = tempdir(), file,
set = "", ext, na = ".", sep = ",", dec = ".",
logtrans = TRUE, regulator = "EMA", ola = FALSE,
print = TRUE, details = FALSE, verbose = FALSE,
ask = FALSE, plot.bxp = FALSE, fence = 2,
data = NULL, option = 2) {
exec <- strftime(Sys.time(), usetz=TRUE)
if (!missing(ext)) ext <- tolower(ext) # case-insensitive
ret <- CV.calc(alpha=alpha, path.in=path.in, path.out=path.out,
file=file, set=set, ext=ext, na=na, sep=sep,
dec=dec, logtrans=logtrans, regulator=regulator, ola=ola,
print=print, verbose=verbose, ask=ask,
plot.bxp=plot.bxp, fence=fence, data=data)
logtrans <- ret$logtrans
# Add description of the degrees of freedom to the result file
if (option == 1) DF.suff <- "Satt"
if (option == 2) DF.suff <- "GL"
if (option == 3) DF.suff <- "KR"
results <- paste0(ret$res.file, "_MethodB_DF_", DF.suff, ".txt")
# generate variables based on the attribute
# 2nd condition: Otherwise, the header from a CSV file will be overwritten
if (!is.null(data) & missing(ext)) {
info <- info.data(data)
file <- info$file
set <- info$set
ref <- info$ref
descr <- info$descr
ext <- ""
}
os <- Sys.info()[[1]] # get OS for line-endings in output (Win: CRLF)
od <- options("digits") # save options
options(digits = 12) # increase digits for anova()
on.exit(od) # ensure that options are reset if an error occurs
if (option == 2) { # by nlme (Residual DF)
oc <- options("contrasts") # save options
options(contrasts=c("contr.treatment", "contr.poly")) # assure defaults
on.exit(oc) # ensure that options are reset if an error occurs
if (logtrans) { # use the raw data and log-transform internally
modB <- lme(log(PK) ~ sequence + period + treatment,
random = ~1|subject, na.action=na.omit,
data = ret$data)
} else { # use the already log-transformed data
modB <- lme(logPK ~ sequence + period + treatment,
random = ~1|subject, na.action=na.omit,
data = ret$data)
}
EMA.B <- summary(modB)
PE <- EMA.B$tTable["treatmentT", "Value"]
CI <- exp(as.numeric(intervals(modB, which="fixed",
level=1-2*alpha)[[1]]["treatmentT", c(1, 3)]))
DF <- EMA.B$tTable["treatmentT", "DF"]
if (verbose) {
name <- paste0(file, set)
cat("\nData set", paste0(name, ": Method B (option = 2) by lme()"),
paste0("\n", paste0(rep("\u2500", 41+nchar(name)), collapse="")), "\n")
if (logtrans) cat("Response: log(PK)\n") else cat("\nResponse: logPK\n")
print(anova(modB), digits = 6, signif.stars = FALSE)
cat("\ntreatment T \u2013 R:\n")
print(signif(EMA.B$tTable["treatmentT", c(1:2, 4:5)], 5))
cat(DF, "Degrees of Freedom (equivalent to SAS\u2019 DDFM=CONTAIN)\n\n")
}
} else { # by lmer/lmerTest (Satterthwaite's or Kenward-Roger DF)
if (logtrans) { # use the raw data and log-transform internally
modB <- lmer(log(PK) ~ sequence + period + treatment + (1|subject),
data = ret$data)
} else { # use the already log-transformed data
modB <- lmer(logPK ~ sequence + period + treatment + (1|subject),
data = ret$data)
}
if (option == 1) {
EMA.B <- summary(modB, ddf="Satterthwaite")
} else {
EMA.B <- summary(modB, ddf="Kenward-Roger")
}
PE <- EMA.B$coefficients["treatmentT", "Estimate"]
CI <- exp(PE + c(-1, +1) *
qt(1-alpha, EMA.B$coef["treatmentT", "df"]) *
EMA.B$coef["treatmentT", "Std. Error"])
DF <- EMA.B$coefficients["treatmentT", "df"]
if (verbose) {
if (option == 1) {
cat("\nData set", paste0(file, set, ": Method B (option = 1) by lmer()"),
paste0("\n", paste0(rep("\u2500", 46), collapse="")), "\n")
if (logtrans) cat("Response: log(PK)\n") else cat("\nResponse: logPK\n")
print(anova(modB, ddf="Satterthwaite"), digits=6, signif.stars=FALSE)
cat("\ntreatment T \u2013 R:\n")
print(signif(EMA.B$coefficients["treatmentT", c(1:2, 4:5)], 5))
cat(signif(DF, 6), "Degrees of Freedom (equivalent to SAS\u2019 DDFM=SATTERTHWAITE)\n\n")
} else {
df.txt <- "3; equivalent to Stata\u2019s dfm=Kenward Roger EIM)"
cat("\nData set", paste0(file, set, ": Method B (option = 3) by lmer()"),
paste0("\n", paste0(rep("\u2500", 46), collapse="")), "\n")
if (logtrans) cat("Response: log(PK)\n") else cat("\nResponse: logPK\n")
print(anova(modB, ddf="Kenward-Roger"), digits=6, signif.stars=FALSE)
cat("\ntreatment T \u2013 R:\n")
print(signif(EMA.B$coefficients["treatmentT", c(1:2, 4:5)], 5))
cat(signif(DF, 6), "Degrees of Freedom (equivalent to Stata\u2019s dfm=Kenward Roger EIM)\n\n")
}
}
} # end of evaluation by option=2 (lme) or option=1/3 (lmer)
PE <- exp(PE)
res <- data.frame(ret$type, "A", ret$n, ret$nTT, ret$nRR,
paste0(ret$Sub.Seq, collapse="|"),
paste0(ret$Miss.seq, collapse="|"),
paste0(ret$Miss.per, collapse="|"), alpha,
DF, ret$CVwT, ret$CVwR, ret$swT, ret$swR, ret$sw.ratio,
ret$sw.ratio.upper, ret$BE1, ret$BE2, CI[1], CI[2],
PE, "fail", "fail", "fail", log(CI[2])-log(PE),
paste0(ret$ol, collapse="|"), ret$CVwR.rec,
ret$swR.rec, ret$sw.ratio.rec,
ret$sw.ratio.rec.upper, ret$BE.rec1,
ret$BE.rec2, "fail", "fail", "fail",
stringsAsFactors=FALSE)
names(res)<- c("Design", "Method", "n", "nTT", "nRR", "Sub/seq",
"Miss/seq", "Miss/per", "alpha", "DF", "CVwT(%)",
"CVwR(%)", "swT", "swR", "sw.ratio", "sw.ratio.CL",
"L(%)", "U(%)", "CL.lo(%)", "CL.hi(%)", "PE(%)",
"CI", "GMR", "BE", "log.half-width", "outlier",
"CVwR.rec(%)", "swR.rec", "sw.ratio.rec",
"sw.ratio.rec.CL", "L.rec(%)", "U.rec(%)",
"CI.rec", "GMR.rec", "BE.rec")
if (ret$BE2 == 1.25) { # change column names if not scaling
colnames(res)[which(names(res) == "L(%)")] <- "BE.lo(%)"
colnames(res)[which(names(res) == "U(%)")] <- "BE.hi(%)"
}
if (regulator == "HC" & alpha == 0.5) {
if (option == 2)
stop("Please use option = 1 or option = 3 instead.")
res[which(names(res) == "L(%)")] <- 0.80
res[which(names(res) == "U(%)")] <- 1.25
colnames(res)[which(names(res) == "L(%)")] <- "BE.lo(%)"
colnames(res)[which(names(res) == "U(%)")] <- "BE.hi(%)"
}
if (!is.na(ret$BE.rec2) & ret$BE.rec2 == 1.25) { # change column names if not scaling
colnames(res)[which(names(res) == "L.rec(%)")] <- "BE.rec.lo(%)"
colnames(res)[which(names(res) == "U.rec(%)")] <- "BE.rec.hi(%)"
}
# Convert CVs, limits, PE, and CI (up to here as fractions) to percent
res$"CVwT(%)" <- 100*res$"CVwT(%)"
res$"CVwR(%)" <- 100*res$"CVwR(%)"
if ("BE.lo(%)" %in% names(res)) { # conventional limits
res$"BE.lo(%)" <- 100*res$"BE.lo(%)"
res$"BE.hi(%)" <- 100*res$"BE.hi(%)"
} else { # expanded limits
res$"L(%)" <- 100*res$"L(%)"
res$"U(%)" <- 100*res$"U(%)"
}
if (!is.na(res$"CVwR.rec(%)")) { # only if recalculated CVwR
res$"CVwR.rec(%)" <- 100*res$"CVwR.rec(%)"
if ("BE.rec.lo(%)" %in% names(res)) { # conventional limits
res$"BE.rec.lo(%)" <- 100*res$"BE.rec.lo(%)"
res$"BE.rec.hi(%)" <- 100*res$"BE.rec.hi(%)"
} else { # expanded limits
res$"L.rec(%)" <- 100*res$"L.rec(%)"
res$"U.rec(%)" <- 100*res$"U.rec(%)"
}
}
res$"PE(%)" <- 100*res$"PE(%)"
res$"CL.lo(%)" <- 100*res$"CL.lo(%)"
res$"CL.hi(%)" <- 100*res$"CL.hi(%)"
if (round(res$"CL.lo(%)", 2) >= 100*ret$BE1 &
round(res$"CL.hi(%)", 2) <= 100*ret$BE2)
res$CI <- "pass" # CI within acceptance range
if (!regulator == "HC") {
if (round(res$"PE(%)", 2) >= 80 & round(res[["PE(%)"]], 2) <= 125)
res$GMR <- "pass" # PE within 80.00-125.00%
} else {
if (round(res$"PE(%)", 1) >= 80 & round(res[["PE(%)"]], 1) <= 125)
res$GMR <- "pass" # PE within 80.00-125.00%
}
if (res$CI == "pass" & res$GMR == "pass")
res$BE <- "pass" # if passing both, conclude BE
if (!is.na(res$"CVwR.rec(%)")) {
if (round(res$"CL.lo(%)", 2) >= 100*ret$BE.rec1 &
round(res$"CL.hi(%)", 2) <= 100*ret$BE.rec2)
res$CI.rec <- "pass" # CI within acceptance range
res$GMR.rec <- res$GMR
if (res$CI.rec == "pass" & res$GMR.rec == "pass")
res$BE.rec <- "pass" # if passing both, conclude BE
}
if (details) { # results in full precision
ret <- res
if (as.character(res$outlier) == "NA") {
# remove superfluous columns if ola=FALSE or ola=TRUE
# and no outlier(s) detected
ret <- ret[, !names(ret) %in% c("outlier", "CVwR.rec(%)",
"swR.rec", "sw.ratio.rec",
"L.rec(%)", "U.rec(%)",
"CI.rec", "GMR.rec", "BE.rec",
"sw.ratio.rec.CL")]
}
if (regulator == "HC" & alpha == 0.5) {
# remove columns for Health Canada (PE of Cmax within 80.0-125.0)
ret <- ret[, !names(ret) %in% c("CL.lo(%)", "CL.hi(%)", "CI", "BE",
"log.half-width")]
}
#class(ret) <- "repBE"
return(ret)
}
# Round percents to two decimals according to the GL
res$"CVwT(%)" <- round(res$"CVwT(%)", 2)
res$"CVwR(%)" <- round(res$"CVwR(%)", 2)
if ("BE.lo(%)" %in% names(res)) { # conventional limits
res$"BE.lo(%)" <- round(res$"BE.lo(%)", 2)
res$"BE.hi(%)" <- round(res$"BE.hi(%)", 2)
} else { # expanded limits
res$"L(%)" <- round(res$"L(%)", 2)
res$"U(%)" <- round(res$"U(%)", 2)
}
res$"PE(%)" <- round(res$"PE(%)", 2)
res$"CL.lo(%)" <- round(res$"CL.lo(%)", 2)
res$"CL.hi(%)" <- round(res$"CL.hi(%)", 2)
if (!is.na(res$"CVwR.rec(%)")) { # only if recalculated CVwR
res$"CVwR.rec(%)" <- round(res$"CVwR.rec(%)", 2)
if ("BE.rec.lo(%)" %in% names(res)) { # conventional limits
res$"BE.rec.lo(%)" <- round(res$"BE.rec.lo(%)", 2)
res$"BE.rec.hi(%)" <- round(res$"BE.rec.hi(%)", 2)
} else { # expanded limits
res$"L.rec(%)" <- round(res$"L.rec(%)", 2)
res$"U.rec(%)" <- round(res$"U.rec(%)", 2)
}
}
overwrite <- TRUE # default
if (print) { # to file in UTF-8
if (ask & file.exists(results)) {
answer <- tolower(readline("Results already exists. Overwrite the file [y|n]? "))
if(answer != "y") overwrite <- FALSE
}
if (overwrite) { # either the file does not exist or should be overwritten
# only binary mode supports UTF-8 and different line endings
res.file <- file(description=results, open="wb")
res.str <- info.env(fun="method.B", option=option, path.in, path.out,
file, set, ext, exec, data)
if (os == "Windows") res.str <- gsub("\n", "\r\n", res.str) # CRLF
if (os == "Darwin") res.str <- gsub("\n", "\r", res.str) # CR
writeBin(charToRaw(res.str), res.file)
close(res.file)
}
}
# insert DF of treatment difference and alpha into the text generated by CV.calc()
cut.pos <- unlist(gregexpr(pattern="Regulator", ret$txt))
left.str <- substr(ret$txt, 1, cut.pos-1)
right.str <- substr(ret$txt, cut.pos, nchar(ret$txt))
insert.str <- "Degrees of freedom : "
if (option == 1) insert.str <- paste0(insert.str, sprintf("%7.3f (Satterthwaite)", DF))
if (option == 2) insert.str <- paste0(insert.str, sprintf("%3i", DF))
if (option == 3) insert.str <- paste0(insert.str, sprintf("%7.3f (Kenward-Roger)", DF))
insert.str <- paste0(insert.str, "\nalpha : ", alpha)
if (!alpha == 0.5) { # Health Canada Cmax
insert.str <- paste0(insert.str," (", 100*(1-2*alpha), "% CI)\n")
} else {
insert.str <- paste0(insert.str," (only PE assessed)\n")
}
txt <- paste0(left.str, insert.str, right.str)
if (!is.na(res$"CVwR.rec(%)")) {
txt1 <- paste0("\n\nAssessment based on original CVwR",
sprintf(" %.2f%%", res$"CVwR(%)"))
txt <- paste0(txt, txt1, "\n", paste0(rep("\u2500", nchar(txt1)-2), collapse=""))
}
if (!alpha == 0.5) {
txt <- paste0(txt,
"\nConfidence interval: ", sprintf("%6.2f%% ... %6.2f%%",
res$"CL.lo(%)", res$"CL.hi(%)"),
" ", res$CI,
"\nPoint estimate : ", sprintf("%6.2f%%", res$"PE(%)"),
" ", res$GMR,
"\nMixed (CI & PE) : ", res$BE, "\n")
} else {
txt <- paste0(txt,
"\nPoint estimate : ", sprintf("%6.2f%%", res$"PE(%)"),
" ", res$GMR, "\n")
}
txt <- paste0(txt,
repBE.draw.line(called.from="ABEL", L=ret$BE1, U=ret$BE2,
lo=CI[1], hi=CI[2], PE=PE), "\n")
if (!is.na(res$"CVwR.rec(%)")) {
txt1 <- paste0("\nAssessment based on recalculated CVwR",
sprintf(" %.2f%%", res$"CVwR.rec(%)"))
txt <- paste0(txt, txt1, "\n", paste0(rep("\u2500", nchar(txt1)-1), collapse=""))
txt <- paste0(txt, "\nConfidence interval: ", res$CI.rec,
"\nPoint estimate : ", res$GMR.rec,
"\nMixed (CI & PE) : ", res$BE.rec, "\n")
txt <- paste0(txt,
repBE.draw.line(called.from="ABEL", L=ret$BE.rec1, U=ret$BE.rec2,
lo=CI[1], hi=CI[2], PE=PE), "\n")
}
if (res$Design == "TRR|RTR")
txt <- paste0(txt, "Note: The extra-reference design assumes lacking period effects. ",
"The treatment\ncomparison will be biased in the presence of a ",
"true period effect.\n")
if (print & overwrite) {
res.file <- file(results, open="ab") # line endings
res.str <- txt # LF (UNIXes, Solaris)
if (os == "Windows") res.str <- gsub("\n", "\r\n", res.str) # CRLF (Windows)
if (os == "Darwin") res.str <- gsub("\n", "\r", res.str) # CR (OSX)
writeBin(charToRaw(res.str), res.file)
close(res.file)
}
} # end of function method.B()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.