linear mixed model required for the U.S. FDA
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>:
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/
______________________________________________ R-help at r-project.org 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.