R version 2.8.0 Under development (unstable) (2008-04-08 r45173)
Copyright (C) 2008 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.

> library(lme4)
Loading required package: Matrix
Loading required package: lattice

Attaching package: 'Matrix'


	The following object(s) are masked from package:stats :

	 xtabs 

> Dyestuff <-
+     data.frame(Batch = gl(6, 5, lab = LETTERS[1:6]),
+                Yield = c(
+                1545, 1440, 1440, 1520, 1580,
+                1540, 1555, 1490, 1560, 1495,
+                1595, 1550, 1605, 1510, 1560,
+                1445, 1440, 1595, 1465, 1545,
+                1595, 1630, 1515, 1635, 1625,
+                1520, 1455, 1450, 1480, 1445))
> str(Dyestuff)
'data.frame':	30 obs. of  2 variables:
 $ Batch: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 2 2 2 2 2 ...
 $ Yield: num  1545 1440 1440 1520 1580 ...
> dotplot(reorder(Batch, Yield) ~ Yield, Dyestuff,
+         ylab = "Batch", jitter.y = TRUE, aspect = 0.3,
+         type = c("p", "a"))
> (fm1 <- lmer(Yield ~ 1|Batch, Dyestuff))
Linear mixed model fit by REML 
Formula: Yield ~ 1 | Batch 
   Data: Dyestuff 
   AIC   BIC logLik deviance REMLdev
 325.7 329.9 -159.8    327.4   319.7
Random effects:
 Groups   Name        Variance Std.Dev.
          (Intercept) 1763.7   41.996  
 Residual             2451.3   49.511  
Number of obs: 30, groups: Batch, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept)  1527.50      19.38   78.81
> set.seed(4376521)
> mcs1 <- mcmcsamp(fm1, 10000)
> xyplot(mcs1)
> HPDinterval(mcs1)
$fixef
               lower    upper
(Intercept) 1485.419 1571.323
attr(,"Probability")
[1] 0.95

$ST
            lower    upper
[1,] 7.663499e-10 1.706135
attr(,"Probability")
[1] 0.95

$sigma
        lower    upper
[1,] 37.87020 74.35672
attr(,"Probability")
[1] 0.95

> 
> proc.time()
   user  system elapsed 
 10.316   0.164  10.512