Bioequivalence and Bioavailability Forum

Hi Nastia,

❝ […] replicateBE was the first thing that I've thought of, but I couldn't dig up the code for get.data in order to modify it for my needs.

I would not try to use this function. It calls others which are not exported. Not worth the efforts to modify it.

❝ […] what is the difference between Dataset<-rds01 and Dataset<-read_excel("rds01.xlsx", sheet = 1)? The initial dataset is the same, but the results of the algo are different. So before using we should somehow prepare the data, but how to do it?

Duno how you accessed the datasets. Try this one:

library(replicateBE)
library(readxl)
library(nlme)
path <- "your path"
file <- "rds01.xlsx"
name <- paste0(path, "/", file)
DS1  <- as.data.frame(read_excel(path = name, sheet = 1,
                                 na = c("NA", "ND", ".", "", "Missing"),
                                 skip = 0, col_names = TRUE))
str(DS1) # show the structure: data.frame!
cols <- c("subject", "period", "sequence", "treatment")
DS1[cols] <- lapply(DS1[cols], factor) # factorize
DS2  <- rds01
str(DS2) # show the structure: named S3-object, factorized data.frame
# add the groups as factors
DS1$group <- factor(ifelse(as.numeric(levels(DS1$subject))[DS1$subject]                            < 31, 1, 2))
DS2$group <- factor(ifelse(as.numeric(levels(DS2$subject))[DS2$subject]                            < 31, 1, 2))
res <- data.frame(origin = c("Excel", "replicateBE"),
                  PE = NA, CL.lo = NA, CL.hi = NA)
ow <- options("contrasts") # save options
options(contrasts = c("contr.treatment", "contr.poly"))
on.exit(ow)                # reset options if an error occurs
for (j in 1:nrow(res)) {
  if (j == 1) data = DS1 else data <- DS2
  modB        <- lme(log(PK) ~ sequence + group + sequence:group + period +
                               period%in%group + treatment,
                               random = ~1 | subject,
                               na.action = na.omit, data = data)
  EMA.B       <- summary(modB)
  PE          <- EMA.B$tTable["treatmentT", "Value"]
  res[j, 3:4] <- 100*exp(PE + c(-1, +1) *
                         qt(1 - 0.05, EMA.B$tTable["treatmentT", "DF"]) *
                         EMA.B$tTable["treatmentT", "Std.Error"])
  res$PE[j]   <- 100*exp(PE)
}
print(res, row.names = FALSE) # should be identical

#      origin       PE    CL.lo   CL.hi
#       Excel 115.7275 107.1136 125.034
# replicateBE 115.7275 107.1136 125.034

---

Edit: Hey Mittyri, you were faster!