[R] confidence interval too small in nlme?

From: Wittner, Ben, Ph.D. <Wittner.Ben_at_mgh.harvard.edu>
Date: Thu, 3 Jan 2008 09:41:39 -0500


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

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.



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Thu 03 Jan 2008 - 14:45:22 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 04 Jan 2008 - 09:30:05 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive