Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-04-20 19:34
(3624 d 01:45 ago)

Posting: # 14718
Views: 26,829
 

 Bear vs. Phoenix & SAS [🇷 for BE/BA]

Post reply
Dear all,

I closed this thread. Let’s continue over here. I would say that nowadays few people perform a replicated study without having reference-scaling in mind. What would we need in bear?
  • A basic model which we use if the conditions for scaling are not fulfilled.
    For the FDA that would be a mixed-effects model as given in the 2001 guidance and the pro­ge­sterone guidance.
    For the EMA we would need an “all effects fixed” model (Method A of the Q&A).
  • FDA – RSABE: The progesterone guidance, again.
    EMA – ABEL: The Q&A, Methods A & B. Crippled model to get CVwR.
We have shown in the past that we get the same results for EMA’s data sets I & II (also when we make #II unbalanced) in PHX and SAS (THX to Shuanghe and Jean-Michel!). I would say, that’s the target.

PS: @ElMaestro. Remember that one of the referees of our reference dataset-MS wanted to discuss the statistical model and we refused? Until somebody shows that what regulators want right now is crap (which will happen – at least partly) we should get the same results independent from the software used. I would be happy if a noncommercial one is amongst them.
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2015-04-20 21:36
(3623 d 23:43 ago)

@ Helmut
Posting: # 14719
Views: 24,229
 

 R vs. Phoenix & SAS?

Post reply
Dear Helmut,

❝ What would we need in bear?


Wow! thank you so much. But that's too much. It must be a typo. May I suggest as "what we can do with R?" There are lots of R users/gurus in this Forum.

❝ For the FDA that would be a mixed-effects model as given in the 2001 guidance and the pro­ge­sterone guidance.


I like this.

❝ For the EMA we would need an “all effects fixed” model (Method A of the Q&A).


One stupid question :-D: if EMA needs an "all effects fixed" model, why do we use a mixed model to do analysis? I don't remember lme() right now, but lmer() must have at least a random variable. I tried lmer() last night and was pretty sure about that. Does EMA imply that we should use lm() to analyze replicate crossover dataset?:confused:

❝ • FDA – RSABE: Again the progesterone guidance.

❝ EMA – ABEL: The Q&A, Methods A & B. Crippled model to get CVwR...


❝ I would say, that’s the target.


You should try to get funded from EMA or FDA for these projects. These are really great projects, indeed.
—
All the best,
-- Yung-jin Lee
bear v2.9.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-04-21 03:02
(3623 d 18:18 ago)

@ yjlee168
Posting: # 14725
Views: 24,133
 

 R vs. Phoenix & SAS?

Post reply
Hi Yung-jin,

❝ ❝ For the EMA we would need an “all effects fixed” model (Method A of the Q&A).

❝ if EMA needs an "all effects fixed" model, why do we use a mixed model to do analysis?


Both the BE-GL and the Q&A state that all effect should be fixed. EMA generally accepts a mixed model, but may ask for recalculation in borderline cases (e.g., extreme imbalance, CI close to the AR). Sometimes assessors ask for recalculation even if the CI was very narrow and the PE was close to 1. Happened to me once. The Q&A specifically mentions SAS’ Proc GLM, not Proc Mixed.

❝ I don't remember lme() right now,…


In lme() it should be possible to specify all effects fixed. Actually the random effect(s) are optional. If is doesn’t work, lm() should do.

❝ Does EMA imply that we should use lm() to analyze replicate crossover dataset?:confused:


No. Both the FDA and the EMA do not mandate any software or routine within. Quote from the Q&A:

SAS (version 9.1, SAS Institute Inc, NC) was used in the previous computations. Results obtained by alternative, validated statistical programs are also acceptable except spread­sheets because outputs of spread­sheets are not suitable for secondary assessment.


❝ You should try to get funded from EMA or FDA for these projects. These are really great projects, indeed.


Gimme a break. I’m a lousy R-coder. Agencies want to get results from validated software. They don’t pay for its development.
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2015-04-22 01:41
(3622 d 19:39 ago)

@ Helmut
Posting: # 14730
Views: 25,012
 

 lme() does not work with all fixed effects

Post reply
Dear Helmut,

❝ ... The Q&A specifically mentions SAS’ Proc GLM, not Proc Mixed.


glm() works in R.

❝ In lme() it should be possible to specify all effects fixed. Actually the random effect(s) are optional. If is doesn’t work, lm() should do.


lme() does not work without random effect. Am I doing anything wrong?

> require(nlme)
Loading required package: nlme
> cnames<-c("subj","drug","seq", "prd","Cmax", "AUC0t", "AUC0INF","lnCmax","lnAUC0t","lnAUC0INF")
> TotalData<-read.csv(file="SingleRep_stat_demo.csv",header=TRUE,row.names=NULL,col.names=cnames, sep=",",dec=".")
> mod.lme<-lme(log(Cmax) ~ subj + seq +  prd + drug, data=TotalData, method="REML")
Error in getGroups.data.frame(dataMix, groups) :
  invalid formula for groups

If I compared lme(), lm(), and glm(), I got
> mod.lme<-lme(log(Cmax) ~ seq +  prd + drug, random=~drug -1|subj, data=TotalData, method="REML")  ###lme in bear
> mod.lm<-lm(log(Cmax) ~ subj + seq +  prd + drug, data=TotalData)    ## lm()
> mod.glm<-glm(log(Cmax) ~ subj + seq +  prd + drug, data=TotalData)  ## glm()
> AIC(mod.lme, mod.lm)
        df       AIC
mod.lme  6 -49.88309
mod.lm   6 -70.88094
Warning message:
In AIC.default(mod.lme, mod.lm) :
  models are not all fitted to the same number of observations
> AIC(mod.lme, mod.glm)
        df       AIC
mod.lme  6 -49.88309
mod.lm   6 -70.88094
Warning message:
In AIC.default(mod.lme, mod.lm) :
  models are not all fitted to the same number of observations

Summary of mod.lme
Linear mixed-effects model fit by REML
 Data: TotalData
        AIC       BIC   logLik
  -49.88309 -38.17563 30.94155

Random effects:
 Formula: ~drug - 1 | subj
              drug  Residual
StdDev: 0.04344759 0.1041244

Fixed effects: log(Cmax) ~ seq + prd + drug
                Value  Std.Error DF  t-value p-value
(Intercept)  7.314670 0.08349877 40 87.60213  0.0000
seq         -0.025701 0.04252039 12 -0.60445  0.5568
prd          0.016233 0.01244525 40  1.30438  0.1996
drug         0.046554 0.03015388 40  1.54387  0.1305
 Correlation:
     (Intr) seq    prd   
seq  -0.764             
prd  -0.373  0.000       
drug -0.461  0.000  0.000

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-2.8081617 -0.6114645  0.1380420  0.5476292  1.8620918

Number of Observations: 56
Number of Groups: 14

Summary of mod.lm
Call:
lm(formula = log(Cmax) ~ subj + seq + prd + drug, data = TotalData)

Residuals:
     Min       1Q   Median       3Q      Max
-0.29483 -0.07058  0.01211  0.07004  0.23080

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)  7.353905   0.087149  84.383   <2e-16 ***
subj        -0.007061   0.004041  -1.747   0.0866 . 
seq         -0.016555   0.032583  -0.508   0.6136   
prd          0.016233   0.014459   1.123   0.2668   
drug         0.046554   0.032332   1.440   0.1560   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.121 on 51 degrees of freedom
Multiple R-squared:  0.1126,    Adjusted R-squared:  0.04302
F-statistic: 1.618 on 4 and 51 DF,  p-value: 0.184

Summary of glm
Call:
glm(formula = log(Cmax) ~ subj + seq + prd + drug, data = TotalData)

Deviance Residuals:
     Min        1Q    Median        3Q       Max 
-0.29483  -0.07058   0.01211   0.07004   0.23080 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)  7.353905   0.087149  84.383   <2e-16 ***
subj        -0.007061   0.004041  -1.747   0.0866 . 
seq         -0.016555   0.032583  -0.508   0.6136   
prd          0.016233   0.014459   1.123   0.2668   
drug         0.046554   0.032332   1.440   0.1560   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.01463458)

    Null deviance: 0.84108  on 55  degrees of freedom
Residual deviance: 0.74636  on 51  degrees of freedom
AIC: -70.881

Number of Fisher Scoring iterations: 2

> AIC(mod.lm, mod.glm)
        df       AIC
mod.lm   6 -70.88094
mod.glm  6 -70.88094
—
All the best,
-- Yung-jin Lee
bear v2.9.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
 
Astea
★★  

Russia,
2016-11-04 01:13
(3060 d 19:06 ago)

@ yjlee168
Posting: # 16771
Views: 20,858
 

 lme() does not work with all fixed effects

Post reply
Dear all!

Trying to understand lm and lme I found a simple and very clear instruction (Lme tutorial). Hope it would be useful for someone not familiar with such models.

I try to evaluate Data set I from QA in R. The model is as follows:

> mod.lme<-lme(log(Cmax) ~ seq +  prd + drug, random=~drug|subj/seq, data=QA)
summary(mod.lme)

Linear mixed-effects model fit by REML
 Data: QA
     AIC    BIC  logLik
  548.86 589.38 -263.43

Random effects:
 Formula: ~drug | subj
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr 
(Intercept) 0.750263 (Intr)
drug        0.044904 -0.885

 Formula: ~drug | seq %in% subj
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr 
(Intercept) 0.530707 (Intr)
drug        0.037234 -0.844
Residual    0.393976       

Fixed effects: log(Cmax) ~ seq + prd + drug
              Value Std.Error  DF t-value p-value
(Intercept)  7.4432   0.32947 219 22.5912  0.0000
seq         -0.0248   0.19682  75 -0.1261  0.9000
prd          0.0468   0.02052 219  2.2811  0.0235
drug         0.1465   0.04627 219 3.1665  0.0018
 Correlation:
     (Intr) seq    prd   
seq  -0.909             
prd  -0.204  0.053       
drug -0.251  0.006  0.000

Standardized Within-Group Residuals:
      Min        Q1       Med        Q3       Max
-3.001657 -0.408221 -0.023744  0.340791  5.038066

Number of Observations: 298
Number of Groups:
         subj seq %in% subj
           77            77


Then I try to interpret the data in order to estimate CI. I took 0.1465 as point estimation (exp(0,1465)=1,1578 and calculate CI as PE+- SE*t(0,1;219), where SE=0.04627.
The result is 115,78: 107,26 - 124,97 (comparing with method B from QA might be 115,73: 107,17- 124,97). Not ideal but seems to be similar. What is the reason for difference: rounding or different model, or am I do something wrong? :confused:

Later I try to recalculate it using bear and get

> Fixed effects: log(Cmax) ~ seq + prd + drug
   Value Std.Error  DF t-value p-value
  7.6332  0.151688 217  50.321  0.0000
 -0.0196  0.197671  75  -0.099  0.9212
  0.0003  0.064043 217   0.004  0.9965
  0.0381  0.062413 217   0.610  0.5425
  0.1474  0.063881 217   2.308  0.0219
  0.1449  0.046986 217   3.083  0.0023

 Point Estimate   CI90 lower   CI90 upper
          115.587      106.955      124.916


And finally if I try to use lme just from yjlee168's (post#14730),
< mod.lme<-lme(log(Cmax) ~ seq +  prd + drug, random=~drug -1|subj, data=TotalData, method="REML")

I got another result:
Fixed effects: log(Cmax) ~ seq + prd + drug
              Value Std.Error  DF t-value p-value
(Intercept)  7.4154  0.225498 219  32.885  0.0000
seq         -0.0007  0.141455  75  -0.005  0.9959
prd          0.0467  0.028185 219   1.658  0.0988
drug          0.1415 0.081819 219    1.730  0.0850

However it seemed to me strange, what model is in bear now?

And I still wonder how does the model affect the Point Estimation? For QA data set 1 I've got manually Tmean=7,830303 and Rmean=7,676036, T-R=0,1543 far from the estimated via lme PE.

Trying to understand afermentioned Phoenix results I see non-integer df. How can they be achieved manually?
—
"Being in minority, even a minority of one, did not make you mad"
 
mittyri
★★  

Russia,
2016-11-05 18:38
(3059 d 01:41 ago)

@ Astea
Posting: # 16777
Views: 20,729
 

 lmer: Method B (PE catched for imbalanced dataset!!!) and Method C

Post reply
Dear Astea!

so many questions ;-)
I'll just put here my code as a solution. Yes, we cannot get SAS results, sorry...
Please let me know if after investigation you still have some questions

library(readxl)
library(lmerTest)
options(contrasts=c("contr.treatment","contr.poly"))
methodC <- function(Dataset){
  muddle.lmer <- (lmer(log(Data)~Sequence+Period+Formulation + (Formulation-1|Subject), data=Dataset))
  cat("\n", "lmer treats Sequence, Period and Formulation as factors, Random is /~Formulation -1|Subject/ (close to Method C), then", "\n")
  FormulationEffect.lmer <- summary(muddle.lmer)$coefficients["FormulationT","Estimate"]   print(FormulationEffect.lmer)
  cat("\n", "PE is", round(exp(as.numeric(FormulationEffect.lmer))*100, digits = 2),
      "\n", "SAS gives 115.66", "\n")
 
  cat("\n", "what about LSMeans?", "\n", "PHX gives for Method C  7.67042954073919 and 7.81589387050252", "\n", "and lmer gives", "\n")
  LSM.lmer <- lmerTest::lsmeans(muddle.lmer, "Formulation")
  print(LSM.lmer$lsmeans.table[c(1,2), c(3,4)])
 
  cat("\n", "\n", "what about CI? SAS gives 107.10-124.89", "\n")
  round(exp(confint(muddle.lmer, "FormulationT", level = 0.90))*100, 2)
}

methodB <- function(Dataset){
  muddle.lmer <- (lmer(log(Data)~Sequence+Period+Formulation + (1|Subject), data=Dataset))
  cat("\n", "lmer treats Sequence, Period and Formulation as factors, Random is /~1|Subject/ (close to Method B), then", "\n")
  FormulationEffect.lmer <- summary(muddle.lmer)$coefficients["FormulationT","Estimate"]   print(FormulationEffect.lmer)
  cat("\n", "PE is", round(exp(as.numeric(FormulationEffect.lmer))*100, digits = 2),
      "\n", "SAS gives 115.73", "\n")
 
  cat("\n", "what about LSMeans?", "\n", "PHX gives for Method B 7.67001367911898 and 7.81610190985527",  "\n", "and lmer gives", "\n")
  LSM.lmer <- lmerTest::lsmeans(muddle.lmer)
  print(LSM.lmer$lsmeans.table[c(7,8), c(3,4)])
 
  cat( "\n", "\n", "what about CI? SAS gives 107.17-124.97", "\n")
  round(exp(confint(muddle.lmer, "FormulationT", level = 0.90))*100, 2)
}

# downloading the reference file from BEBAC; that's the most convenient way to get the file
download.file('http://bebac.at/downloads/Validation Replicate Design EMA.xls', "Dataset.xls", cacheOK = FALSE, mode="wb")
# read it without Perl and Java; thanks, Hadley!
Dataset<-read_excel("Dataset.xls", sheet = 1)
# prepare dataset

Dataset$Formulation<-factor(Dataset$Formulation, levels = c("R", "T")) # for LSMeans in lmer
Dataset$Sequence<-factor(Dataset$Sequence, levels = c("TRTR", "RTRT")) # for LSMeans in lmer
Dataset$Period<-factor(Dataset$Period) # for LSMeans in lmer

methodC(Dataset)
methodB(Dataset)

—
Kind regards,
Mittyri
 
Astea
★★  

Russia,
2016-11-05 20:27
(3058 d 23:53 ago)

@ mittyri
Posting: # 16778
Views: 20,600
 

 lmer: Method B (PE catched for imbalanced dataset!!!) and Method C

Post reply
mittyri!!!
May I hug you?! This is what I actually wanted to get spending two nights in a desperate attempt to apply Alexandra Kuznetsova's lmerTest to BE!

And the questions.. Still there are too many questions.

Why didn't you provide the result of the code (that's the way of secret initiation?)

First start with method B. My PE coincides fully with QA result (115.73), the CI calculated by your method slightly differs (107.26-124.88). But! What is that function - confint? May be something wrong with it? Cause if I use PE+-SE*t(0.1;219) where SE is 4.65e-02 I get 1,0717-1,2797 totally the target!!!

Please explain how to get Std.Error and Value with more digits?
I tried to add "digits=8" but nothing happened...

Why does the model have (1|Subject) as random? I also tried (1|Subject/Sequence) and the result (PE and CI) remains the same. Can we use it instead?

The other question is what from the output we can use to estimate R-variance in order to apply scABE?

We are in front of the final step - from now we can make R tutorial to hold scABEL as did Helmut in WinNonLin, am I right? If it is so, Es ist fantastisch!:ok:
—
"Being in minority, even a minority of one, did not make you mad"
 
mittyri
★★  

Russia,
2016-11-05 21:01
(3058 d 23:18 ago)

@ Astea
Posting: # 16779
Views: 20,615
 

 lmer: Method B is ready for scaling

Post reply

❝ mittyri!!!

❝ May I hug you?!

:flower::party:

❝ Why didn't you provide the result of the code (that's the way of secret initiation?)


I was so surprized by the results that was in hurry to share my code :-D

❝ Please explain how to get Std.Error and Value with more digits?

Oh, you are right, please substitute confint string to
  FormulationT<-summary(muddle.lmer)$coefficients["FormulationT",]   cat(round(exp(FormulationT[1]+FormulationT[2]*qt(0.05, FormulationT[3]))*100, 2), "-",
      round(exp(FormulationT[1]-FormulationT[2]*qt(0.05, FormulationT[3]))*100, 2))
confint is not the right thing here. Here we go:

lmer treats Sequence, Period and Formulation as factors, Random is /~1|Subject/ (close to Method B), then
[1] 0.1460882
PE is 115.73
SAS gives 115.73
what about LSMeans?
PHX gives for Method B 7.67001367911898 and 7.81610190985527
and lmer gives
               Formulation Estimate
Formulation  R           R   7.6700
Formulation  T           T   7.8161
what about CI? SAS gives 107.17-124.97
107.17 - 124.97

❝ Why does the model have (1|Subject) as random? I also tried (1|Subject/Sequence) and the result (PE and CI) remains the same. Can we use it instead?


Yes, we can. Subject in Sequence is SAS legacy. By the way as ElMaestro wrote, that's meaningless, because the reference grid (G-matrix) doesn't change. It will be useful only if we are enumerating the subjects in each sequence from 1 (Seq ABAB: 1, 2, 3...; Seq BABA: 1, 2, 3...)

❝ The other question is what from the output we can use to estimate R-variance in order to apply scABE?


No, we cannot. Sorry, I must confess I don't have time now to implement this, but I would do it as follows:
- exclude all T data first and then
- build a linear model with lm() using only reference data.

If nobody does it during the next week I'll try to code it by myself.

❝ We are in front of the final step - from now we can make R tutorial to hold scABEL as did Helmut in WinNonLin, am I right? If it is so, Es ist fantastisch!:ok:


Ja, new feature to implement in BEAR, Yung-jin will be happy :cool:
—
Kind regards,
Mittyri
 
Astea
★★  

Russia,
2016-11-06 12:50
(3058 d 07:30 ago)

@ mittyri
Posting: # 16781
Views: 20,659
 

 lmer: Method B is ready for scaling

Post reply

❝ If nobody does it during the next week I'll try to code it by myself.


Challenge accepted!

Here is my rough code:

library(readxl)
library(lmerTest)
options(contrasts=c("contr.treatment","contr.poly"))
Dataset<-read_excel("Dataset.xls", sheet = 1)

methodB.scABEL <- function(Dataset){
Dataset$Formulation<-factor(Dataset$Formulation, levels = c("R", "T"))
Dataset$Sequence<-factor(Dataset$Sequence, levels = c("TRTR", "RTRT"))
Dataset$Period<-factor(Dataset$Period)
mod.lmer <- (lmer(log(Data)~Sequence+Period+Formulation + (1|Subject/Sequence), data=Dataset))
FormulationEffect.lmer <- summary(mod.lmer)$coefficients["FormulationT","Estimate"]   cat("\n", "PE is", round(exp(as.numeric(FormulationEffect.lmer))*100, digits = 2), "\n")
FormulationT<-summary(mod.lmer)$coefficients["FormulationT",]  
CI_Lr<-round(exp(FormulationT[1]+FormulationT[2]*qt(0.05, FormulationT[3]))*100, 2)
CI_Ur<-round(exp(FormulationT[1]-FormulationT[2]*qt(0.05, FormulationT[3]))*100, 2)
cat("\n", "CI is", CI_Lr, "-", CI_Ur,"\n")

Dataset2<-Dataset[order(Dataset$Formulation), ]
n<-length(which(Dataset2$Formulation=="R"))

Dataset3<-Dataset2[-c((n+1):nrow(Dataset2)),]
for (i in 1:nrow(Dataset3)) {if(Dataset3$Period[i]==3) Dataset3$Formulation[i]<-"T" else if(Dataset3$Period[i]==2) Dataset3$Formulation[i]<-"T"}
for (i in 1:nrow(Dataset3)) {if(Dataset3$Period[i]==3) Dataset3$Period[i]<-2 else if(Dataset3$Period[i]==4) Dataset3$Period[i]<-2 else if(Dataset3$Period[i]==2) Dataset3$Period[i]<-1}

Dataset3$Formulation<-factor(Dataset3$Formulation, levels = c("R", "T"))
Dataset3$Sequence<-factor(Dataset3$Sequence, levels = c("RTRT", "TRTR"))
Dataset3$Period<-factor(Dataset3$Period)
Dataset3$Subject<-factor(Dataset3$Subject)
muddle   <- lm(log(Data)~Formulation+Period+Sequence+Subject, data=Dataset3)
CV_R<-(exp(anova(muddle)[5,3])-1)^(1/2)*100
cat("\n", "CV_R is", round(CV_R, 2),"\n")
{if (CV_R>30)
 CI_L<-exp(-anova(muddle)[5,3]^(1/2)*0.760)*100
 CI_U<-exp(anova(muddle)[5,3]^(1/2)*0.760)*100
cat("Wide CI:", round(CI_L,2),"-", round(CI_U,2))}
{if ((CI_Lr<CI_L)|(CI_Ur>CI_U)) cat("BE failed")
else cat("=> BE success")}
}

methodB.scABEL(Dataset)

gives
PE is 115.73

 CI is 107.17 - 124.97

 CV_R is 46.96
Wide CI: 71.23 - 140.4=> BE success>
 
as in Phoenix (see Helmut's RSABE tutorial). The code needs to be tested via another data...
—
"Being in minority, even a minority of one, did not make you mad"
 
mittyri
★★  

Russia,
2016-11-07 07:07
(3057 d 13:13 ago)

@ Astea
Posting: # 16782
Views: 20,785
 

 lmer: Method B is ready for scaling

Post reply
Hi Astea,

❝ Here is my rough code:


Nice work! :ok:
Here's mine. I tried to prepare the dataframe for validation (Results)

library(readxl)
library(lmerTest)
library(dplyr)
options(contrasts=c("contr.treatment", "contr.poly"))
methodB <- function(Dataset){
  muddle.lmer <- (lmer(log(Data)~Sequence+Period+Formulation + (1|Subject), data=Dataset))
  cat("\n", "lmer treats Subject, Sequence, Period and Formulation as factors, Random is /~1|Subject/ (Method B), then")
  FormulationT.lmer <- summary(muddle.lmer)$coefficients["FormulationT",]
  Results<-data.frame(exp(as.numeric(FormulationT.lmer[1])),
    exp(FormulationT.lmer["Estimate"] + FormulationT.lmer["Std. Error"]*qt(0.05, FormulationT.lmer["df"])),
    exp(FormulationT.lmer["Estimate"] - FormulationT.lmer["Std. Error"]*qt(0.05, FormulationT.lmer["df"])))
  colnames(Results) <- c("PE", "Lower_CI", "Upper_CI")
  cat("\n", "PE is", sprintf("%.2f", Results$PE*100), "%\n")
  cat("CI is", sprintf("%.2f", Results$Lower_CI*100), "-", sprintf("%.2f", Results$Upper_CI*100), "%\n")
  # get reference dataset ready
  Ref<-Dataset %>%
    filter(Formulation=="R") %>%
    group_by(Subject) %>%
    filter(max(n())==2)
  Ref.lm <- lm(log(Data)~Sequence+Period+Subject, data=Ref)
  MSerror <-sum(Ref.lm$residuals^2)/Ref.lm$df.residual
  # get Stud Residuals ready
  Ref$StudResiduals <- rstandard(Ref.lm)
  StudResidualsFirst <- Ref %>%
    filter(Period=="1"|Period=="2")
  # plotting
  bxpdat <- boxplot(StudResidualsFirst$StudResiduals, range = 3, main ="Studentized Residuals Boxplot, IQR = 3")
  text(bxpdat$group,                                              # the x locations
       bxpdat$out,                                                # the y values
       StudResidualsFirst[match(StudResidualsFirst$StudResiduals, bxpdat$out, nomatch=0)>0, "Subject"]$Subject,  # the labels
       pos = 4) 
  qqnorm(StudResiduals)
  # estimating CVr
  Results$CV_R <- (exp(MSerror)-1)^(1/2)*100
  cat("CV_R is ", sprintf("%.2f",Results$CV_R), "%\n")
  # decision tree for CV_Rs
  if(Results$CV_R>30)  {
    if(Results$CV_R<=50)    {
      cat("certainly we need scaling\n")
      Results$Ext_L_Bound<-exp(-MSerror^(1/2)*0.760)
      Results$Ext_U_Bound<-exp(+MSerror^(1/2)*0.760) 
      cat("the scaled bounds are", sprintf("%.2f", Results$Ext_L_Bound*100), "-", sprintf("%.2f", Results$Ext_U_Bound*100), "%\n" )
    } else  {
      cat("certainly we need scaling, but the bounds should be fixed\n")
      Results$Ext_L_Bound<-0.6984
      Results$Ext_U_Bound<-1.4319
      cat("the scaled bounds are", sprintf("%.2f", Results$Ext_L_Bound*100), "-", sprintf("%.2f", Results$Ext_U_Bound*100), "%\n" )
    }
  }  else  {
    cat("we do not need scaling (bounds are 80.00-125.00 %)\n")
    Results$Ext_L_Bound<-80.00
    Results$Ext_U_Bound<-125.00
  }
  # print results
  if(Results$Lower_CI>Results$Ext_L_Bound & Results$Upper_CI<Results$Ext_U_Bound) {
    cat("BE is met, congrats")
  } else {
    cat("BE is not met, sorry")
  }
  return(Results)
}

# downloading the reference file from BEBAC; that's the most convenient way to get the file
download.file('http://bebac.at/downloads/Validation Replicate Design EMA.xls', "Dataset.xls", cacheOK = FALSE, mode="wb")
# read it without Perl and Java; thanks, Hadley!
Dataset<-read_excel("Dataset.xls", sheet = 1)
# prepare dataset
Dataset$Formulation<-factor(Dataset$Formulation, levels = c("R", "T"))
Dataset$Sequence<-factor(Dataset$Sequence, levels = c("TRTR", "RTRT"))
Dataset$Period<-factor(Dataset$Period)
Dataset$Subject<-factor(Dataset$Subject)
Results<- methodB(Dataset)

❝ as in Phoenix (see Helmut's RSABE tutorial). The code needs to be tested via another data...


I suspect deviations from SAS, because some values are not equal in full precision...
OK, that would be my next task :-D
—
Kind regards,
Mittyri
 
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2015-04-20 23:34
(3623 d 21:46 ago)

@ Helmut
Posting: # 14720
Views: 24,218
 

 info for lsmeans

Post reply
Dear Helmut,

❝ I closed this thread.


Since you have closed that thread, I like to thank you for this post first (very nice codes). But do you think that we should also present df, upper/lower CI (90% or 95%), and p values stuffs for lsmean? Or just keep it simple as you suggest? Thanks in advanced.
—
All the best,
-- Yung-jin Lee
bear v2.9.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-04-21 03:15
(3623 d 18:04 ago)

@ yjlee168
Posting: # 14726
Views: 24,164
 

 info for lsmeans

Post reply
Hi Yung-jin,

❝ I like to thank you for this post first


Welcome.

❝ (very nice codes).


I would never (!) use such a code given the nice vectorization possible in R. It was only intended as an example to get the means without a model.

❝ But do you think that we should also present df, upper/lower CI (90% or 95%), and p values stuffs for lsmean?


IMHO, the df makes sense since it tells us whether we pulled the right stuff from the model.
The CI is nice to know but not necessary. If the CIs overlap there is no significant difference between treatments. So what? No agency asks for it.
I never understood the purpose of the t- and p-values. What is tested here? A difference to zero or what?

❝ Or just keep it simple as you suggest?


Please don’t think about using this awful code!
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
Astea
★★  

Russia,
2016-11-03 00:43
(3061 d 19:37 ago)

@ Helmut
Posting: # 16767
Views: 20,719
 

 once more about R and replicate designes

Post reply
Dear all!

While reading tonns of posts I stil can't understand how do SAS has got such a unique position? The method C for replicate studies (that is random effect with interactions PROC MIXED) in R (not only in bear but also in other packages) doesn't exist, does it?

As I understand only method A is available in bear and only for partial replicate designs 2x2x4, 2x2x6.. But what about fully replicated 2x4x4? VStus (in that post) mentioned that

❝ bear's lm.mod() was not confused by having more than 2 periods and 2 sequences


Also a question appears what should we use as a reference variance in order to apply scaling of confidence intervals (fully replicated design allows us even to estimate intrasubject variability of reference drug for two pairs of R that is we can desintegrate 2x4x4 for two 2x2x4 and estimate R variance independently for both. How will choosing one of them affect the scaling and possibility to fail?)

And the thing that is totally out of my mind: why should point estimation be affected by the way of analysis? As it was once shown by Helmut for balanced studies PE is just a means over periods and sequences. Why can't we do analogous in the case of inbalanced designs?
—
"Being in minority, even a minority of one, did not make you mad"
 
VStus
★    

Poland,
2016-11-06 12:34
(3058 d 07:45 ago)

@ Astea
Posting: # 16780
Views: 20,481
 

 once more about R and replicate designes

Post reply
Dear Astea,

❝ As I understand only method A is available in bear and only for partial replicate designs 2x2x4, 2x2x6.. But what about fully replicated 2x4x4? VStus (in that post) mentioned that


I meant 3- and 4-arm non-replicate studies in my post, exploratory pilot comparing 2 or 3 pilot formulations versus the Reference. lm() or PROC GLM is enough in this case.

Regards, VStus
 
StatR
☆    

Turkey,
2017-02-03 14:53
(2969 d 05:27 ago)

@ VStus
Posting: # 17007
Views: 19,626
 

 Getting variance components from the lmer output

Post reply
Hi dear all,

First of all, I am sorry if my question is too stupid.

I could not find the answer no matter how much I investigate. I am new in the area of bioequivalence, and recently I am interested in individual bioequivalence. Also, I am not so familiar with the lme4. I am trying to extract the variance components (sigmaBT^2, sigmaBR^2, sigmaWT^2, sigmaWR^2 and rho) of the model (Restricted Maximum Likelihood Estimation). My data is as follows:
head(mydata)
  Subject Formula Period Sequence      Yijk
1       1       1      1        1 100.11145
2       1       0      2        1 100.00666
3       1       1      3        1 100.06645
4       1       0      4        1  99.96603
5       2       1      1        1 100.04678
6       2       0      2        1 100.03281

I am trying to do it by lmer function as follows:

lmer(Yijk~Formula+Period+Sequence+(1|Sequence))
Linear mixed model fit by REML ['lmerMod']
Formula: Yijk ~ Formula + Period + Sequence + (1 | Sequence)
REML criterion at convergence: 60.4953
Random effects:
 Groups   Name        Std.Dev.
 Sequence (Intercept) 0.002309
 Residual             0.356224
Number of obs: 64, groups:  Sequence, 2
Fixed Effects:
(Intercept)      Formula       Period     Sequence 
  99.787798     0.178645     0.080047    -0.004724 
model <- lmer(Yijk~Formula+Period+Sequence+(1|Sequence))
summary(model)
Linear mixed model fit by REML ['lmerMod']
Formula: Yijk ~ Formula + Period + Sequence + (1 | Sequence)

REML criterion at convergence: 60.5

Scaled residuals:
    Min      1Q  Median      3Q     Max
-2.5750 -0.4585  0.1577  0.5746  2.2238

Random effects:
 Groups   Name        Variance  Std.Dev.
 Sequence (Intercept) 5.331e-06 0.002309
 Residual             1.269e-01 0.356224
Number of obs: 64, groups:  Sequence, 2

Fixed effects:
             Estimate Std. Error t value
(Intercept) 99.787798   0.178187   560.0
Formula      0.178645   0.089056     2.0
Period       0.080047   0.039827     2.0
Sequence    -0.004724   0.089116    -0.1

Correlation of Fixed Effects:
         (Intr) Formul Period
Formula  -0.250             
Period   -0.559  0.000       
Sequence -0.750  0.000  0.000

I know, the estimated value of sigmaWR^2 = 1.269e-01 and the mean difference Ft-Fr = 0.178645.

But how can I get the estimated values of sigmaBT^2, sigmaBR^2, sigmaWT^2 and rho ? I appreciate any help.
 
VStus
★    

Poland,
2017-02-03 16:47
(2969 d 03:33 ago)

@ StatR
Posting: # 17009
Views: 19,214
 

 Getting variance components from the lmer output

Post reply
Dear StatR,

Sorry for confusion, deleted my incorrect reponse. summary(lmerdata) seems to be already provided.

What is Yijk?

Yijk seems to be a PK paramenter of interest... Residual you've extracted seems to be pooled both formulations.

Regards, VStus

PS: Helmut, thank you for options(digits=12)!
 
StatR
☆    

Turkey,
2017-02-03 18:12
(2969 d 02:08 ago)

@ VStus
Posting: # 17011
Views: 19,135
 

 Getting variance components from the lmer output

Post reply

❝ Yijk seems to be a PK paramenter of interest... Residual you've extracted seems to be pooled both formulations.


Dear VStur,

Thanks for your reply. You are right, Yijk is PK characteristic of interest. The design is 2 sequence 2 treatment and 4 period design(TRTR/RTRT). More clearly, I want to extract between subject variances (sigmaBT^2, sigmaBR^2), within subject variances (sigmaWT^2, sigmaWR^2) under the test and reference formulation, also subject-by-formulation interaction (sigmaD^2 = sigmaBT^2 + sigmaBR^2 - 2*rho*sigmaBT*sigmaBR). Is it possible to obtain all of these estimates from the summary(lmer_model), or do I need to do an extra step to obtain those values?

Thanks in advance

Edit: Full quote removed. Please delete everything from the text of the original poster which is not necessary in understanding your answer; see also this post #5! [Helmut]
 
d_labes
★★★

Berlin, Germany,
2017-02-07 12:16
(2965 d 08:03 ago)

@ StatR
Posting: # 17029
Views: 19,020
 

 Getting variance components

Post reply
Dear StatR,

❝ ... More clearly, I want to extract between subject variances (sigmaBT^2, sigmaBR^2), within subject variances (sigmaWT^2, sigmaWR^2) under the test and reference formulation, also subject-by-formulation interaction (sigmaD^2 = sigmaBT^2 + sigmaBR^2 - 2*rho*sigmaBT*sigmaBR).


What is needed for that goal is to formulate a model which have these variance-covariance parameters.
Your model is one with assuming sigmaWT = sigmaWR and sigmaBT=sigmaBR, rho=1.
Thus you have only two variance-covariance parameters in your output above:
Random effects:
 Groups   Name        Std.Dev.
 Sequence (Intercept) 0.002309 # sigmaBT = sigmaBR
 Residual             0.356224 # sigmaWT = sigmaWR

Unfortunately there is currently no way to formulate a model similar to the FDA code for replicate designs within package lme4 / lmer(). See f.i. this post for more details on the why.

In the older package nlme / lme() you may formulate the model similar to the FDA code. See that post.
But nevertheless you can't get the 'correct' CI from that R solution since the degrees of freedom for the contrast T-R are different to that of SAS Proc MIXED.

Hope this helps.
—
Regards,

Detlew
 
StatR
☆    

Turkey,
2017-02-07 12:36
(2965 d 07:43 ago)

@ d_labes
Posting: # 17030
Views: 18,945
 

 Getting variance components

Post reply

❝ What is needed for that goal is to formulate a model which have these variance-covariance parameters.

❝ Your model is one with assuming sigmaWT = sigmaWR and sigmaBT=sigmaBR, rho=1.

❝ Thus you have only two variance-covariance parameters


Dear d_labes,

Thank you very much for the explanation. At least, now, I know what should I do.

Best.
 
StatR
☆    

Turkey,
2017-02-08 09:41
(2964 d 10:39 ago)

@ d_labes
Posting: # 17033
Views: 18,909
 

 Getting variance components

Post reply
Hi again d_labes,

I tried to find the correct model by the post you referred.

❝ In the older package nlme / lme() you may formulate the model similar to the FDA code. See that post.


Unfortunately, I am confused after reading tons of post, and I could not figure out what is "trt", "tmt" ... Simply, my data is as follows

 Subject Formula Period Sequence     Yijk
       1       1      1        1 4.933014
       1       0      2        1 5.006829
       1       1      3        1 4.854162
       1       0      4        1 5.081755
       2       1      1        1 4.854738
       2       0      2        1 4.690311

Where Yijk is the PK characteristic. The first sequence represents "TRTR" and the second sequence: "RTRT". In the formula 1 = T and 0 = R. In the referred post you say that

model2 <- lme(y~ drug + prd + seq,
# this random statement fits a symmetric positive definite covariance matrix,
# identical to UN in SAS???
                random= ~ tmt-1|subj,
                #different within variabilities                 
                weights=varIdent(form = ~ 1 | drug),
                data=BlaBla, method="REML")

In your model y = Yijk, drug=Formula, prd=Period, seq=Sequence,
but what is tmt? Also, is the my data format correct for the lme? Or, Should I change Formula as 1 = T, 2 = R and Sequence as 1 = TRTR, 2 = RTRT? Could you please make me clear.

Best.
 
d_labes
★★★

Berlin, Germany,
2017-02-08 11:13
(2964 d 09:07 ago)

@ StatR
Posting: # 17034
Views: 19,181
 

 Getting variance components

Post reply
Dear StatR,

sorry for confusing you. I was a little bit sloppy in the post and mixed the factor names in Bear with mine.

To clarify my naming convention and coding:
tmt is treatment coded as "T" or "R"
period is period no
sequence is sequence coded as "TRTR" or "RTRT" (or whatever you have)
subject is subject no.
y is log-transformed PK metric (Cmax or AUC or whatever you like to evaluate)

Don't forget to make them as.factor()!

Then call
muddle <- lme(y ~ tmt + period + sequence,
              # random variance-covariance matrix
              random= ~ tmt-1|subject,
              #different within variabilities                 
              weights=varIdent(form = ~ 1 | tmt),
              data=EMAsetII, method="REML")


If you prefer your numeric coding of tmt = Formula and Sequence:
Don't forget to make them as.factor()!
Then call
muddle <- lme(y ~ Formula + Period + Sequence,
              random= ~ Formula-1|subject,
              weights=varIdent(form = ~ 1 | Formula),
              data=EMAsetII, method="REML")

Using the dataset II from the EMA Q&A you will get the following output with my convention:
Linear mixed-effects model fit by REML
  Data: EMAsetII
  Log-restricted-likelihood: 15.33728
  Fixed: y ~ tmt + period + sequence
(Intercept)        tmtT     period2     period3 sequenceRTR sequenceTRR
7.904258602 0.022391427 0.001296055 0.048118876 0.054776131 0.050923729

Random effects:
 Formula: ~tmt - 1 | subject
 Structure: General positive-definite, Log-Cholesky parametrization
         StdDev     Corr
tmtR     0.19030500 tmtR
tmtT     0.24869920 0.964
Residual 0.09296929     

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | tmt
 Parameter estimates:
       T        R
1.000000 1.237973
Number of Observations: 72
Number of Groups: 24

From that it is read (note that lme() works with SD instead of variances):
s2wT = 0.09296929^2
s2wR = (0.09296929*1.237973)^2 = 0.1150935^2
s2bT = 0.24869920^2
s2bR = 0.19030500^2
rho = 0.964 (unfortunately no more decimals)

Hope this helps.
Homework: Identify the variance-covariance terms with your coding :cool:.
—
Regards,

Detlew
 
StatR
☆    

Turkey,
2017-02-08 11:19
(2964 d 09:00 ago)

@ d_labes
Posting: # 17035
Views: 18,883
 

 Getting variance components

Post reply
Dear d_labes,

The only thing that I say: Wow! I do not know how I thanks for this reply. Thank you!

❝ sorry for confusing you. I was a little bit sloppy in the post and mixed the factor names in Bear with mine.


❝ Homework: Identify the variance-covariance terms with your coding :cool:.


:cool:
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2017-02-08 11:33
(2964 d 08:47 ago)

@ StatR
Posting: # 17036
Views: 19,078
 

 Data structure

Post reply
Hi StatR,

❝  Subject Formula Period Sequence     Yijk

❝        1       1      1        1 4.933014

❝        1       0      2        1 5.006829

❝        1       1      3        1 4.854162

❝        1       0      4        1 5.081755

❝ ...

❝

❝ Where Yijk is the PK characteristic. The first sequence represents "TRTR" and the second sequence: "RTRT". In the formula 1 = T and 0 = R.


Only a general remark about data structures in R. I suggest to code categorial (classification) variables as characters instead of numeric (Treatment, Sequence). It not only helps others to understand what you have meant without requiring an explanation (1=TRTR, 2=RTRT and 1=T, 0=R) but also models (e.g., lm(), lme(), lmer()) will automatically treat them as factors. Avoid Formula and use Treatment or Formulation instead. Personally I prefer Treatment over Formulation for generality (imagine a food-effect study where we have the same formulation).
Convert numeric variables (except the PK-response) to factors as well:

data$Subject <- as.factor(data$Subject)
data$Period  <- as.factor(data$Period)


❝ Unfortunately, I am confused after reading tons of post, and I could not figure out what is "trt", "tmt" ...


Always the treatment is meant. To add confusion in package bear Yung-jin uses drug. ;-)
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
StatR
☆    

Turkey,
2017-02-08 11:49
(2964 d 08:30 ago)

@ Helmut
Posting: # 17037
Views: 18,932
 

 Data structure

Post reply
Hi Helmut,

Many thans for your suggestions. I think you are right, otherwise it's confusing. I got the message ;-)

Best

❝ Only a general remark about data structures in R. I suggest to code categorial (classification) variables as characters instead of numeric (Treatment, Sequence). It not only helps others to understand what you have meant without requiring an explanation (1=TRTR, 2=RTRT and 1=T, 0=R) but also models (e.g., lm(), lme(), lmer()) will automatically treat them as factors. Avoid Formula and use Treatment or Formulation instead. Personally I prefer Treatment over Formulation for generality...

 
UA Flag
Activity
 Admin contact
23,409 posts in 4,921 threads, 1,669 registered users;
18 visitors (0 registered, 18 guests [including 6 identified bots]).
Forum time: 20:20 CET (Europe/Vienna)

Facts are meaningless. You could use facts to prove
anything that’s even remotely true!    Homer Simpson

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5