Bioequivalence and Bioavailability Forum

Hi BEQool,

❝ Would anyone know what are the equations to calculate standard deviation of a study for different study designs (2x2, 2x3x3, 2x2x3, 2x2x4)?

❝ For example, standard error of a study for a given parameter is given in Phoenix Winnonlin with "Diff_SE" in the "Average Bioequivalence" tab but I cannot either find or calculate its standard deviation.

❝ I assume it has something to do with coefficients and number of subjects in sequences like in known.designs() but I cant get my head around it

If I am not mistaken, you are looking for the standard error of the treatment difference, which is used to construct a 1-2*alpha Confidence interval. I agree, this is not all that easy.
One thing is to apply an easy equation from Chow & Liu which works well in the case of a simple 222BE design, but it can get more obscure if the design gets a little complex.

For a general approach, you can derive the covariance matrix, V, for the fixed effects. If X is your model matrix (=design matrix) and you have fitted your model and extracted your MSE from that model, then you simple have
V=MSE*(X'X)^-1

If treatment A is indicated by the column a in X, and treatment B is in the column b of X, then V[a,a] is the SE of A and V[b,b] is the SE of B. V[a,b]=V[b,a] is the covariance of the two, hence:
SE_diff=SE_A+SE_B-2Cov_AB

If you program in R here is a little snippet for you to play with, and it may be handy for you to have the emmeans package installed.
R turns out to have a very handy function called vcov which directly gives us V without us having to do the underlying matrix algebra. I hope the following is illustrative:


library(emmeans)
set.seed(148923)
Seq=rep(sample(rep(c("RT", "TR"), 5)),2)
Per=c(rep(1, 10), rep(2,10))
Trt=substr(Seq,Per,Per)
Subj=c(rep(c(1:10),2))
lnCmax=runif(20, 100,150)
data.frame(Subj, Per, Seq, Trt, lnCmax)  #just some invented data for a 222BE trial
M=lm(lnCmax~0+Trt+factor(Per)+Seq+factor(Subj)) #I am fitting without intercept and Trt first!
V=vcov(M)
V[1:2, 1:2]  #show me the first two columns/rows
SEd=sqrt(V[1,1] + V[2,2] - 2*V[1,2])
confint(pairs(emmeans(M, "Trt"), reverse =T))
# Is the SEd the same as the SE from the confint/pairs object?
SEd
# Yes it is. Perhaps ElMaestro occasionally is right after all?


You can extend this as necessary.
Be careful about intercepts and the order of factors in your lm.