Bioequivalence and Bioavailability Forum

Getting variance components [🇷 for BE/BA]

posted by d_labes  – Berlin, Germany, 2017-02-08 11:13 (3054 d 11:08 ago) – Posting: # 17034
Views: 19,768

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 .

Complete thread:

• Bear vs. Phoenix & SAS Helmut 2015-04-20 17:34 
  • 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 effects Astea 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 componentsd_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