Re: [R] confidence interval too small in nlme?

From: Dimitris Rizopoulos <dimitris.rizopoulos_at_med.kuleuven.be>
Date: Fri, 4 Jan 2008 09:29:40 +0100

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



> 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}}



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 Fri 04 Jan 2008 - 08:33:03 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 - 10: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