Bioequivalence and Bioavailability Forum • lme() does not work with all fixed effects

lme() does not work with all fixed effects [🇷 for BE/BA]

posted by Astea – Russia, 2016-11-04 01:13 (3060 d 18:12 ago) – Posting: # 16771
Views: 20,856

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 T_mean=7,830303 and R_mean=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"

Complete thread:

Bear vs. Phoenix & SAS Helmut 2015-04-20 17:34 [🇷 for BE/BA]
- R vs. Phoenix & SAS? yjlee168 2015-04-20 19:36
  - R vs. Phoenix & SAS? Helmut 2015-04-21 01:02
    - lme() does not work with all fixed effects yjlee168 2015-04-21 23:41
      - lme() does not work with all fixed effectsAstea 2016-11-04 00:13
        
        lmer: Method B (PE catched for imbalanced dataset!!!) and Method C mittyri 2016-11-05 17:38
        
        lmer: Method B (PE catched for imbalanced dataset!!!) and Method C Astea 2016-11-05 19:27
        
        lmer: Method B is ready for scaling mittyri 2016-11-05 20:01
        
        lmer: Method B is ready for scaling Astea 2016-11-06 11:50
        
        lmer: Method B is ready for scaling mittyri 2016-11-07 06:07
- info for lsmeans yjlee168 2015-04-20 21:34
  - info for lsmeans Helmut 2015-04-21 01:15
    - once more about R and replicate designes Astea 2016-11-02 23:43
      - once more about R and replicate designes VStus 2016-11-06 11:34
        
        Getting variance components from the lmer output StatR 2017-02-03 13:53
        
        Getting variance components from the lmer output VStus 2017-02-03 15:47
        
        Getting variance components from the lmer output StatR 2017-02-03 17:12
        
        Getting variance components d_labes 2017-02-07 11:16
        
        Getting variance components StatR 2017-02-07 11:36
        Getting variance components StatR 2017-02-08 08:41
        
        Getting variance components d_labes 2017-02-08 10:13
        
        Getting variance components StatR 2017-02-08 10:19
        
        Data structure Helmut 2017-02-08 10:33
        
        Data structure StatR 2017-02-08 10:49