[R] confidence interval too small in nlme?

[Dimitris Rizopoulos]

-- sorry but the code I posted yesterday wasn't self-contained; here 
it is again --

well, these are *approximate* confidence intervals (i.e., big enough 
sample sizes are required for the asympotics to work), check Section 
2.4.3 in Pinheiro and Bates (2000), and also the code below

library('nlme')

M <- 6
n <- 3
beta <- 67
sigma.b <- 25
sigma <- 4
Rail <- rep(1:M, each=n)
set.seed(56820)
B <- 10000
tvals <- numeric(B)
num.wrong <- 0
for (K in 1:B) {
    travel <- beta + rep(rnorm(M, sd = sigma.b), each = n) + 
rnorm(M*n, sd = sigma)
    fm1Rail.lme <- lme(travel ~ 1, random = ~ 1 | Rail)
    tvals[K] <- (fixef(fm1Rail.lme) - beta) / sqrt(fm1Rail.lme$varFix)
    CI <- intervals(fm1Rail.lme, which = "fixed")$fixed
    if (CI[1, "lower"] > beta || CI[1, "upper"] < beta)
        num.wrong <- num.wrong + 1
}

num.wrong / B
# this is based on the empirical distribution
quantile(tvals, c(0.025, 0.975))
# this is based on the asympotic distribution
qt(c(0.025, 0.975), 12)


I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm



----- Original Message ----- 
From: "Wittner, Ben, Ph.D." <Wittner.Ben at mgh.harvard.edu>
To: <r-help at r-project.org>
Sent: Thursday, January 03, 2008 3:41 PM
Subject: [R] confidence interval too small in nlme?


> Hello,
>
> I am interested in using nlme to model repeated measurements, but I 
> don't seem
> to get good CIs.
>
> With the code below I tried to generate data sets according to the 
> model given
> by equations (1.4) and (1.5) on pages 7 and 8 of Pinheiro and Bates 
> 2000 (having
> chosen values for beta, sigma.b and sigma similar to those estimated 
> in the
> text).
> For each data set I used lme() to fit a model, used intervals() to 
> get a 95% CI
> for beta, and then checked whether the the CI contained beta.
> The rate at which the CI did not contain beta was 8%, which was 
> greater than the
> 5% I was expecting.
> This may seem like a small difference, but in the lab in which I 
> work M would
> more likely be 2 or 3. When I re-ran with M = 3 I got 13% of the CIs 
> not
> containing beta and when I re-ran with M = 2, I got 21%.
>
> Am I calculating the CIs incorrectly?
> Am I interpreting them incorrectly?
> Am I doing anything else wrong?
>
> Output of packageDescription('nlme') and version given below the 
> code.
>
> Any help will be greatly appreciated. Thanks very much in advance.
> -Ben
>
> #########################################################################
> ##
> ##  Code to test intervals() based on equations (1.4) and (1.5) of
> ##  Pinheiro and Bates
> ##
>
> library('nlme')
>
> M <- 6
> n <- 3
> beta <- 67
> sigma.b <- 25
> sigma <- 4
>
> Rail <- rep(1:M, each=n)
>
> set.seed(56820)
> B <- 10000
> num.wrong <- 0
> error.fraction <- Ks <- c()
> for (K in 1:B) {
>  travel <- beta + rep(rnorm(M, sd=sigma.b), each=n) + rnorm(M*n, 
> sd=sigma)
>  fm1Rail.lme <- lme(travel ~ 1, random = ~ 1 | Rail)
>  CI <- intervals(fm1Rail.lme, which='fixed')$fixed
>  if ((CI[1, 'lower'] > beta) || (CI[1, 'upper'] < beta))
>    num.wrong <- num.wrong + 1
>  if (K %% 200 == 0) {
>    error.fraction <- c(error.fraction, num.wrong/K)
>    Ks <- c(Ks, K)
>    plot(Ks, error.fraction, type='b',
>         ylim=range(c(0, 0.05, error.fraction)))
>    abline(h=0.05, lty=3)
>  }
> }
>
> num.wrong/B
>
> #########################################################################
> ##
> ##  version information
> ##
>
>> packageDescription('nlme')
> Package: nlme
> Version: 3.1-86
> Date: 2007-10-04
> Priority: recommended
> Title: Linear and Nonlinear Mixed Effects Models
> Author: Jose Pinheiro <Jose.Pinheiro at pharma.novartis.com>, Douglas
>        Bates <bates at stat.wisc.edu>, Saikat DebRoy
>        <saikat at stat.wisc.edu>, Deepayan Sarkar
>        <Deepayan.Sarkar at R-project.org> the R Core team.
> Maintainer: R-core <R-core at R-project.org>
> Description: Fit and compare Gaussian linear and nonlinear
>        mixed-effects models.
> Depends: graphics, stats, R (>= 2.4.0)
> Imports: lattice
> LazyLoad: yes
> LazyData: yes
> License: GPL (>=2)
> Packaged: Thu Oct 4 23:25:21 2007; hornik
> Built: R 2.6.0; i686-pc-linux-gnu; 2007-12-26 15:48:00; unix
>
> -- File: /home/bwittner/R-2.6.0/library/nlme/DESCRIPTION
>> version
>               _
> platform       i686-pc-linux-gnu
> arch           i686
> os             linux-gnu
> system         i686, linux-gnu
> status
> major          2
> minor          6.0
> year           2007
> month          10
> day            03
> svn rev        43063
> language       R
> version.string R version 2.6.0 (2007-10-03)
>
> The information transmitted in this electronic communi...{{dropped:25}}