linear mixed model required for the U.S. FDA

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

linear mixed model required for the U.S. FDA

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

--
Ing. Helmut Schütz
BEBAC – Consultancy Services for
W https://bebac.at/
C https://bebac.at/Contact.htm
F https://forum.bebac.at/

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: linear mixed model required for the U.S. FDA

R help mailing list-2
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
[hidden email]
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 <[hidden email]>:

> 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
>
> --
> Ing. Helmut Schütz
> BEBAC – Consultancy Services for
> W https://bebac.at/
> C https://bebac.at/Contact.htm
> F https://forum.bebac.at/
>
> ______________________________________________
> [hidden email] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.