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.
>
[[alternative HTML version deleted]]