R version 2.6.0 Patched (2007-10-17 r43203)
Copyright (C) 2007 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require(lme4)
Loading required package: lme4
Loading required package: Matrix
Loading required package: lattice
> if (!interactive())
+     pdf(file = "profiled.pdf", width = 10, height = 5)
> 
> sessionInfo()
R version 2.6.0 Patched (2007-10-17 r43203) 
x86_64-unknown-linux-gnu 

locale:
LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] lme4_0.999375-0   Matrix_0.999375-3 lattice_0.17-1   

loaded via a namespace (and not attached):
[1] grid_2.6.0
> 
> ## Illustrate problem with mcmcsamp for two variance components mixed-effects model.
> 
> make.data <- function(N.batch)
+ {
+   ## Create example data set
+   ## N.batch Batches, 15 units per batch = (5 Trt A, 5 Trt B, 5 Trt C)
+ 
+ 
+   ## Trt means, Batch SD, and residual SD
+   theta <- c(5, 10, 15)
+   sd.batch <- 0.5
+   sd.resid <- 1.5
+ 
+   d <- expand.grid( Batch=factor(1:N.batch), Trt=factor(c("A", "B", "C")), Units=1:5 )
+   d <- d[order(d$Batch, d$Trt),][,-3]
+   d$y <- theta[ unclass(d$Trt) ] + rep( rnorm( N.batch, sd=sd.batch ), each=15 ) + rnorm( nrow(d), sd=sd.resid )
+   
+   d
+ }
> 
> ## Values of theta
> th1 <- exp(seq(-6.1, 0.8, len = 201))
> 
> ## Create the response variable and fit the linear mixed model
> set.seed(991); d04 <- make.data(N.batch=4)
> fit04 <- lmer(y ~ Trt + (1|Batch), d04, method = "ML")
> ML04 <- devmat(fit04, th1)[["ML"]] - deviance(fit04)
> 
> ## Repeat for 10 batches
> set.seed(991); d10 <- make.data(N.batch=10)
> fit10 <- lmer(y ~ Trt + (1|Batch), d10, method = "ML")
> ML10 <- devmat(fit10, th1)[["ML"]] - deviance(fit10)
> 
> ## Repeat for 25 batches
> set.seed(991); d25 <- make.data(N.batch=25)
> fit25 <- lmer(y ~ Trt + (1|Batch), d25, method = "ML")
> ML25 <- devmat(fit25, th1)[["ML"]] - deviance(fit25)
> 
> ## Plot all three curves.  The vertical scale is based on the 0.99
> ## quantile of a chisquared distribution with 1 degree of freedom.
> ## Plateaus below this level can be expected to be reached with a high
> ## probability in a long chain.
> print(xyplot(ML04 + ML10 + ML25 ~ log(th1), type = c("g", "l"),
+              aspect = 0.4, xlab = "Logarithm of relative standard deviation",
+              ylab = "Change in the profiled deviance",
+              ylim = qchisq(0.99, 1) * c(-0.035, 1.035),
+              scales = list(x = list(axs = 'i')),
+              auto.key = list(text = as.character(c(4, 10, 25)),
+              lines = TRUE, points = FALSE,
+              space = "top", columns = 3)))
> 
> 
> proc.time()
   user  system elapsed 
  8.068   0.204   8.420