Skip to content

linear mixed model required for the U.S. FDA

2 messages · Helmut Schütz, Thierry Onkelinx

Original
Helmut Schütz
# 
Dear all,

I?m struggling to set up a model required for the FDA (haha, and the 
Chinese agency). The closest I could get given at the end (which matches 
the one preferred by other regulatory agencies worldwide). The FDA is 
happy with R but "close" is not close /enough/.

Don't hit me. I'm well aware of the community's attitudes towards SAS. 
I'm not a SASian myself (software agnostic) but that's not related to 
SAS; one could set up this model in other (commercial...) software as well.

The FDA?s model allows different subject effects for each treatment 
(i.e., a subject-by-treatment interaction), and therefore, has 5 
variance terms:
 ? 1. within subject for T
 ? 2. within subject for R
 ? 3. between subject for T
 ? 4. between subject for R
 ? 5. covariance for between subject Test and Reference
The last three are combined to give the subject ? formulation 
interaction variance component.

The code provides a lot of significant digits only for comparison.

# FDA 2001 (APPENDIX E)
# https://www.fda.gov/media/70958/download
# FDA 2011 (p. 8)
# 
https://www.accessdata.fda.gov/drugsatfda_docs/psg/Progesterone_caps_19781_RC02-11.pdf
###############################################
# PROC MIXED;???????????????????????????????? #
# CLASSES SEQ SUBJ PER TRT;?????????????????? #
# MODEL Y = SEQ PER TRT/ DDFM = SATTERTH;???? #
# RANDOM TRT/TYPE = FA0(2) SUB = SUBJ G;????? #
# REPEATED/GRP=TRT SUB = SUBJ;??????????????? #
# ESTIMATE 'T vs. R' TRT 1 -1/CL ALPHA = 0.1; #
###############################################
# Example data set (EMA)
# 
https://www.ema.europa.eu/en/documents/other/31-annex-ii-statistical-analysis-bioequivalence-study-example-data-set_en.pdf
library(RCurl)
library(lme4)
library(emmeans)
url? <- getURL("https://bebac.at/downloads/ds01.csv")
data <- read.csv(text = url, colClasses=c(rep("factor", 4), "numeric"))
mod? <- lmer(log(PK) ~ period + sequence + treatment + (1|subject),
 ?????????????????????? data = data)
res1 <- test(pairs(emmeans(mod, ~ treatment, mode = "satterth"),
 ?????????????????? reverse = TRUE), delta = log(0.8))
res2 <- confint(emmeans(mod, pairwise ~ treatment, mode = "satterth"),
 ??????????????? level = 0.9)
# Workaround at the end because of lexical order
#?? I tried relevel(data$treatment, ref = "R") /before/ the model
#?? However, is not observed by confint(...)
cat(paste0("\nEMA Example data set 1",
 ?????????? "\nAnalysis of log-transformed data",
 ?????????? "\nSatterthwaite\u2019s degrees of freedom, 90% CI",
 ?????????? "\n\n? SAS 9.4, Phoenix/WinNonlin 8.1",
 ?????????? "\n?????????????????? mean???????? SE?????? df p.value",
 ?????????? "\n??? R????? :? 7.6704296 0.10396421? 74.762420",
 ?????????? "\n??? T????? :? 7.8158939 0.09860609? 74.926384",
 ?????????? "\n??? T vs. R:? 0.1454643 0.04650124 207.734958 0.00201129",
 ?????????? "\n???????????????????? PE? lower.CL? upper.CL",
 ?????????? "\n??? antilog:? 1.1565764 1.0710440 1.2489393",
 ?????????? "\n\n? lmer (lme 1.1-21), emmeans 1.4",
 ?????????? "\n?????????????????? mean???????? SE?????? df p.value",
 ?????????? "\n??? R????? :? ", sprintf("%.7f %.8f? %3.6f",
 ?????????????????????????????????????? res2$emmeans$emmean[1],
 ?????????????????????????????????????? res2$emmeans$SE[1],
 ?????????????????????????????????????? res2$emmeans$df[1]),
 ?????????? "\n??? T????? :? ", sprintf("%.7f %.8f? %3.6f",
 ?????????????????????????????????????? res2$emmeans$emmean[2],
 ?????????????????????????????????????? res2$emmeans$SE[2],
 ?????????????????????????????????????? res2$emmeans$df[2]),
 ?????????? "\n??? T vs. R:? ", sprintf("%.7f %.8f %3.6f %.8f",
 ?????????????????????????????????????? res1$estimate, res1$SE, res1$df,
 ?????????????????????????????????????? res1$p.value),
 ?????????? "\n???????????????????? PE? lower.CL? upper.CL",
 ?????????? "\n??? antilog:? ", sprintf("%.7f %.7f %.7f",
exp(-res2$contrasts$estimate),
exp(-res2$contrasts$upper.CL),
exp(-res2$contrasts$lower.CL)), "\n"))

Cheers,
Helmut
Thierry Onkelinx
# 
Dear Helmut,

The mixed models list is more suitable for this kind of question. I'm
forwarding it to that list. Please send any follow-up to that list instead
of the general R help list.

If I understand correctly, you'll need a different variance term for both
treatments (the within subject for T and R). I don't think you can do that
with lmer(). However, you can with nlme::lme() by using the weights
argument. The model does not converge on my machine.

library(nlme)
model2 <- lme(log(PK) ~ period + sequence + treatment , random = ~
treatment | subject, data = data, weights = varIdent(~treatment))

Another option is to go Bayesian with the INLA package (r-inla.org). Note
that the data needs some preparing. And the summary returns the precision
(1/var).

data$lPK_T <- ifelse(data$treatment == "T", log(data$PK), NA)
data$lPK_R <- ifelse(data$treatment == "R", log(data$PK), NA)
data$subject_T <- as.integer(data$subject)
n_subject <- max(data$subject_T)
data$subject_R <- ifelse(data$treatment == "R", data$subject_T + n_subject,
NA)
data$subject_T[data$treatment == "R"] <- NA

library(INLA)
model3 <- inla(
  cbind(lPK_T, lPK_R) ~ period + sequence + treatment +
    f(subject_T, model = "iid2d", n = 2 * n_subject) +
    f(subject_R, copy = "subject_T"),
  data = data,
  family = c("gaussian", "gaussian")
)
summary(model3)

Fixed effects:
                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
(Intercept)   7.6501 0.1529     7.3492   7.6501     7.9507  7.6501   0
period2       0.0423 0.0729    -0.1011   0.0423     0.1854  0.0423   0
period3       0.0057 0.0613    -0.1148   0.0057     0.1262  0.0057   0
period4       0.0718 0.0731    -0.0718   0.0718     0.2153  0.0718   0
sequenceTRTR -0.0218 0.1960    -0.4076  -0.0218     0.3636 -0.0217   0
treatmentT    0.1462 0.0597     0.0288   0.1462     0.2636  0.1462   0

Random effects:
Name  Model
 subject_T   IID2D model
subject_R   Copy

Model hyperparameters:
                                             mean     sd 0.025quant
0.5quant 0.975quant   mode
Precision for the Gaussian observations    9.4943 1.4716     6.8915
9.3972    12.6699 9.2192
Precision for the Gaussian observations[2] 5.7145 0.8390     4.2257
5.6602     7.5228 5.5594
Precision for subject_T (component 1)      1.4670 0.2541     1.0265
1.4471     2.0243 1.4092
Precision for subject_T (component 2)      1.3545 0.2436     0.9350
1.3345     1.8913 1.2962
Rho1:2 for subject_T                       0.9176 0.0236     0.8631
0.9205     0.9551 0.9261

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx at inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

<https://www.inbo.be>


Op ma 19 aug. 2019 om 12:29 schreef Helmut Sch?tz <helmut.schuetz at bebac.at>: