Bioequivalence and Bioavailability Forum • extra-reference design: R code

extra-reference design: R code [RSABE / ABEL]

posted by mittyri – Russia, 2025-12-23 12:58 (176 d 22:25 ago) – Posting: # 24537
Views: 1,896
(edited on 2025-12-23 15:33)

Dear BEQool!

❝ This seems unreasonable, because the additional reference data contain some additional information about reference treatment and should therefore probably influence the PE? :confused: So my question would be what is the correct statistical ANOVA model for an extra-reference design?

The model is correct. See explanations in the code below

library(PowerTOST)

library(emmeans)



set.seed(1234)



seqs  <- c("TRR", "RTR")

nseq  <- c(12, 12)

ldiff <- log(0.95)

s2wT  <- 0.2^2

s2wR  <- 0.1^2



d0 <- PowerTOST:::prep_data2(

  seqs = seqs,

  nseq = nseq,

  muR  = log(10),

  ldiff = ldiff,

  s2wT  = s2wT,

  s2wR  = s2wR

)



d0$subject  <- factor(d0$subject)

d0$sequence <- factor(d0$sequence)

d0$period   <- factor(d0$period)

d0$tmt      <- factor(d0$tmt, levels = c("R", "T"))



# Say you want period 3 R values about 50% lower

eff.per <- c("1" = 0,

             "2" = 0,

             "3" = -log(2)) 



d0$logCmax <- d0$logval + eff.per[d0$period]


sd.subj <- 0.3 # BSW



n.subj <- nlevels(d0$subject)

b.subj <- rnorm(n.subj, mean = 0, sd = sd.subj)



d0$logCmax <- d0$logCmax + b.subj[d0$subject]


d0$Cmax <- exp(d0$logCmax)



# in period 3, the model can only estimate the sum of Period3 and TreatmentR

m.full <- lm(logCmax ~ sequence + subject %in% sequence + period + tmt, data = d0)

# There is no information at all to distinguish Period and Reference

# (no T to distinguish)

# due to residual DF changes F for trt effects is also changed

anova(m.full)

# Analysis of Variance Table

# 

# Response: logCmax

# Df Sum Sq Mean Sq  F value    Pr(>F)    

# sequence          1 0.0010  0.0010   0.0345    0.8534    

# period            2 8.6272  4.3136 146.8938 < 2.2e-16 ***

#   tmt               1 0.0673  0.0673   2.2932    0.1369    

# sequence:subject 22 7.1322  0.3242  11.0398 1.124e-11 ***

#   Residuals        45 1.3215  0.0294                       

# ---

#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



emm.full <- emmeans(m.full, ~ tmt)

# The least-squares estimate of the treatment contrast (T–R) uses only periods 1 and 2

# Any change in period 3’s reference mean simply changes the period 3 effect, not the treatment effect

# Therefore the PE (T/R ratio) is exactly the same as in the pure 2-period 2×2 core:

PE.full  <- summary(contrast(emm.full, "revpairwise"), infer = TRUE)

PE.full

# contrast estimate     SE df lower.CL upper.CL t.ratio p.value

# T - R     -0.0749 0.0495 45   -0.175   0.0247  -1.514  0.1369

# 

# Results are averaged over the levels of: subject, period, sequence 

# Confidence level used: 0.95 



emm.period.full <- emmeans(m.full, ~ period)

summary(contrast(emm.period.full, "revpairwise"), infer = TRUE)

# contrast          estimate     SE df lower.CL upper.CL t.ratio p.value

# period2 - period1 -0.00261 0.0495 45   -0.123    0.117  -0.053  0.9985

# period3 - period1 -0.77306 0.0553 45   -0.907   -0.639 -13.978 <0.0001

# period3 - period2 -0.77045 0.0553 45   -0.904   -0.636 -13.930 <0.0001

# 

# Results are averaged over the levels of: subject, tmt, sequence 

# Confidence level used: 0.95 

# Conf-level adjustment: tukey method for comparing a family of 3 estimates 

# P value adjustment: tukey method for comparing a family of 3 estimates

# Periods comparison:

dtwo <- d0[d0$period != 3,]


m.two <- lm(logCmax ~ sequence + subject %in% sequence + period + tmt, data = dtwo)

anova(m.two)

# Analysis of Variance Table

# 

# Response: logCmax

# Df Sum Sq  Mean Sq F value    Pr(>F)    

# sequence          1 0.0040 0.004008  0.1146    0.7382    

# period            1 0.0001 0.000082  0.0023    0.9619    

# tmt               1 0.0673 0.067342  1.9249    0.1792    

# sequence:subject 22 5.1212 0.232782  6.6538 1.877e-05 ***

#   Residuals        22 0.7697 0.034985                      

# ---

#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

emm.two <- emmeans(m.two, ~ tmt)

PE.two  <- summary(contrast(emm.two, "revpairwise"), infer = TRUE)

PE.two

# contrast estimate    SE df lower.CL upper.CL t.ratio p.value

# T - R     -0.0749 0.054 22   -0.187   0.0371  -1.387  0.1792

# 

# Results are averaged over the levels of: subject, period, sequence 

# Confidence level used: 0.95 



# just for reference, period comparison:

emm.period.two <- emmeans(m.two, ~ period)

summary(contrast(emm.period.two, "revpairwise"), infer = TRUE)

# contrast          estimate    SE df lower.CL upper.CL t.ratio p.value

# period2 - period1 -0.00261 0.054 22   -0.115    0.109  -0.048  0.9619

# 

# Results are averaged over the levels of: subject, tmt, sequence 

# Confidence level used: 0.95 



#          *

#         ***

#        *****

#       *******

#      *********

#     ***********

#    *************

#   ***************

#  *****************

# *******************

#      *********

#     ***********

#    *************

#   ***************

#  *****************

# *******************

#    *************

#   ***************

#  *****************

# *******************

#         |||

#         |||

#         |||

#    *************

#   ***************

#  *****************

# *******************

#

#  Merry Christmas!

—
Kind regards,
Mittyri

Complete thread:

Partial replicate design Silva 2013-12-31 12:17 [RSABE / ABEL]
- extra-reference design Helmut 2013-12-31 14:28
  - extra-reference design BEQool 2025-12-22 18:27
    - extra-reference design: R codemittyri 2025-12-23 11:58
      - extra-reference design: R code BEQool 2025-12-29 14:17
        
        extra-reference design: no period effect assumption mittyri 2025-12-29 22:16
        
        extra-reference design: no period effect assumption BEQool 2025-12-31 16:42