R/test2x2.R

In BE: Bioequivalence Study Data Analysis

test2x2 = function(bedata, Var)
{
  if(!assert(bedata)) {
    cat("\n Subject count should be balanced!\n");
    return(NULL);
  }
  bedata = TrimData(bedata)
  bedata = af(bedata, c("GRP", "SUBJ", "PRD", "TRT"))

  TrtLevel = levels(bedata$TRT)
  if (!setequal(TrtLevel, c("R", "T"))) stop("TRT column should have only 'R' and 'T'!")

  Formula = as.formula(paste0("log(`", Var, "`) ~ GRP/SUBJ + PRD + TRT"))

# Type III sums of squares as SAS PROC GLM. GROUP (sequence) is tested against
# SUBJECT(GROUP); every other term is tested against the residual mean square.
  aovT3 = T3test(Formula, bedata, E="GRP:SUBJ")
  Grp   = aovT3[["GRP"]]
  Rest  = aovT3[["The Rest Terms"]]

  mse = Rest["RESIDUALS", "Mean Sq"]
  dfe = Rest["RESIDUALS", "Df"]

  dfs = Grp["GRP", "Df"] + Rest["GRP:SUBJ", "Df"]
  sss = Grp["GRP", "Sum Sq"] + Rest["GRP:SUBJ", "Sum Sq"]
  mss = sss/dfs
  SUBJ = c(dfs, sss, mss, mss/mse, 1 - pf(mss/mse, dfs, dfe))

  y = model.frame(Formula, bedata)[, 1]
  TOT = c(length(y) - 1, sum((y - mean(y))^2), NA, NA, NA)

  ANOVA = rbind(SUBJ, Grp["GRP", ], Rest["GRP:SUBJ", ], Rest["PRD", ],
                Rest["TRT", ], Rest["RESIDUALS", ], TOT)
  dimnames(ANOVA) = list(c("SUBJECT", "GROUP", "SUBJECT(GROUP)", "PERIOD", "DRUG",
                           "ERROR", "TOTAL"),
                         c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)"))
  class(ANOVA) = "anova"

  sig2b = max(0, (Rest["GRP:SUBJ", "Mean Sq"] - mse)/2) # See pp 62-63 of Chow & Liu
  sig2w = mse
  cvs = rbind(cbind(sig2b, sig2w), cbind(mse2cv(sig2b), mse2cv(sig2w)))
  dimnames(cvs) = list(cbind("Variance Estimate", "Coefficient of Variation, CV(%)"),
                       cbind("Between Subject", "Within Subject"))

  lsmeans = LSM(Formula, bedata, "TRT")
  lsm = cbind(exp(lsmeans["R", "LSmean"]), exp(lsmeans["T", "LSmean"]))
  dimnames(lsm) = list("Geometric Means", cbind("Reference Drug", "Test Drug"))

  Contrast = ifelse(TrtLevel == "T", 1, -1) # Test - Reference
  ci0 = CIest(Formula, bedata, "TRT", Contrast, conf.level=0.90)
  pe = ci0[1, "Estimate"]
  ci = exp(ci0[1, c("Lower CL", "Estimate", "Upper CL"), drop=FALSE])
  dimnames(ci) = list("90% CI for Ratio", c("Lower Limit", "Point Estimate", "Upper Limit"))

  ss = cbind(ssmse(mse), ssmse(mse, True.R=exp(pe)))
  dimnames(ss) = list("80% Power Sample Size", c("True Ratio=1", "True Ratio=Point Estimate"))

  result = list(ANOVA, cvs, lsm, ci, ss);
  names(result) = c("Analysis of Variance (log scale)", "Between and Within Subject Variability", "Least Square Means (geometric mean)", "90% Confidence Interval of Geometric Mean Ratio (T/R)", "Sample Size")

  return(result);
}