- On 4 Aug 2015 at 16:24:34, Dennis Fisher (fisher.at.plessthan.com) sent the message

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 - On 4 Aug 2015 at 22:38:02, Leonid Gibiansky (lgibiansky.at.quantpharm.com) sent the message

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 - On 5 Aug 2015 at 12:54:11, Boucher, Joseph F. (joseph.f.boucher.-at-.zoetis.com) sent the message

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 - On 5 Aug 2015 at 12:57:21, Leonid Gibiansky (lgibiansky.-at-.quantpharm.com) sent the message

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 - On 5 Aug 2015 at 12:59:04, Dennis Fisher (fisher.at.plessthan.com) sent the message

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) - On 5 Aug 2015 at 13:21:12, Dennis Fisher (fisher.aaa.plessthan.com) sent the message

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) - On 6 Aug 2015 at 14:38:49, Helmut Schutz (helmut.schuetz.at.bebac.at) sent the message

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)