[R] confidence interval too small in nlme?

[Wittner, Ben, Ph.D.]

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 communication is intended only
for the person or entity to whom it is addressed and may contain confidential
and/or privileged material. Any review, retransmission, dissemination or other
use of or taking of any action in reliance upon this information by persons or
entities other than the intended recipient is prohibited. If you received this
information in error, please contact the Compliance HelpLine at 800-856-1983 and
properly dispose of this information.