R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), 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.

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 a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("lmm-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('lmm')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "ecmeml.lmm"
> 
> ### * ecmeml.lmm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ecmeml.lmm
> ### Title: ECME algorithm for maximum-likelihood (ML) estimation in linear
> ###   mixed models
> ### Aliases: ecmeml.lmm
> ### Keywords: models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D For a detailed example, see the file "example.lmm.R" distributed
> ##D with this library.
> ## End(Not run)
> 
> 
> cleanEx(); ..nameEx <- "ecmerml.lmm"
> 
> ### * ecmerml.lmm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ecmerml.lmm
> ### Title: ECME algorithm for restricted maximum-likelihood (RML)
> ###   estimation in linear mixed models
> ### Aliases: ecmerml.lmm
> ### Keywords: models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D For a detailed example, see the file "example.lmm.R" distributed
> ##D with this library.
> ##D 
> ## End(Not run)
> 
> 
> cleanEx(); ..nameEx <- "example.lmm"
> 
> ### * example.lmm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: example.lmm
> ### Title: lmm library example command file
> ### Aliases: example.lmm
> ### Keywords: models
> 
> ### ** Examples
> 
> #
> #   -----------------------------------------------------------------
> #                         15 minutes                 90 minutes
> #                   ----------------------     ----------------------
> #                   Placebo   Low   High      Placebo   Low   High
> #   -----------------------------------------------------------------
> #   Subject 1          16     20     16          20     -6     -4       
> #           2          12     24     12          -6      4     -8
> #           3           8      8     26          -4      4      8
> #           4          20      8     NA          NA     20     -4
> #           5           8      4     -8          NA     22     -8
> #           6          10     20     28         -20     -4     -4
> #           7           4     28     24          12      8     18
> #           8          -8     20     24          -3      8    -24
> #           9          NA     20     24           8     12     NA
> #   -----------------------------------------------------------------
> #
> ########################################################################
> # Below we show how to fit a traditional compound symmetry model
> # with a fixed effect for each column (occasion) and a random 
> # intercept for each subject. First we enter the data.
> #
> y <- c(16,20,16,20,-6,-4,
+     12,24,12,-6,4,-8,
+     8,8,26,-4,4,8,
+     20,8,20,-4,
+     8,4,-8,22,-8,
+     10,20,28,-20,-4,-4,
+     4,28,24,12,8,18,
+     -8,20,24,-3,8,-24,
+     20,24,8,12)
> occ <- c(1,2,3,4,5,6,
+       1,2,3,4,5,6,
+       1,2,3,4,5,6,
+       1,2,5,6,
+       1,2,3,5,6,
+       1,2,3,4,5,6,
+       1,2,3,4,5,6,
+       1,2,3,4,5,6,
+       2,3,4,5)
> subj <- c(1,1,1,1,1,1,
+        2,2,2,2,2,2,
+        3,3,3,3,3,3,
+        4,4,4,4,
+        5,5,5,5,5,
+        6,6,6,6,6,6,
+        7,7,7,7,7,7,
+        8,8,8,8,8,8,
+        9,9,9,9)
> ########################################################################
> # Now we must specify the model. 
> # If the six measurements per subject were ordered in time, we might
> # consider using a model with time of measurement entered with linear
> # (or perhaps higher-order polynomial) effects. But because the
> # six measurements are not clearly ordered, let's use a model that has
> # an intercept and five dummy codes to allow the population means for
> # the six occasions to be estimated freely. We will also allow the
> # intercept to randomly vary by subject. For a subject i with no
> # missing values, the covariate matrices will be
> #
> #                   1 1 0 0 0 0              1
> #                   1 0 1 0 0 0              1
> #           Xi =    1 0 0 1 0 0       Zi =   1
> #                   1 0 0 0 1 0              1
> #                   1 0 0 0 0 1              1
> #                   1 0 0 0 0 0              1
> #
> # The Xi's and Zi's are combined into a single matrix called
> # pred. The pred matrix has length(y) rows. Each column of Xi and Zi
> # must be represented in pred. Because Zi is merely the first column
> # of Xi, we do not need to enter that column twice. So pred is simply
> # the matrices Xi (i=1,...,9), stacked upon each other.
> #
> pred <- cbind(int=rep(1,49),dummy1=1*(occ==1),dummy2=1*(occ==2),
+     dummy3=1*(occ==3),dummy4=1*(occ==4),dummy5=1*(occ==5))
> xcol <- 1:6
> zcol <- 1
> ########################################################################
> # Now find ML estimates using the ECME procedure and the faster
> # scoring algorithm
> #
> ecmeml.result <- ecmeml.lmm(y,subj,pred,xcol,zcol)
Performing ECME...
> fastml.result <- fastml.lmm(y,subj,pred,xcol,zcol)
Performing FAST-ML...
> #
> # In this example, ECME converged in 212 cycles, but the fast
> # algorithm took only 8. The results can be viewed by printing the
> # various components of "ecmeml.result" and "fastml.result". 
> # For example, extract the ML estimate of the fixed effects beta.
> #
> beta.hat <- fastml.result$beta
> #
> # Because of the dummy codes used in the Xi's, the first element of
> # beta (the intercept) estimates the mean for the last occasion,
> # and the other elements of beta estimate the differences in means
> # between the first five occasions and the last one. So we can find
> # the estimated means for the six occasions like this:
> #
> muhat <- c(beta.hat[2]+beta.hat[1], beta.hat[3]+beta.hat[1],
+    beta.hat[4]+beta.hat[1], beta.hat[5]+beta.hat[1],
+    beta.hat[6]+beta.hat[1], beta.hat[1]) 
> #
> # The functions for RML estimation work exactly the same way:
> #
> ecmerml.result <- ecmerml.lmm(y,subj,pred,xcol,zcol)
Performing ECME for RML estimation...
> fastrml.result <- fastrml.lmm(y,subj,pred,xcol,zcol)
Performing FAST-RML...
> #
> #######################################################################
> # The function "fastrml.lmm" calculates the improved variance
> # estimates for random effects described in Section 4 of Schafer
> # (1998). The code below reproduces Table 2, which compares 
> # 95% interval estimates under the new method to conventional
> # empirical Bayes intervals.
> #
> b.hat <- as.vector(fastrml.result$b.hat)
> se.new <- sqrt(as.vector(fastrml.result$cov.b.new))
> se.old <- sqrt(as.vector(fastrml.result$cov.b))
> table2 <- cbind(round(b.hat,3),
+    round(cbind(b.hat-2*se.old,b.hat+2*se.old,
+       b.hat-2*se.new,b.hat+2*se.new),2),
+    round(100*(se.new-se.old)/se.old))
> dimnames(table2) <- list(paste("Subject",format(1:9)),
+    c("Est.","Lower.old","Upper.old","Lower.new","Upper.new",
+    "Increase (%)"))
> print(table2)
            Est. Lower.old Upper.old Lower.new Upper.new Increase (%)
Subject 1  0.190     -2.39      2.77     -3.21      3.59           32
Subject 2 -0.168     -2.75      2.41     -3.43      3.09           26
Subject 3  0.011     -2.57      2.59     -2.59      2.61            1
Subject 4  0.215     -2.40      2.83     -3.46      3.89           40
Subject 5 -0.459     -3.06      2.14     -6.50      5.58          133
Subject 6 -0.288     -2.87      2.29     -4.54      3.96           65
Subject 7  0.668     -1.91      3.25     -7.52      8.86          218
Subject 8 -0.482     -3.06      2.10     -6.68      5.71          140
Subject 9  0.313     -2.31      2.93     -4.27      4.90           75
> #
> #######################################################################
> # The functions "mgibbs.lmm" and "fastmcmc.lmm" perform the MCMC
> # procedures described in Section 5. The code below runs each
> # algorithm for 5,000 cycles, and then displays autocorrelation
> # plots like those of Figure 1.
> #
> prior <- list(a=3*100,b=3,c=3,Dinv=3*5)
> gibbs.result <- mgibbs.lmm(y,subj,pred,xcol,zcol,prior=prior,
+    seed=1234,iter=5000)
Performing modified Gibbs...
> fmcmc.result <- fastmcmc.lmm(y,subj,pred,xcol,zcol,prior=prior,
+    seed=2345,iter=5000)
Performing FAST-MODE...
Performing FAST-MCMC...
> # 
> # Before doing this, make sure that a graphics device is active:
> # library(ts)
> par(mfrow=c(2,1))
> #
> # there were problems with debian Linux; so add ylim
> acf(log(gibbs.result$psi.series[1,1,]),lag.max=10, ylim=0:1)
> acf(log(fmcmc.result$psi.series[1,1,]),lag.max=10, ylim=0:1)
> 
> #######################################################################
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "fastmcmc.lmm"
> 
> ### * fastmcmc.lmm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fastmcmc.lmm
> ### Title: Rapidly converging Markov chain Monte Carlo algorithm for
> ###   Bayesian inference in linear mixed models
> ### Aliases: fastmcmc.lmm
> ### Keywords: models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D For a detailed example, see the file "example.lmm.R" distributed
> ##D with this library.
> ##D 
> ## End(Not run)
> 
> 
> cleanEx(); ..nameEx <- "fastml.lmm"
> 
> ### * fastml.lmm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fastml.lmm
> ### Title: Rapidly converging algorithm for maximum-likelihood (ML)
> ###   estimation in linear mixed models
> ### Aliases: fastml.lmm
> ### Keywords: models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D For a detailed example, see the file "example.lmm.R" distributed
> ##D with this library.
> ##D 
> ## End(Not run)
> 
> 
> cleanEx(); ..nameEx <- "fastmode.lmm"
> 
> ### * fastmode.lmm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fastmode.lmm
> ### Title: Rapidly converging algorithm for calculating posterior modes in
> ###   linear mixed models
> ### Aliases: fastmode.lmm
> ### Keywords: models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D For a detailed example, see the file "example.lmm.R" distributed
> ##D with this library.
> ##D 
> ## End(Not run)
> 
> 
> cleanEx(); ..nameEx <- "fastrml.lmm"
> 
> ### * fastrml.lmm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fastrml.lmm
> ### Title: Rapidly converging algorithm for restricted maximum-likelihood
> ###   (RML) estimation in linear mixed models
> ### Aliases: fastrml.lmm
> ### Keywords: models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D For a detailed example, see the file "example.lmm.R" distributed
> ##D with this library.
> ##D 
> ## End(Not run)
> 
> 
> cleanEx(); ..nameEx <- "mgibbs.lmm"
> 
> ### * mgibbs.lmm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mgibbs.lmm
> ### Title: Modified Gibbs sampler for Bayesian inference in linear mixed
> ###   models
> ### Aliases: mgibbs.lmm
> ### Keywords: models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D For a detailed example, see the file "example.lmm.R" distributed
> ##D with this library.
> ##D 
> ## End(Not run)
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  0.75 0.01 0.77 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')