Diagnostics: R [Study As­sess­ment]

posted by zizou – Plzeň, Czech Republic, 2016-05-22 21:07 (2888 d 14:58 ago) – Posting: # 16348
Views: 29,558

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:

UA Flag
Activity
 Admin contact
22,987 posts in 4,824 threads, 1,662 registered users;
86 visitors (1 registered, 85 guests [including 6 identified bots]).
Forum time: 12:06 CEST (Europe/Vienna)

The only way to comprehend what mathematicians mean by Infinity
is to contemplate the extent of human stupidity.    Voltaire

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