Diagnostics: R [Study As­sess­ment]

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

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,138 posts in 4,246 threads, 1,387 registered users;
online 55 (0 registered, 55 guests [including 5 identified bots]).
Forum time (Europe/Vienna): 15:50 CET

The fundamental cause of trouble in the world today is
that the stupid are cocksure
while the intelligent are full of doubt.    Bertrand Russell

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