Diagnostics: R [Study As­sess­ment]

posted by zizou – Plzeň, Czech Republic, 2016-05-22 17:07  – Posting: # 16348
Views: 18,406

Dear Helmut.
» Estimates and their SEs are exactly the same. CIs are not (due to different DFs?).
Exactly, different DFs.

# equivalent of your code:
library(lme4)
Subj   <- c(1, 2, 4, 5, 6, 4, 5, 6, 7, 8, 9, 7, 8, 9)
Dose   <- c(25, 25, 50, 50, 50, 250, 250, 250, 75, 75, 75, 250, 250, 250)
AUC    <- c(326.40, 437.82, 557.47, 764.85, 943.59, 2040.84, 2989.29,
            4107.58, 1562.42, 982.02, 1359.68, 3848.86, 4333.10, 3685.55)
Cmax   <- c(64.82, 67.35, 104.15, 143.12, 243.63, 451.44, 393.45,
            796.57, 145.13, 166.77, 296.90, 313.00, 387.00, 843.00)
resp   <- data.frame(Subj, Dose, Cmax, AUC)
muddle <- lmer(log(Cmax) ~ log(Dose) + (1|Subj), data=resp)
summary(muddle)


DF, p-values, CIs not reported in lmer for the reasons stated there.

For more decimal places:
print(muddle, digits=7, ranef.comp=c("Var","Std.Dev."))
# REML criterion at convergence: 6.2435 (same as in SPSS 6.243548)
# note: -2 (Restricted) Log Likelihood from R lme: (-2)*(-3.121774) = 6.243548


Some reference for lmer can be found in this PDF.

For DFs, p-values, ...: library lmerTest can be used, see this PDF (page 2 - description, details, references to SAS).
library(lmerTest) # Library includes modification of function lmer (if I am not mistaken).
muddle <- lmer(log(Cmax) ~ log(Dose) + (1|Subj), data=resp)
summary(muddle)


Results:
Linear mixed model fit by REML t-tests use Satterthwaite approximations to
  degrees of freedom [lmerMod] Formula: log(Cmax) ~ log(Dose) + (1 | Subj)
   Data: resp

REML criterion at convergence: 6.2

Scaled residuals:
     Min       1Q   Median       3Q      Max
-1.07548 -0.35579 -0.03301  0.45089  0.91854

Random effects:
 Groups   Name        Variance Std.Dev.
 Subj     (Intercept) 0.11205  0.3347 
 Residual             0.01456  0.1207 
Number of obs: 14, groups:  Subj, 8

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)   
(Intercept)  1.94139    0.24314 9.19600   7.985 1.98e-05 ***
log(Dose)    0.76174    0.04728 5.89600  16.111 4.24e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr)
log(Dose) -0.863

I think there is a little bug in df visualization - rounded on 3 decimal places, visible on 5 (actually I have not the latest version of R).
summary(muddle)$coefficients["(Intercept)","df"]  # [1] 9.195607

For calculation of 90% CIs: (below is only example for Intercept lower and upper limit)
alpha=0.05
summary(muddle)$coefficients["(Intercept)","Estimate"] - qt(1-alpha,summary(muddle)$coefficients["(Intercept)","df"]) * summary(muddle)$coefficients["(Intercept)","Std. Error"]
summary(muddle)$coefficients["(Intercept)","Estimate"] + qt(1-alpha,summary(muddle)$coefficients["(Intercept)","df"]) * summary(muddle)$coefficients["(Intercept)","Std. Error"]


Best regards,
zizou

Complete thread:

Activity
 Admin contact
20,122 posts in 4,245 threads, 1,383 registered users;
online 19 (1 registered, 18 guests [including 7 identified bots]).
Forum time (Europe/Vienna): 07:33 UTC

Genius is that which forces
the inertia of humanity to learn.    Henri Bergson

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