Back to the Top
Hi everyone,
While I was looking for a quick, easy, and free way of running the
typical bioequivalence ANOVA, I stumbled upon the following very
useful post:
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/16552.html
"R", similar to S and S+, is an open source statistical package that
has gained much utility in academia (and industry), although it
promises no validation or guarantees whatsoever. Weeks ago I didn't
know R existed, and now I am its newest fan. R is available for
download for most computer platforms here: http://www.r-project.org/
To take you through how to run a non-linear mixed model in R (similar
to PROC MIXED in SAS), consider the example 2-way crossover data
provided in the appendix of the TPD's Part A guidance for
bioequivalence:
Subject Sequence Period Treatment AUCt
1 AB 1 A 364.74595
1 AB 2 B 375.426
2 BA 1 B 595.03875
2 BA 2 A 404.9456
3 BA 1 B 471.1644
3 BA 2 A 702.8314
4 AB 1 A 233.2503
4 AB 2 B 190.39375
5 BA 1 B 257.45585
5 BA 2 A 247.41105
6 AB 1 A 178.19155
6 AB 2 B 175.37495
7 BA 1 B 381.824
7 BA 2 A 246.3878
8 AB 1 A 407.99185
8 AB 2 B 360.8275
9 BA 1 B 218.47035
9 BA 2 A 315.4761
10 AB 1 A 140.1254
10 AB 2 B 91.80655
11 AB 1 A 165.365
11 AB 2 B 269.01575
12 BA 1 B 105.5625
12 BA 2 A 87.9882
13 BA 1 B 290.142
13 BA 2 A 182.77325
14 AB 1 A 122.47815
14 AB 2 B 230.48885
15 BA 1 B 143.5487
15 BA 2 A 67.9815
16 AB 1 A 274.5791
16 AB 2 B 344.4828
Highlight and copy the entire range of data (from here, or Excel, or
another program), including the headers. In R Console window, type the
following command, which will import the data from the clipboard into
a data set with the name of your choosing (say, TPDdata):
TPDdata <- read.table(file="clipboard", header=TRUE)
To see the data you've imported, simply type "TPDdata" and it will
display the dataset.
Now to run the ANOVA on ln(AUCt), copy and paste the following lines
into the R Console window:
library(nlme)
lm.lnauct <- lme(-log(AUCt)~Treatment+Period+Sequence, data=TPDdata,
random=~1|Subject,na.action=na.exclude)
summary(lm.lnauct)
intervals(lm.lnauct,0.90)
VarCorr(lm.lnauct)
The ANOVA is run on the -log because the default comparison is
actually Treatment B/A rather than A/B, and the negative will provide
the results the way we are used to seeing them, Test/Reference. Just
be careful in interpreting estimates as they will be negative when in
fact they should be positive. The code "random=~1|Subject" will treat
Subject(Sequence) as a random effect, and na.action=exclude will throw
out missing data from the analysis. Please note that in R, missing
data is indicated with an "NA" rather than a "." like in SAS.
The command "summary(lm.auct)" will provide an assortment of output,
including the p-values for each effect (treatment, period, and
sequence).
The command "intervals(lm.lnauct,0.90)" will calculate the 90%
confidence intervals for ln-AUCt, providing the following output:
Approximate 90% confidence intervals
Fixed effects:
lower est. upper
(Intercept) -5.8653512 -5.43223224 -4.99911330
TreatmentB -0.2992754 -0.13105891 0.03715754
Period -0.1133948 0.05482169 0.22303814
SequenceBA -0.5652860 -0.08186053 0.40156492
attr(,"label")
[1] "Fixed effects:"
Random Effects:
Level: Subject
lower est. upper
sd((Intercept)) 0.3603975 0.5146334 0.7348761
Within-group standard error:
lower est. upper
0.1979603 0.2701330 0.3686185
Remember that it's on the ln-scale. The 90% confidence intervals
transformed to the normal scale will be the lower and upper intervals
calculated on the TreatmentB line: exp(-0.2992754) to exp(0.0371574),
or 74.13 - 103.79% (failing to demonstrate bioequivalence). The point
estimate will be exp(-0.13105891) = 87.71%. These match the results
reported by the TPD Part A guidance.
The command "VarCorr(lm.lnauct)" produces the following output:
Subject = pdLogChol(1)
Variance StdDev
(Intercept) 0.26484753 0.5146334
Residual 0.07297183 0.2701330
The intrasubject variability is calculated by sqrt(exp(MSR)-1) =
sqrt(exp(0.07297183)-1) = 27.5%. The intersubject variability is
calculated by sqrt(exp(0.26484753)-1) = 55.1%. The TPD guidance uses
slightly different formula to calculate Intra and InterCV, however the
numbers are comparable (27% and 51%). If you wish to use the TPD's
formulae (they use sqrt(MSR) as the intraCV, etc.), the results in R
will produce the same results.
Although this program wouldn't be suitable for a regulatory drug
submission, it provides a quick answer (1-line ANOVA!) that matches
PROC MIXED and could be a useful educational tool (reading between the
lines, I wanted my students to be able to run the ANOVA without having
to buy SAS).
I've run the same data through PROC MIXED in SAS using the following
code (adapted from Scott Patterson and Byron Jones's new book,
"Bioequivalence and Statistics in Clinical Pharmacology" - a great
read by the way):
data TPDdata;
input Subject Sequence$ Period Treatment$ AUCt;
logauct=log(AUCt);
datalines;
1 AB 1 A 364.74595
...
;
run;
proc mixed data=TPDdata;
class Sequence Subject Period Treatment;
model logauct=Sequence Period Treatment/ddfm=kenwardroger;
random Subject(Sequence);
lsmeans Treatment/pdiff cl alpha=0.1;
estimate 'Average Bioequivalence for log(AUCt)' Treatment 1 -1;
run;
The results match R.
The one thing I'm missing in the R analysis is how to extract the
least square means from the ANOVA output for the separate treatments.
If anyone figures this out, please let me know.
Best,
-David Dubins
--
David Dubins, Ph.D., B.A.Sc.
Global Bioequivalence Consulting
Assistant Professor, Leslie Dan Faculty of Pharmacy
University of Toronto
Back to the Top
The following message was posted to: PharmPK
Dear Dave,
I am a fan of R for a few years (I made almost complete PhD in R) and
I share your opinion about R.
It is very fine resource. I can not help you with the problem you
mentioned but you are probably
aware of good reference manual. If you use some of the add-in packages
try to contact the developer
by e-mail. There are also some other possibilities for advanced users.
Enjoy using R.
Zeljko Debeljak, PhD
Medical Biochemistry Specialist
Osijek Clinical Hospital
CROATIA
Want to post a follow-up message on this topic?
If this link does not work with your browser send a follow-up message to PharmPK@boomer.org with "Average Bioequivalence ANOVA in R" as the subject | Support PharmPK by using the |
Copyright 1995-2011 David W. A. Bourne (david@boomer.org)