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
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." <[EMAIL PROTECTED]>
To: <r-help@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 <[EMAIL PROTECTED]>, Douglas
> Bates <[EMAIL PROTECTED]>, Saikat DebRoy
> <[EMAIL PROTECTED]>, Deepayan Sarkar
> <[EMAIL PROTECTED]> the R Core team.
> Maintainer: R-core <[EMAIL PROTECTED]>
> 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@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.
