R version 2.6.0 (2007-10-03)
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.

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.

> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## This plot.new() patch has no effect yet for persp();
> ## layout() & filled.contour() are now ok
> assign("plot.new", function() { .Internal(plot.new())
+ 		       pp <- par(c("mfg","mfcol","oma","mar"))
+ 		       if(all(pp$mfg[1:2] == c(1, pp$mfcol[2]))) {
+ 			 outer <- (oma4 <- pp$oma[4]) > 0; mar4 <- pp$mar[4]
+ 			 mtext(paste("help(",..nameEx,")"), side = 4,
+ 			       line = if(outer)max(1, oma4 - 1) else min(1, mar4 - 1),
+ 			       outer = outer, adj=1, cex= .8, col="orchid")} },
+        env = .CheckExEnv)
> assign("cleanEx", function(env = .GlobalEnv) {
+ 	rm(list = ls(envir = env, all.names = TRUE), envir = env)
+ 	RNGkind("Wichmann-Hill", "Kinderman-Ramage")
+ 	set.seed(290875)
+ 	#	assign(".Random.seed", c(0,rep(7654,3)), pos=1)
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) #-- for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> postscript("multcomp-Examples.ps")
> assign("par.postscript", par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> library('multcomp')
Loading required package: mvtnorm
> cleanEx(); ..nameEx <- "MultipleEndpoints"
> ###--- >>> `MultipleEndpoints' <<<----- Multiple Endpoints Data Set
> 
> 	## alias	 help(MultipleEndpoints)
> 
> ##___ Examples ___:
> 
> cleanEx(); ..nameEx <- "angina"
> ###--- >>> `angina' <<<----- Dose Response Data Set
> 
> 	## alias	 help(angina)
> 
> ##___ Examples ___:
> 
> load("angina.rda")
> 
> # perform a dose-response analysis using simultaneous confidence 
> # intervals for Willimas' contrasts
> summary(simint(response~dose, data=angina, alternative="greater",
+                type="Williams"))

	 Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: Williams Contrasts


Fit: lm(formula = response ~ dose, data = angina)

Estimated Quantile = 1.978

Linear Hypotheses:
         Estimate lwr     upr    
C 1 <= 0 10.4990   7.4351     Inf
C 2 <= 0  7.7470   5.0936     Inf
C 3 <= 0  6.2970   3.7953     Inf
C 4 <= 0  5.2465   2.8243     Inf

95% family-wise confidence level
 

Warning message:
'simint.default' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> 
> # compute now adjusted p-values for McDermott's test on trend
> summary(simtest(response~dose, data=angina, type="McDermott",
+                 alternative="greater",ttype="logical"))

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: McDermott Contrasts


Fit: lm(formula = response ~ dose, data = angina)

Linear Hypotheses:
         Estimate Std. Error t value p value    
C 1 <= 0    2.095      1.549   1.353  0.3161    
C 2 <= 0    2.349      1.341   1.751  0.1609    
C 3 <= 0    3.164      1.265   2.502  0.0315 *  
C 4 <= 0    7.877      1.225   6.433  <1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported)

Warning message:
'simtest.default' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> 
> ## Keywords: 'datasets'.
> 
> 
> cleanEx(); ..nameEx <- "cholesterol"
> 
> ### * cholesterol
> 
> ### Name: cholesterol
> ### Title: Cholesterol Reduction Data Set
> ### Aliases: cholesterol
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(cholesterol)
> 
> # adjusted p-values for all-pairwise comparisons in a one-way layout 
> # tests for restricted combinations
> simtest(response ~ trt, data=cholesterol, type="Tukey",
+         ttype="logical")

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = response ~ trt, data = cholesterol)

Linear Hypotheses:
                     Estimate Std. Error t value p value    
2times - 1time == 0     3.443      1.443   2.385 0.13816    
4times - 1time == 0     6.593      1.443   4.568 < 0.001 ***
drugD - 1time == 0      9.579      1.443   6.637 < 0.001 ***
drugE - 1time == 0     15.166      1.443  10.507 < 0.001 ***
4times - 2times == 0    3.150      1.443   2.182 0.20503    
drugD - 2times == 0     6.136      1.443   4.251 < 0.001 ***
drugE - 2times == 0    11.723      1.443   8.122 < 0.001 ***
drugD - 4times == 0     2.986      1.443   2.069 0.25124    
drugE - 4times == 0     8.573      1.443   5.939 < 0.001 ***
drugE - drugD == 0      5.586      1.443   3.870 0.00307 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported)

Warning message:
'simtest.default' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> 
> # adjusted p-values all-pairwise comparisons in a one-way layout 
> # (tests for free combinations -> p-values will be larger)
> simtest(response ~ trt, data=cholesterol, type="Tukey",
+         ttype="free")

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = response ~ trt, data = cholesterol)

Linear Hypotheses:
                     Estimate Std. Error t value p value    
2times - 1time == 0     3.443      1.443   2.385 0.13805    
4times - 1time == 0     6.593      1.443   4.568 < 0.001 ***
drugD - 1time == 0      9.579      1.443   6.637 < 0.001 ***
drugE - 1time == 0     15.166      1.443  10.507 < 0.001 ***
4times - 2times == 0    3.150      1.443   2.182 0.20504    
drugD - 2times == 0     6.136      1.443   4.251 < 0.001 ***
drugE - 2times == 0    11.723      1.443   8.122 < 0.001 ***
drugD - 4times == 0     2.986      1.443   2.069 0.25129    
drugE - 4times == 0     8.573      1.443   5.939 < 0.001 ***
drugE - drugD == 0      5.586      1.443   3.870 0.00313 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported)

Warning message:
'simtest.default' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> 
> # the following lines illustrate the basic principles of
> # parameter estimation used in all functions in this package
> # and how the low-level functions can be used with raw parameter
> # estimates.
> 
> # the full design matrix (with reduced rank!)
> x <- cbind(1, 
+            matrix(c(rep(c(rep(1,10), rep(0,50)), 4), 
+                     rep(1, 10)), nrow = 50))
> y <- cholesterol$response
> 
> xpxi   <- multcomp:::MPinv(t(x) %*% x)$MPinv
> rankx  <- sum(diag((xpxi %*% (t(x) %*% x))))
> n      <- nrow(x)
> p      <- ncol(x)
> df     <- round(n-rankx)
> 
> # parameter estimates and their correlation
> parm   <- xpxi %*% t(x) %*% y
> mse    <- t(y-x %*% parm) %*% (y-x %*% parm)/df
> covm   <- mse[1,1]*xpxi
> 
> # the contrast matrix
> contrast <- contrMat(table(cholesterol$trt), type="Tukey")
> 
> # use the work-horse directly (and add zero column for the intercept)
> 
> csimint(estpar=parm, df=df, covm=covm, cmatrix=cbind(0, contrast))    

	 Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: user-defined Contrasts


Fit: NULL

Estimated Quantile = 2.8416

Linear Hypotheses:
                     Estimate lwr     upr    
2times - 1time == 0   3.4430  -0.6584  7.5444
4times - 1time == 0   6.5928   2.4914 10.6943
drugD - 1time == 0    9.5792   5.4778 13.6806
drugE - 1time == 0   15.1655  11.0641 19.2670
4times - 2times == 0  3.1498  -0.9516  7.2513
drugD - 2times == 0   6.1362   2.0348 10.2376
drugE - 2times == 0  11.7226   7.6211 15.8240
drugD - 4times == 0   2.9864  -1.1151  7.0878
drugE - 4times == 0   8.5727   4.4713 12.6742
drugE - drugD == 0    5.5864   1.4849  9.6878

95% family-wise confidence level
 

Warning message:
'csimint' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> csimtest(estpar=parm, df=df, covm=covm, cmatrix=cbind(0, contrast),
+          ttype="logical")      

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: user-defined Contrasts


Fit: NULL

Linear Hypotheses:
                     Estimate Std. Error t value  p value    
2times - 1time == 0     3.443      1.443   2.385 0.138175    
4times - 1time == 0     6.593      1.443   4.568 0.000342 ***
drugD - 1time == 0      9.579      1.443   6.637  < 1e-04 ***
drugE - 1time == 0     15.166      1.443  10.507  < 1e-04 ***
4times - 2times == 0    3.150      1.443   2.182 0.205041    
drugD - 2times == 0     6.136      1.443   4.251 0.000977 ***
drugE - 2times == 0    11.723      1.443   8.122  < 1e-04 ***
drugD - 4times == 0     2.986      1.443   2.069 0.251301    
drugE - 4times == 0     8.573      1.443   5.939  < 1e-04 ***
drugE - drugD == 0      5.586      1.443   3.870 0.003060 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported)

Warning message:
'csimtest' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> 
> cleanEx(); ..nameEx <- "contrMat"
> 
> ## Keywords: 'datasets'.
> 
> 
> data(detergent)
> 
> N <- rep(2, 5)
> 
> # BIBD: prepare the contrast matrix = all-pair comparisons for 
> # the 5 levels of detergent
> C <- contrMat(N, type="Tukey")
> # the additional 10 columns of are for the 10 blocks
> C <- cbind( matrix(0, ncol=10, nrow=10), C )
> # numerate the contrasts
> colnames(C) <- NULL
> rownames(C) <- paste("C", 1:nrow(C), sep="")
> 
> # adjusted p-values
> summary(simtest(plates ~ block+detergent, data=detergent,
+         cmatrix = list(detergent = contrMat(table(detergent$detergent), type = "Tukey"))))

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: lm(formula = plates ~ block + detergent, data = detergent)

Linear Hypotheses:
           Estimate Std. Error t value p value    
B - A == 0  -2.1333     0.8679  -2.458 0.15020    
C - A == 0   3.6000     0.8679   4.148 0.00586 ** 
D - A == 0   2.2000     0.8679   2.535 0.13155    
E - A == 0  -4.3333     0.8679  -4.993 0.00111 ** 
C - B == 0   5.7333     0.8679   6.606 < 0.001 ***
D - B == 0   4.3333     0.8679   4.993 0.00103 ** 
E - B == 0  -2.2000     0.8679  -2.535 0.13165    
D - C == 0  -1.4000     0.8679  -1.613 0.51052    
E - C == 0  -7.9333     0.8679  -9.140 < 0.001 ***
E - D == 0  -6.5333     0.8679  -7.527 < 0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported)

Warning message:
'simtest.default' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> # whichf="detergent", type="Tukey", ttype="logical")) # , cmatrix=C))
> 
> 
> ## Keywords: 'datasets'.
> 
> 
> cleanEx(); ..nameEx <- "recovery"
> ###--- >>> `recovery' <<<----- Recovery Time Data Set
> 
> 	## alias	 help(recovery)
> 
> ##___ Examples ___:
> 
> data(recovery)
> 
> # one-sided simultaneous confidence intervals for Dunnett 
> # in the one-way layout
> simint(minutes~blanket, data=recovery, conf.level=0.9, 
+        alternative="less",eps=0.0001)

	 Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = minutes ~ blanket, data = recovery)

Estimated Quantile = -1.843

Linear Hypotheses:
             Estimate lwr      upr     
b1 - b0 >= 0 -2.13333     -Inf  0.82239
b2 - b0 >= 0 -7.46667     -Inf -4.51094
b3 - b0 >= 0 -1.66667     -Inf -0.03606

90% family-wise confidence level
 

Warning message:
'simint.default' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> 
> # same results, but specifying the contrast matrix by hand
> C <- c(0, 0, 0, -1, -1, -1, 1, 0, 0, 0, 1, 0, 0, 0, 1)
> C <- matrix(C, ncol=5)
> # numerate the contrasts
> rownames(C) <- paste("C", 1:nrow(C), sep="")
> test <- simint(minutes~blanket, data=recovery, conf.level=0.9, 
+                alternative="less",eps=0.0001, cmatrix=C[,-1])
Warning message:
'simint.default' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> print(test)

	 Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: lm(formula = minutes ~ blanket, data = recovery)

Estimated Quantile = -1.8431

Linear Hypotheses:
        Estimate lwr     upr    
C1 >= 0 -2.1333     -Inf  0.8227
C2 >= 0 -7.4667     -Inf -4.5107
C3 >= 0 -1.6667     -Inf -0.0359

90% family-wise confidence level
 

> 
> # same results, but more detailed information using the summary method
> summary(test)

	 Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: lm(formula = minutes ~ blanket, data = recovery)

Estimated Quantile = -1.8431

Linear Hypotheses:
        Estimate lwr     upr    
C1 >= 0 -2.1333     -Inf  0.8227
C2 >= 0 -7.4667     -Inf -4.5107
C3 >= 0 -1.6667     -Inf -0.0359

90% family-wise confidence level
 

> 
> ## Keywords: 'datasets'.
> 
> 
> cleanEx(); ..nameEx <- "simint"
> ###--- >>> `simint' <<<----- Simultaneous Intervals
> 
> 	## alias	 help(simint)
> 	## alias	 help(simint.default)
> 	## alias	 help(simint.formula)
> 
> ##___ Examples ___:
> 
> data(recovery)
> 
> # one-sided simultaneous confidence intervals for Dunnett 
> # in the one-way layout
> summary(simint(minutes~blanket, data=recovery, type="Dunnett", conf.level=0.9, 
+        alternative="less",eps=0.0001))

	 Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = minutes ~ blanket, data = recovery)

Estimated Quantile = -1.843

Linear Hypotheses:
             Estimate lwr      upr     
b1 - b0 >= 0 -2.13333     -Inf  0.82239
b2 - b0 >= 0 -7.46667     -Inf -4.51094
b3 - b0 >= 0 -1.66667     -Inf -0.03606

90% family-wise confidence level
 

Warning message:
'simint.default' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> 
> 
> ## Keywords: 'htest'.
> 
> 
> cleanEx(); ..nameEx <- "simtest"
> ###--- >>> `simtest' <<<----- Simultaneous comparisons
> 
> 	## alias	 help(simtest.default)
> 	## alias	 help(simtest.formula)
> 	## alias	 help(simtest)
> 
> ##___ Examples ___:
> 
> data(cholesterol)
> 
> # adjusted p-values for all-pairwise comparisons in a onw-way 
> # layout (tests for restricted combinations)
> simtest(response ~ trt, data=cholesterol, type="Tukey", ttype="logical")

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = response ~ trt, data = cholesterol)

Linear Hypotheses:
                     Estimate Std. Error t value p value    
2times - 1time == 0     3.443      1.443   2.385 0.13816    
4times - 1time == 0     6.593      1.443   4.568 < 0.001 ***
drugD - 1time == 0      9.579      1.443   6.637 < 0.001 ***
drugE - 1time == 0     15.166      1.443  10.507 < 0.001 ***
4times - 2times == 0    3.150      1.443   2.182 0.20503    
drugD - 2times == 0     6.136      1.443   4.251 < 0.001 ***
drugE - 2times == 0    11.723      1.443   8.122 < 0.001 ***
drugD - 4times == 0     2.986      1.443   2.069 0.25124    
drugE - 4times == 0     8.573      1.443   5.939 < 0.001 ***
drugE - drugD == 0      5.586      1.443   3.870 0.00307 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported)

Warning message:
'simtest.default' is deprecated.
Use 'glht' instead.
See help("Deprecated") and help("multcomp-deprecated"). 
> 
> 
> ## Keywords: 'htest'.
> 
> 
> cleanEx(); ..nameEx <- "tire"
> ###--- >>> `tire' <<<----- Tire Wear Data Set
> 
> 	## alias	 help(tire)
> 
> ##___ Examples ___:
> 
> #tire <- read.csv("tire.csv", header = TRUE)
> #C <- c(0,1,-1,0,10,-10)
> #for ( x in seq(15,70,5) ) { C <- rbind( C,c(0,1,-1,0,x,-x) ) }
> ## numerate the contrasts
> #rownames(C) <- paste("C", 1:nrow(C), sep="")
> #
> ## simultaneous confidence intervals of two regression functions
> #summary(simint(cost ~ make + mph + make:mph, data=tire,
> #               cmatrix=C, eps=0.001, whichf = NULL))
> 
> ## Keywords: 'datasets'.
> 
> 
> cleanEx(); ..nameEx <- "waste"
> ###--- >>> `waste' <<<----- Industrial Waste Data Set
> 
> 	## alias	 help(waste)
> 
> ##___ Examples ___:
> 
> data(waste)
> summary(aov(waste ~ envir + temp + envir*temp, data=waste))
            Df  Sum Sq Mean Sq F value    Pr(>F)    
envir        4 24.6854  6.1713  5.2532 0.0075463 ** 
temp         2 30.6928 15.3464 13.0632 0.0005185 ***
envir:temp   8 22.9116  2.8639  2.4378 0.0651340 .  
Residuals   15 17.6217  1.1748                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> #summary(simint(waste ~ envir:temp, data=waste,
> #               type="Tetrade", eps = 0.01))
> 
> ## Keywords: 'datasets'.
> 
> 
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  11.652 0.388 12.217 0 0 
> dev.off(); quit('no')
null device 
          1