Back to the Top
Colleagues
I am attempting to use R to replicate FDA’s code for reference-scaled bioequivalence, based on the
SAS code provided in the Draft Guidance on Progesterone (Progesterone_caps_19781_RC02-11.pdf)
I am able to generate results nearly identical to the SAS output (generated by a colleague) with one
exception. In the calculation of dout1 for dat, the SAS code is:
class seq
model dlat=seq/ddfm=satterth
estimate 'average' intercept 1 seq 0.5 0.5/e cl alpha=0.1;
ods output CovParms=dout1;
I attempted to replicate this with:
require(“lme4”)
lmer(DLAT ~ 0 + offset(rep(0, nrow(DATA))) + 1 | SEQ, data=DATA, REML=T)
where SEQ is a factor.
It is possible that my formula is entirely bogus — the intent was to force an intercept of 1 by
fixing the intercept to zero and adding an offset of 1.
The output from R does not provide anything named anything like CovParms. There are several
variance estimates, one of which is within 1% of the value reported by SAS.
Is anyone sufficiently familiar with SAS to recommend the correct implementation in R. I am not wed
to lmer.
Any assistance will be greatly appreciate.
Dennis
Dennis Fisher MD
Back to the Top
I am out of office so can not check, but if you need intercept 1, I would use
1+y ~ -1+ ...
One on the left will give you intercept 1 (or you may need to rename 1+y as Y1) and -1 on the right
will tell R that intercept is omitted (or equal to 0)
Leonid
Back to the Top
Dennis,
You're making this way more complicated than it is. The sas code is a simple fixed effect one way
analysis of variance and the CovParm output will output the Residual variance.
In R there is no need the use lmer(), lm() would works just fine here:
aov1<-lm(DLAT ~ SEQ, DATA)
sum.aov1<-summary(aov1)
residual<-aov1$sigma**2 ## will equal estimate in dout2 from sas
It seems you may be misinterpreting the 1 next to the intercept, there is no fixing of the intercept
here, the 1 is simply a coefficient to multiply the intercept by as in:
Average = 1*intercept + 0.5*seq1 + 0.5*seq2
I think all you need from this statement is the associated denominator degree of freedom.
dfd<-aov1$df.residual
I don't believe the satterth option in sas will ever have an effect on the one-way fixed effect
anova.
--Joe Boucher
Back to the Top
Correction: there was a typo in the R code that was sent yesterday
To fix intercept of the y vs x dependence to 1 one Can use
y1 ~ -1+ x
Where y1= y-1, and -1 on the right will tell R that intercept is omitted (or equal to 0)
Leonid
Back to the Top
Leonid
Thanks for the suggestion. I tried it (and a bunch of variants) without success — output does not
match SAS.
Dennis
Dennis Fisher MD
P < (The "P Less Than" Company)
Back to the Top
Joe
Thanks so much. With a minor modification, that appears to work:
residual<-aov1$sigma**2 ## will equal estimate in dout2 from sas
should be:
residual <- sum.aov1$sigma**2 ## i.e., sum.aov1, not aov1
Dennis
Dennis Fisher MD
P < (The "P Less Than" Company)
Back to the Top
Dear all,
I'm sorry, it's not that easy to reproduce FDA's SAS-code in R. /Some/
parts are possible but one has to reproduce the code in its /entirety/.
Please stop posting half-backed solutions.
If you think that you have a working code, start testing it with data
sets from the public domain (e.g. the ones in EMA's Q&A document,
Section 10:
http://www.ema.europa.eu/docs/en_GB/document_library/Scientific_guideline/2009/09/WC500002963.pdf).
Results should agree for all designs, balanced or imbalanced sequences,
complete and incomplete data. AFAIK, such an agreement was only shown
for Phoenix/WinNonlin (http://dx.doi.org/10.13140/RG.2.1.1008.0800). I
would be happy to see a working code in R!
The devil is in the details.
Package nlme, function lme():
Allows to specify the variance-covariance-structure, but /not/ SAS'
FA0(2) (i.e. Banded No-Diagonal Factor Analytic f=2). nlme is not
developed any more; little hope for a solution.
Package lme4, function lmer():
Does not allow to specify weights at all. According to the man-pages the
authors do not intend to implement them.
Given that I don't see how we could reproduce FDA's code in R right now.
Remember that if reference-scaling is /not/ allowed (Swr <0.294) one has
to use FDA's model for ABE (also given in the 2001 guidance, Appendix
E). Quote: "TYPE=FA0(2) could possibly be replaced by TYPE=CSH". This
structure is supported by lme(). Problems start with the partial
replicate design (TRR | RTR | RRT). FDA's model is over-specified (since
the test-product is /not/ repeated) for this design. Therefore, one
cannot get a meaningful estimate for the within-subject variance for the
test-product. Yes, you get a number in most cases but it is nonsense (by
definition). In the worst case you don't get convergence and the
optimizer stops (doesn't matter whether you specify FA0(2) or CSH). You
invested a lot of money in a BE study and are not able to evaluate it –
independent from the software used! IIRC I have at least three such data
sets in my files. One of my clients send a request to the FDA and
/never/ got an answer…
BTW, FA0(1) in SAS (or its analogue in Phoenix/WinNonlin) leads to
convergence in /all/ cases I have seen so far. My suggestion (for
partial replicate designs) is to state in the SAP that you intend to
deviate from FDA's code or – IMHO better – opt for a fully replicated
3-period design (i.e. TRT | RTR) instead.
Another story. In none of FDA's product-specific guidances a minimum
sample size is mentioned. However, at the August 5, 2009 ACPS meeting
Dale Conner stated "At least 24 subjects should be enrolled"
(http://www.fda.gov/downloads/AdvisoryCommittees/CommitteesMeetingMaterials/Drugs/
AdvisoryCommitteeforPharmaceuticalScienceandClinicalPharmacology/UCM178927.pdf
slide 18). The TRT | RTR design is less efficient in estimating Swr than
the 4-period full replicate and the partial replicate since Swr is
estimated from only half of the subjects (i.e. in the RTR sequence).
Recently the EMA stated in the Q&A document that 12 eligible subjects
are required in a TRT | RTR design. I have no clue what's behind the
minimum 24 for the FDA but if it is the 'reliability' of the Swr, in
some cases one has to increase the sample size. Try this code:
library(PowerTOST)
des <- c("2x2x4", "2x3x3", "2x2x3")
# type known.designs() for information
GMR <- 0.90
pwr <- 0.80
CV <- seq(0.3, 1, 0.02)
pch <- c(1, 2, 0)
cex <- c(1.5, 1.10, 1.25)
col <- c("black", "blue", "red")
l1 <- paste0("\nIntra-subject CV of R: ",
min(CV), " to ", max(CV), " (step size ",
unique(signif(diff(CV), 2)), ")",
"\nExpected GMR: ", GMR, ", target power: \u2265",
pwr, "\n")
l2 <- vector("character", length=length(des))
for(j in seq_along(des)) {
n <- vector("numeric", length=length(CV))
for(k in seq_along(CV)) {
n[k] <- sampleN.RSABE(CV=CV[k], theta0=GMR,
targetpower=pwr, design=des[j],
print=F, details=F)[["Sample size"]]
n.min <- min(n)
CVs <- CV[which(n == n.min, arr.ind=T, useNames=T)]
l2[j] <- paste0("\nDesign: ", des[j],
"\nMinimum sample size ", n.min, " at CVs: ",
min(CVs), " to ", max(CVs), "\n", collapse="")
}
if(j == 1) {
plot(100*CV, n, las=1, xlab="CV (%)", ylim=c(0, 60),
pch=pch[j], cex=cex[j], col=col[j],
main=paste0("RSABE (GMR ", GMR, ", target power \u2265 ", pwr, ")"))
abline(h=24, lty=3)
legend("topleft", title="design", legend=des, pch=pch, pt.cex=cex,
col=col)
} else {
points(100*CV, n, pch=pch[j], cex=cex[j], col=col[j])
}
}}
cat(l1, l2, "\n")
All the best,
Helmut
--
Ing. Helmut Schütz
BEBAC - Consultancy Services for
Bioequivalence and Bioavailability Studies
Neubaugasse 36/11
1070 Vienna, Austria
| PharmPK Discussion List Archive Index page |
Copyright 1995-2014 David W. A. Bourne (david@boomer.org)