# R/test2x2.R In BE: Bioequivalence Study Data Analysis

#### Documented in test2x2

```test2x2 = function(bedata, Var)
{
if(!assert(bedata)) {
cat("\n Subject count should be balanced!\n");
return(NULL);
}
bedata = TrimData(bedata)

Yijk = log(bedata[, Var])
Yi11 = log(bedata[bedata\$GRP=="RT" & bedata\$PRD==1, Var])
Yi21 = log(bedata[bedata\$GRP=="RT" & bedata\$PRD==2, Var])
Yi12 = log(bedata[bedata\$GRP=="TR" & bedata\$PRD==1, Var])
Yi22 = log(bedata[bedata\$GRP=="TR" & bedata\$PRD==2, Var])

n1 = length(Yi11)
n2 = length(Yi12)

Y... = mean(Yijk)
SStotal = sum((Yijk - Y...)^2)

Y.11 = mean(Yi11)
Y.21 = mean(Yi21)
Y.12 = mean(Yi12)
Y.22 = mean(Yi22)

Yi.1 = (Yi11 + Yi21) / 2
Yi.2 = (Yi12 + Yi22) / 2

Y..1 = mean(Yi.1)
Y..2 = mean(Yi.2)

mu.r = (Y.11 + Y.22) / 2
mu.t = (Y.21 + Y.12) / 2

di1 = (Yi21 - Yi11) / 2
di2 = (Yi22 - Yi12) / 2

d.1 = mean(di1)
d.2 = mean(di2)

#  Chat = Y.12 + Y.22 - Y.11 - Y.21
#  Fhat = mu.t - mu.r
#  Phat = (Y.21 - Y.11 - Y.12 + Y.22)/2

SScarry   = 2*n1*n2/(n1 + n2)*(Y.12 + Y.22 - Y.11 - Y.21)^2 / 4
SSinter   = (sum((Yi.1 - Y..1)^2) + sum((Yi.2 - Y..2)^2)) * 2
SSbetween = SScarry + SSinter

SSperiod  = 2*n1*n2/(n1+n2)*(Y.21 + Y.22 - Y.11 - Y.12)^2 / 4
SSdrug    = 2*n1*n2/(n1+n2)*(Y.21 + Y.12 - Y.11 - Y.22)^2 / 4

SSintra   = 2*(sum((di1 - d.1)^2) + sum((di2 - d.2)^2))
#  SSmodel = SStotal - SSintra

Source = c("SUBJECT", "GROUP", "SUBJECT(GROUP)", "PERIOD", "DRUG", "ERROR", "TOTAL");
SS     = c(SSbetween, SScarry, SSinter, SSperiod, SSdrug, SSintra, SStotal);
DF     = c(n1 + n2 - 1, 1, n1 + n2 - 2, 1, 1, n1 + n2 - 2, 2*n1 + 2*n2 - 1);
MS     = SS / DF
mse    = SSintra/(n1 + n2 - 2)
Fv      = MS / c(mse, MS[3], mse, mse, mse, mse, mse);
p1 = 1 - pf(Fv[1], n1 + n2 - 1, n1 + n2 - 2)
p2 = 1 - pf(Fv[2], 1, n1 + n2 - 2);
p3 = 1 - pf(Fv[3], n1 + n2 - 2, n1 + n2 - 2);
p4 = 1 - pf(Fv[4], 1, n1 + n2 - 2);
p5 = 1 - pf(Fv[5], 1, n1 + n2 - 2);
p  = c(p1, p2, p3, p4, p5, NA, NA)
Fv[6] = Fv[7] = MS[7] = NA

ANOVA = cbind(SS, DF, MS, Fv, p)
dimnames(ANOVA) = list(Source,c("Sum Sq", "Df", "Mean Sq", "F value", "Pr(>F)"))
class(ANOVA) = "anova"

pe = mu.t - mu.r
se = sqrt(mse/2 * (1/n1 + 1/n2))   # See pp 62-63 of Chow & Liu
t0 = qt(0.95, n1 + n2 - 2);
ci0 = cbind(pe - t0 * se, pe, pe + t0 * se)

sig2b = max(0, (MS[3] - MS[6])/2)
sig2w = MS[6]

cvs = cbind(sqrt(exp(sig2b) - 1), sqrt(exp(sig2w) - 1)) * 100

lsm = cbind(exp(mu.r), exp(mu.t))
dimnames(lsm) = list("Geometric Means", cbind("Reference Drug", "Test Drug"))

ci = exp(ci0);
dimnames(ci) = list("90% CI for Ratio", c("Lower Limit", "Point Estimate", "Upper Limit"));

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

cvs = rbind(cbind(sig2b, sig2w), cvs)
dimnames(cvs) = list(cbind("Variance Estimate", "Coefficient of Variation, CV(%)"), cbind("Between Subject", "Within Subject"))

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);
}
```

## Try the BE package in your browser

Any scripts or data that you put into this service are public.

BE documentation built on Jan. 7, 2023, 1:11 a.m.