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.
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: med.kuleuven.be/biostat 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}}
-- 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: med.kuleuven.be/biostat 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}}