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("gmodels-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('gmodels')
Loading required package: MASS
Loading required package: gdata
Loading required package: gtools
Loading required package: gplots

Attaching package: 'gplots'


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

	 lowess 

> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "CrossTable"
> 
> ### * CrossTable
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: CrossTable
> ### Title: Cross Tabulation with Tests for Factor Independence
> ### Aliases: CrossTable
> ### Keywords: category univar
> 
> ### ** Examples
> 
> 
> # Simple cross tabulation of education versus prior induced abortions
> # using infertility data
> data(infert, package = "datasets")
> CrossTable(infert$education, infert$induced, expected = TRUE)
Warning in chisq.test(t, correct = FALSE, ...) : 
	 Chi-squared approximation may be incorrect

 
   Cell Contents
|-------------------------|
|                       N |
|              Expected N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  248 

 
                 | infert$induced 
infert$education |         0 |         1 |         2 | Row Total | 
-----------------|-----------|-----------|-----------|-----------|
          0-5yrs |         4 |         2 |         6 |        12 | 
                 |     6.919 |     3.290 |     1.790 |           | 
                 |     1.232 |     0.506 |     9.898 |           | 
                 |     0.333 |     0.167 |     0.500 |     0.048 | 
                 |     0.028 |     0.029 |     0.162 |           | 
                 |     0.016 |     0.008 |     0.024 |           | 
-----------------|-----------|-----------|-----------|-----------|
         6-11yrs |        78 |        27 |        15 |       120 | 
                 |    69.194 |    32.903 |    17.903 |           | 
                 |     1.121 |     1.059 |     0.471 |           | 
                 |     0.650 |     0.225 |     0.125 |     0.484 | 
                 |     0.545 |     0.397 |     0.405 |           | 
                 |     0.315 |     0.109 |     0.060 |           | 
-----------------|-----------|-----------|-----------|-----------|
         12+ yrs |        61 |        39 |        16 |       116 | 
                 |    66.887 |    31.806 |    17.306 |           | 
                 |     0.518 |     1.627 |     0.099 |           | 
                 |     0.526 |     0.336 |     0.138 |     0.468 | 
                 |     0.427 |     0.574 |     0.432 |           | 
                 |     0.246 |     0.157 |     0.065 |           | 
-----------------|-----------|-----------|-----------|-----------|
    Column Total |       143 |        68 |        37 |       248 | 
                 |     0.577 |     0.274 |     0.149 |           | 
-----------------|-----------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  16.53059     d.f. =  4     p =  0.002383898 


 
> CrossTable(infert$education, infert$induced, expected = TRUE, format="SAS")
Warning in chisq.test(t, correct = FALSE, ...) : 
	 Chi-squared approximation may be incorrect

 
   Cell Contents
|-------------------------|
|                       N |
|              Expected N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  248 

 
                 | infert$induced 
infert$education |         0 |         1 |         2 | Row Total | 
-----------------|-----------|-----------|-----------|-----------|
          0-5yrs |         4 |         2 |         6 |        12 | 
                 |     6.919 |     3.290 |     1.790 |           | 
                 |     1.232 |     0.506 |     9.898 |           | 
                 |     0.333 |     0.167 |     0.500 |     0.048 | 
                 |     0.028 |     0.029 |     0.162 |           | 
                 |     0.016 |     0.008 |     0.024 |           | 
-----------------|-----------|-----------|-----------|-----------|
         6-11yrs |        78 |        27 |        15 |       120 | 
                 |    69.194 |    32.903 |    17.903 |           | 
                 |     1.121 |     1.059 |     0.471 |           | 
                 |     0.650 |     0.225 |     0.125 |     0.484 | 
                 |     0.545 |     0.397 |     0.405 |           | 
                 |     0.315 |     0.109 |     0.060 |           | 
-----------------|-----------|-----------|-----------|-----------|
         12+ yrs |        61 |        39 |        16 |       116 | 
                 |    66.887 |    31.806 |    17.306 |           | 
                 |     0.518 |     1.627 |     0.099 |           | 
                 |     0.526 |     0.336 |     0.138 |     0.468 | 
                 |     0.427 |     0.574 |     0.432 |           | 
                 |     0.246 |     0.157 |     0.065 |           | 
-----------------|-----------|-----------|-----------|-----------|
    Column Total |       143 |        68 |        37 |       248 | 
                 |     0.577 |     0.274 |     0.149 |           | 
-----------------|-----------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  16.53059     d.f. =  4     p =  0.002383898 


 
> CrossTable(infert$education, infert$induced, expected = TRUE, format="SPSS")
Warning in chisq.test(t, correct = FALSE, ...) : 
	 Chi-squared approximation may be incorrect

   Cell Contents
|-------------------------|
|                   Count |
|         Expected Values |
| Chi-square contribution |
|             Row Percent |
|          Column Percent |
|           Total Percent |
|-------------------------|

Total Observations in Table:  248 

                 | infert$induced 
infert$education |        0  |        1  |        2  | Row Total | 
-----------------|-----------|-----------|-----------|-----------|
          0-5yrs |        4  |        2  |        6  |       12  | 
                 |    6.919  |    3.290  |    1.790  |           | 
                 |    1.232  |    0.506  |    9.898  |           | 
                 |   33.333% |   16.667% |   50.000% |    4.839% | 
                 |    2.797% |    2.941% |   16.216% |           | 
                 |    1.613% |    0.806% |    2.419% |           | 
-----------------|-----------|-----------|-----------|-----------|
         6-11yrs |       78  |       27  |       15  |      120  | 
                 |   69.194  |   32.903  |   17.903  |           | 
                 |    1.121  |    1.059  |    0.471  |           | 
                 |   65.000% |   22.500% |   12.500% |   48.387% | 
                 |   54.545% |   39.706% |   40.541% |           | 
                 |   31.452% |   10.887% |    6.048% |           | 
-----------------|-----------|-----------|-----------|-----------|
         12+ yrs |       61  |       39  |       16  |      116  | 
                 |   66.887  |   31.806  |   17.306  |           | 
                 |    0.518  |    1.627  |    0.099  |           | 
                 |   52.586% |   33.621% |   13.793% |   46.774% | 
                 |   42.657% |   57.353% |   43.243% |           | 
                 |   24.597% |   15.726% |    6.452% |           | 
-----------------|-----------|-----------|-----------|-----------|
    Column Total |      143  |       68  |       37  |      248  | 
                 |   57.661% |   27.419% |   14.919% |           | 
-----------------|-----------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  16.53059     d.f. =  4     p =  0.002383898 


 
       Minimum expected frequency: 1.790323 
Cells with Expected Frequency < 5: 2 of 9 (22.22222%)

> 
> 
> 
> cleanEx(); ..nameEx <- "ci"
> 
> ### * ci
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ci
> ### Title: Compute Confidence Intervals
> ### Aliases: ci ci.default ci.lm ci.lme
> ### Keywords: regression
> 
> ### ** Examples
> 
> 
> data(state)
> reg  <-  lm(Area ~ Population, data=as.data.frame(state.x77))
> ci(reg) 
                Estimate     CI lower     CI upper   Std. Error      p-value
(Intercept) 6.890624e+04 34919.405723 102893.06460 16903.532236 0.0001710224
Population  4.308676e-01    -5.114345      5.97608     2.757941 0.8765085381
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "coefFrame"
> 
> ### * coefFrame
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: coefFrame
> ### Title: Return model parameters in a data frame
> ### Aliases: coefFrame
> ### Keywords: models
> 
> ### ** Examples
> 
> # load example data
> library(gtools)
> data(ELISA)
> 
> # Coefficients for four parameter logistic fits:
> coefFrame(log(Signal) ~ SSfpl(log(Concentration), A, B, xmid, scal),
+            data = ELISA, fitfun = nls,
+            by = c("PlateDay", "Read"),
+            fit.on = Description == "Standard" & Concentration != 0)
Warning: no finite arguments to max; returning -Inf
Warning: no finite arguments to max; returning -Inf
Warning in "[<-.factor"(`*tmp*`, ri, value = c(1, 1, 1, 2, 2, 2, 3, 3, 3,  : 
	 invalid factor level, NAs generated
           PlateDay Read  A  B xmid scal
1   Plate 1 (Day 1)    1 NA NA   NA   NA
43  Plate 1 (Day 1)    2 NA NA   NA   NA
85  Plate 1 (Day 1)    3 NA NA   NA   NA
127 Plate 2 (Day 1)    1 NA NA   NA   NA
169 Plate 2 (Day 1)    2 NA NA   NA   NA
211 Plate 2 (Day 1)    3 NA NA   NA   NA
253 Plate 3 (Day 2)    1 NA NA   NA   NA
295 Plate 3 (Day 2)    2 NA NA   NA   NA
337 Plate 3 (Day 2)    3 NA NA   NA   NA
379 Plate 4 (Day 2)    1 NA NA   NA   NA
421 Plate 4 (Day 2)    2 NA NA   NA   NA
463 Plate 4 (Day 2)    3 NA NA   NA   NA
> 
> # Coefficients for linear fits:
> coefFrame(log(Signal) ~ log(Concentration), 
+            data = ELISA, fitfun = lm, 
+            by = c("PlateDay", "Read"),
+            fit.on = Description == "Standard" & Concentration != 0 )
Warning: no finite arguments to max; returning -Inf
Warning: no finite arguments to max; returning -Inf
Warning in "[<-.factor"(`*tmp*`, ri, value = c(1, 1, 1, 2, 2, 2, 3, 3, 3,  : 
	 invalid factor level, NAs generated
           PlateDay Read X.Intercept. log.Concentration.
1   Plate 1 (Day 1)    1           NA                 NA
43  Plate 1 (Day 1)    2           NA                 NA
85  Plate 1 (Day 1)    3           NA                 NA
127 Plate 2 (Day 1)    1           NA                 NA
169 Plate 2 (Day 1)    2           NA                 NA
211 Plate 2 (Day 1)    3           NA                 NA
253 Plate 3 (Day 2)    1           NA                 NA
295 Plate 3 (Day 2)    2           NA                 NA
337 Plate 3 (Day 2)    3           NA                 NA
379 Plate 4 (Day 2)    1           NA                 NA
421 Plate 4 (Day 2)    2           NA                 NA
463 Plate 4 (Day 2)    3           NA                 NA
> 
> # Example passing arguments to fitfun, and example of
> # error handling during model fitting:
> ELISA$Signal[1] <- NA
> coefFrame(log(Signal) ~ log(Concentration), 
+            data = ELISA, fitfun = lm, na.action = na.fail,
+            by = c("PlateDay", "Read"),
+            fit.on = Description == "Standard" & Concentration != 0 )
Warning: no finite arguments to max; returning -Inf
Warning: no finite arguments to max; returning -Inf
Warning in "[<-.factor"(`*tmp*`, ri, value = c(1, 1, 1, 2, 2, 2, 3, 3, 3,  : 
	 invalid factor level, NAs generated
           PlateDay Read X.Intercept. log.Concentration.
1   Plate 1 (Day 1)    1           NA                 NA
43  Plate 1 (Day 1)    2           NA                 NA
85  Plate 1 (Day 1)    3           NA                 NA
127 Plate 2 (Day 1)    1           NA                 NA
169 Plate 2 (Day 1)    2           NA                 NA
211 Plate 2 (Day 1)    3           NA                 NA
253 Plate 3 (Day 2)    1           NA                 NA
295 Plate 3 (Day 2)    2           NA                 NA
337 Plate 3 (Day 2)    3           NA                 NA
379 Plate 4 (Day 2)    1           NA                 NA
421 Plate 4 (Day 2)    2           NA                 NA
463 Plate 4 (Day 2)    3           NA                 NA
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "estimable"
> 
> ### * estimable
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: estimable
> ### Title: Contrasts and estimable linear functions of model coefficients
> ###   for lm, glm, lme, and geese objects
> ### Aliases: estimable .wald
> ### Keywords: models regression
> 
> ### ** Examples
> 
> # simple contrast and confidence interval
> y <- rnorm(100)
> x <-  cut(rnorm(100, mean=y, sd=0.25),c(-4,-1.5,0,1.5,4))
> reg <- lm(y ~ x)
> estimable(reg, c(    0,   1,    0,   -1) )
            Estimate Std. Error   t value DF Pr(>|t|)
(0 1 0 -1) -2.631369  0.2298629 -11.44756 96        0
> 
> # Fit a spline with a single knot at 0.5 and plot the *pointwise*
> # confidence intervals
> library(gplots)
> x2 <- rnorm(100,mean=y,sd=0.5)
> reg2 <- lm(y ~ x + x2 + pmax(x2-0.5,0) )
> range <- seq(-2,2,,50)
> tmp <- estimable(reg2,cbind(1,0,0,1,range,pmax(range-0.5,0)), conf.int=0.95)
> plotCI(x=range,y=tmp[,1],li=tmp[,6],ui=tmp[,7])
> 
> # Fit both linear and quasi-Poisson models to iris data, then compute
> # conficence intervals on a contrast.
> data(iris)
> lm1  <- lm(Sepal.Length~Sepal.Width+Species+Sepal.Width:Species, data=iris)
> glm1 <- glm(Sepal.Length~Sepal.Width+Species+Sepal.Width:Species, data=iris,
+             family=quasipoisson("identity"))
> cm <- rbind( lambda1 = c(1,0,1,0,0,0),
+              lambda2 = c(1,0,0,1,0,0))
> estimable(lm1,cm)
        Estimate Std. Error  t value  DF     Pr(>|t|)
lambda1 3.539735  0.5580189 6.343396 144 2.735194e-09
lambda2 3.906836  0.5826525 6.705260 144 4.247154e-10
> estimable(glm1,cm)
        Estimate Std. Error  t value  DF     Pr(>|t|)
lambda1 3.549904  0.5422795 6.546263 144 9.688996e-10
lambda2 3.921773  0.6056852 6.474937 144 1.398190e-09
> 
> 
> 
> cleanEx(); ..nameEx <- "fast.prcomp"
> 
> ### * fast.prcomp
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fast.prcomp
> ### Title: Efficient computation of principal components and singular value
> ###   decompositions.
> ### Aliases: fast.prcomp fast.svd
> ### Keywords: multivariate algebra array
> 
> ### ** Examples
> 
> 
>   # create test matrix
>   set.seed(4943546)
>   nr <- 50
>   nc <- 2000
>   x  <- matrix( rnorm( nr*nc), nrow=nr, ncol=nc )
>   tx <- t(x)
> 
>   # SVD directly on matrix is SLOW:
>   system.time( val.x <- svd(x)$u )
[1] 0.12 0.06 0.18 0.00 0.00
> 
>   # SVD on t(matrix) is FAST:
>   system.time( val.tx <- svd(tx)$v )
[1] 0.11 0.03 0.14 0.00 0.00
> 
>   # and the results are equivalent:
>   max( abs(val.x) - abs(val.tx) )
[1] 1.560002e-13
> 
>   # Time gap dissapears using fast.svd:
>   system.time( val.x <- fast.svd(x)$u )
[1] 0.31 0.02 0.33 0.00 0.00
>   system.time( val.tx <- fast.svd(tx)$v )
[1] 0.11 0.03 0.14 0.00 0.00
>   max( abs(val.x) - abs(val.tx) )
[1] 1.560002e-13
> 
>   library(stats)
> 
>   # prcomp directly on matrix is SLOW:
>   system.time( pr.x <- prcomp(x) )
[1] 0.21 0.06 0.27 0.00 0.00
> 
>   # prcomp.fast is much faster
>   system.time( fast.pr.x <- fast.prcomp(x) )
[1] 0.30 0.05 0.36 0.00 0.00
> 
>   # and the results are equivalent
>   max( pr.x$sdev - fast.pr.x$sdev )
[1] 0
>   max( abs(pr.x$rotation[,1:49]) - abs(fast.pr.x$rotation[,1:49]) )
[1] 0
>   max( abs(pr.x$x) - abs(fast.pr.x$x)  )
[1] 0
> 
>   # (except for the last and least significant component):
>   max( abs(pr.x$rotation[,50]) - abs(fast.pr.x$rotation[,50]) )
[1] 0
> 
> 
> 
> cleanEx(); ..nameEx <- "fit.contrast"
> 
> ### * fit.contrast
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fit.contrast
> ### Title: Compute and test arbitrary contrasts for regression objects
> ### Aliases: fit.contrast fit.contrast.lm fit.contrast.lme
> ### Keywords: models regression
> 
> ### ** Examples
> 
> y <- rnorm(100)
> x <-  cut(rnorm(100, mean=y, sd=0.25),c(-4,-1.5,0,1.5,4))
> reg <- lm(y ~ x)
> summary(reg)

Call:
lm(formula = y ~ x)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.891668 -0.254190 -0.003064  0.290988  1.000554 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.8831     0.2184  -8.624 1.34e-13 ***
x(-1.5,0]     1.2893     0.2299   5.609 1.96e-07 ***
x(0,1.5]      2.4694     0.2262  10.919  < 2e-16 ***
x(1.5,4]      3.9206     0.3088  12.696  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.4367 on 96 degrees of freedom
Multiple R-Squared: 0.7708,	Adjusted R-squared: 0.7636 
F-statistic: 107.6 on 3 and 96 DF,  p-value: < 2.2e-16 

> 
> # look at the group means
> gm <- sapply(split(y,x),mean)
> gm
 (-4,-1.5]   (-1.5,0]    (0,1.5]    (1.5,4] 
-1.8831443 -0.5938913  0.5862796  2.0374775 
> 
> # mean of 1st group vs mean of 4th group
> fit.contrast(reg, x, c(    1,    0,    0,   -1) )
                  Estimate Std. Error   t value     Pr(>|t|)
x c=( 1 0 0 -1 ) -3.920622  0.3088111 -12.69586 2.922641e-22
> # estimate should be equal to:
> gm[1] - gm[4]
(-4,-1.5] 
-3.920622 
> 
> # mean of 1st and 2nd groups vs mean of 3rd and 4th groups
> fit.contrast(reg, x, c( -1/2, -1/2,  1/2,  1/2) )
                          Estimate Std. Error  t value     Pr(>|t|)
x c=( -0.5 -0.5 0.5 0.5 ) 2.550396   0.161235 15.81788 1.735806e-28
> # estimate should be equal to:
> sum(-1/2*gm[1], -1/2*gm[2], 1/2*gm[3], 1/2*gm[4])
[1] 2.550396
> 
> # mean of 1st group vs mean of 2nd, 3rd and 4th groups
> fit.contrast(reg, x, c( -3/3,  1/3,  1/3,  1/3) ) 
                                                                 Estimate
x c=( -1 0.333333333333333 0.333333333333333 0.333333333333333 ) 2.559766
                                                                 Std. Error
x c=( -1 0.333333333333333 0.333333333333333 0.333333333333333 )  0.2322460
                                                                  t value
x c=( -1 0.333333333333333 0.333333333333333 0.333333333333333 ) 11.02179
                                                                     Pr(>|t|)
x c=( -1 0.333333333333333 0.333333333333333 0.333333333333333 ) 9.674065e-19
> # estimate should be equal to:
> sum(-3/3*gm[1], 1/3*gm[2], 1/3*gm[3], 1/3*gm[4])
[1] 2.559766
> 
> # all at once
> cmat <- rbind( "1 vs 4"    =c(-1, 0, 0, 1),
+                "1+2 vs 3+4"=c(-1/2,-1/2, 1/2, 1/2),
+                "1 vs 2+3+4"=c(-3/3, 1/3, 1/3, 1/3))
> fit.contrast(reg,x,cmat)
            Estimate Std. Error  t value     Pr(>|t|)
x1 vs 4     3.920622  0.3088111 12.69586 2.922641e-22
x1+2 vs 3+4 2.550396  0.1612350 15.81788 1.735806e-28
x1 vs 2+3+4 2.559766  0.2322460 11.02179 9.674065e-19
> 
> #
> x2 <- rnorm(100,mean=y,sd=0.5)
> reg2 <- lm(y ~ x + x2 )
> fit.contrast(reg2,x,c(-1,0,0,1))
                 Estimate Std. Error t value     Pr(>|t|)
x c=( -1 0 0 1 ) 2.134428   0.342401 6.23371 1.248135e-08
> 
> #
> # Example for Analysis of Variance
> #
> 
> set.seed(03215)
> Genotype <- sample(c("WT","KO"), 1000, replace=TRUE)
> Time <- factor(sample(1:3, 1000, replace=TRUE))
> y <- rnorm(1000)
> data <- data.frame(y, Genotype, Time)
> 
> # Compute Contrasts & obtain 95% confidence intervals
> 
> model <- aov( y ~ Genotype + Time + Genotype:Time, data=data )
> 
> fit.contrast( model, "Genotype", rbind("KO vs WT"=c(-1,1) ), conf=0.95 )
                   Estimate Std. Error   t value Pr(>|t|)   lower CI  upper CI
GenotypeKO vs WT 0.01683876  0.1095764 0.1536714   0.8779 -0.1981888 0.2318664
> 
> fit.contrast( model, "Time",
+             rbind("1 vs 2"=c(-1,1,0),
+                   "2 vs 3"=c(0,-1,1)
+                   ),
+             conf=0.95 )
               Estimate Std. Error     t value  Pr(>|t|)   lower CI   upper CI
Time1 vs 2 -0.003321336  0.1150931 -0.02885782 0.9769838 -0.2291746 0.22253196
Time2 vs 3 -0.131231895  0.1152710 -1.13846378 0.2552012 -0.3574344 0.09497061
> 
> cm.G <- rbind("KO vs WT"=c(-1,1) )
> cm.T <- rbind("1 vs 2"=c(-1,1,0),
+               "2 vs 3"=c(0,-1,1) )
> 
> # Compute contrasts and show SSQ decompositions
> 
> model <- aov( y ~ Genotype + Time + Genotype:Time, data=data,
+               contrasts=list(Genotype=make.contrasts(cm.G),
+                              Time=make.contrasts(cm.T) )
+             )
> 
> summary(model, split=list( Genotype=list( "KO vs WT"=1 ),
+                            Time = list( "1 vs 2" = 1,
+                                         "2 vs 3" = 2 ) ) )
                                  Df  Sum Sq Mean Sq F value Pr(>F)
Genotype                           1    1.17    1.17  1.1207 0.2900
  Genotype: KO vs WT               1    1.17    1.17  1.1207 0.2900
Time                               2    0.74    0.37  0.3526 0.7029
  Time: 1 vs 2                     1    0.18    0.18  0.1710 0.6793
  Time: 2 vs 3                     1    0.56    0.56  0.5342 0.4650
Genotype:Time                      2    1.15    0.58  0.5524 0.5758
  Genotype:Time: KO vs WT.1 vs 2   1    0.33    0.33  0.3131 0.5759
  Genotype:Time: KO vs WT.2 vs 3   1    0.83    0.83  0.7916 0.3738
Residuals                        994 1036.64    1.04               
> 
> # example for lme 
> library(nlme)
> data(Orthodont)
> fm1 <- lme(distance ~ Sex, data = Orthodont,random=~1|Subject) 
> 
> # Contrast for sex.  This example is equivalent to standard treatment
> # contrast.
> #
> fit.contrast(fm1, "Sex", c(-1,1), conf.int=0.95 )
                Estimate Std. Error   t-value    Pr(>|t|)  lower CI   upper CI
Sex c=( -1 1 ) -2.321023  0.7614192 -3.048285 0.005375176 -3.889195 -0.7528506
> #
> # and identical results can be obtained using lme built-in 'intervals'
> #
> intervals(fm1)
Approximate 95% confidence intervals

 Fixed effects:
                lower      est.      upper
(Intercept) 24.001758 24.968750 25.9357420
SexFemale   -3.889195 -2.321023 -0.7528506
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: Subject 
                   lower     est.    upper
sd((Intercept)) 1.050702 1.595847 2.423834

 Within-group standard error:
   lower     est.    upper 
1.903425 2.220310 2.589949 
> 
> # Cut age into quantile groups & compute some contrasts
> Orthodont$AgeGroup <- gtools::quantcut(Orthodont$age)
> fm2 <- lme(distance ~ Sex + AgeGroup, data = Orthodont,random=~1|Subject)
> #
> fit.contrast(fm2, "AgeGroup", rbind("Linear"=c(-2,-1,1,2),
+                                     "U-Shaped"=c(-1,1,1,-1),
+                                     "Change-Point at 11"=c(-1,-1,1,1)),
+                               conf.int=0.95)
                            Estimate Std. Error    t-value     Pr(>|t|)
AgeGroupLinear              9.296296  0.8773838 10.5954731 0.000000e+00
AgeGroupU-Shaped           -0.462963  0.5549062 -0.8343085 4.066551e-01
AgeGroupChange-Point at 11  5.388889  0.5549062  9.7113505 4.440892e-15
                            lower CI   upper CI
AgeGroupLinear              7.549559 11.0430332
AgeGroupU-Shaped           -1.567696  0.6417705
AgeGroupChange-Point at 11  4.284155  6.4936223
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "glh.test"
> 
> ### * glh.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: glh.test
> ### Title: Test a General Linear Hypothesis for a Regression Model
> ### Aliases: glh.test print.glh.test summary.glh.test
> ### Keywords: models regression
> 
> ### ** Examples
> 
> 
> # fit a simple model
> y <- rnorm(100)
> x <-  cut(rnorm(100, mean=y, sd=0.25),c(-4,-1.5,0,1.5,4))
> reg <- lm(y ~ x)
> summary(reg)

Call:
lm(formula = y ~ x)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.891668 -0.254190 -0.003064  0.290988  1.000554 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.8831     0.2184  -8.624 1.34e-13 ***
x(-1.5,0]     1.2893     0.2299   5.609 1.96e-07 ***
x(0,1.5]      2.4694     0.2262  10.919  < 2e-16 ***
x(1.5,4]      3.9206     0.3088  12.696  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.4367 on 96 degrees of freedom
Multiple R-Squared: 0.7708,	Adjusted R-squared: 0.7636 
F-statistic: 107.6 on 3 and 96 DF,  p-value: < 2.2e-16 

> 
> # test both group 1 = group 2  and group 3 = group 4
> # *Note the 0 in the column for the intercept term*
> 
> C <- rbind( c(0,-1,0,0), c(0,0,-1,1) )
> ret <- glh.test(reg, C)
> ret  # same as 'print(ret) '

	 Test of General Linear Hypothesis 	
Call:
glh.test(reg = reg, cm = C)
F = 36.3156, df1 =  2, df2 = 96, p-value = 1.803e-12 

> summary(ret)

	 Test of General Linear Hypothesis 	
Regression:  reg 

Null Hypothesis: C %*% Beta-hat = d 

C matrix: 
     (Intercept) x(-1.5,0] x(0,1.5] x(1.5,4]
[1,]           0        -1        0        0
[2,]           0         0       -1        1

d vector: 
[1] 0 0

C %*% Beta-hat: 
[1] -1.289253  1.451198

F = 36.3156, df1 =  2, df2 = 96, p-value = 1.803e-12 

> 
> # To compute a contrast between the first and second level of the factor
> # 'x' using 'fit.contrast' gives:
> 
> fit.contrast( reg, x,c(1,-1,0,0) )
                  Estimate Std. Error  t value     Pr(>|t|)
x c=( 1 -1 0 0 ) -1.289253  0.2298629 -5.60879 1.959651e-07
>         
> # To test this same contrast using 'glh.test', use a contrast matrix
> # with a zero coefficient for the intercept term.  See the Note section,
> # above, for an explanation.
> 
> C <- rbind( c(0,-1,0,0) )
> glh.test( reg, C )

	 Test of General Linear Hypothesis 	
Call:
glh.test(reg = reg, cm = C)
F = 31.4585, df1 =  1, df2 = 96, p-value = 1.960e-07 

> 
> 
> 
> 
> cleanEx(); ..nameEx <- "make.contrasts"
> 
> ### * make.contrasts
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: make.contrasts
> ### Title: Construct a User-Specified Contrast Matrix
> ### Aliases: make.contrasts
> ### Keywords: models regression
> 
> ### ** Examples
> 
> set.seed(4684)
> y <- rnorm(100)
> x.true <- rnorm(100, mean=y, sd=0.25)
> x <-  factor(cut(x.true,c(-4,-1.5,0,1.5,4)))
> reg <- lm(y ~ x)
> summary(reg)

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.24953 -0.32750  0.01022  0.39324  0.79245 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.7287     0.1907  -9.066 1.51e-14 ***
x(-1.5,0]     1.1677     0.2027   5.760 1.01e-07 ***
x(0,1.5]      2.2522     0.2035  11.065  < 2e-16 ***
x(1.5,4]      3.5063     0.2828  12.397  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.4671 on 96 degrees of freedom
Multiple R-Squared: 0.7405,	Adjusted R-squared: 0.7324 
F-statistic: 91.31 on 3 and 96 DF,  p-value: < 2.2e-16 

> 
> # Mirror default treatment contrasts
> test <- make.contrasts(rbind( c(-1,1,0,0), c(-1,0,1,0), c(-1,0,0,1) ))
> lm( y ~ x, contrasts=list(x = test ))

Call:
lm(formula = y ~ x, contrasts = list(x = test))

Coefficients:
(Intercept)          xC1          xC2          xC3  
   0.002878     1.167728     2.252186     3.506275  

> 
> # Specify some more complicated contrasts
> #   - mean of 1st group vs mean of 4th group
> #   - mean of 1st and 2nd groups vs mean of 3rd and 4th groups
> #   - mean of 1st group vs mean of 2nd, 3rd and 4th groups
> cmat <- rbind( "1 vs 4"    =c(-1, 0, 0, 1),
+                "1+2 vs 3+4"=c(-1/2,-1/2, 1/2, 1/2),
+                "1 vs 2+3+4"=c(-3/3, 1/3, 1/3, 1/3))
> 
> summary(lm( y ~ x, contrasts=list(x=make.contrasts(cmat) )))

Call:
lm(formula = y ~ x, contrasts = list(x = make.contrasts(cmat)))

Residuals:
     Min       1Q   Median       3Q      Max 
-1.24953 -0.32750  0.01022  0.39324  0.79245 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.002878   0.074919   0.038     0.97    
x1 vs 4     3.506275   0.282824  12.397   <2e-16 ***
x1+2 vs 3+4 2.295366   0.149838  15.319   <2e-16 ***
x1 vs 2+3+4 2.308729   0.205663  11.226   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.4671 on 96 degrees of freedom
Multiple R-Squared: 0.7405,	Adjusted R-squared: 0.7324 
F-statistic: 91.31 on 3 and 96 DF,  p-value: < 2.2e-16 

> # or
> contrasts(x) <- make.contrasts(cmat)
> summary(lm( y ~ x ) )

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.24953 -0.32750  0.01022  0.39324  0.79245 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.002878   0.074919   0.038     0.97    
x1 vs 4     3.506275   0.282824  12.397   <2e-16 ***
x1+2 vs 3+4 2.295366   0.149838  15.319   <2e-16 ***
x1 vs 2+3+4 2.308729   0.205663  11.226   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.4671 on 96 degrees of freedom
Multiple R-Squared: 0.7405,	Adjusted R-squared: 0.7324 
F-statistic: 91.31 on 3 and 96 DF,  p-value: < 2.2e-16 

> 
> # or use contrasts.lm
> reg <- lm(y ~ x)
> fit.contrast( reg, "x", cmat )
            Estimate Std. Error  t value     Pr(>|t|)
x1 vs 4     3.506275  0.2828243 12.39736 1.218181e-21
x1+2 vs 3+4 2.295366  0.1498378 15.31901 1.585924e-27
x1 vs 2+3+4 2.308729  0.2056633 11.22577 3.562316e-19
> 
> # compare with values computed directly using group means
> gm <- sapply(split(y,x),mean)
> gm 
 (-4,-1.5]   (-1.5,0]    (0,1.5]    (1.5,4] 
-1.7286691 -0.5609414  0.5235166  1.7776055 
> 
> #
> # Example for Analysis of Variance
> #
> 
> set.seed(03215)
> Genotype <- sample(c("WT","KO"), 1000, replace=TRUE)
> Time <- factor(sample(1:3, 1000, replace=TRUE))
> data <- data.frame(y, Genotype, Time)
> y <- rnorm(1000)
> 
> data <- data.frame(y, Genotype, as.factor(Time))
> 
> # Compute Contrasts & obtain 95% confidence intervals
> 
> model <- aov( y ~ Genotype + Time + Genotype:Time, data=data )
> 
> fit.contrast( model, "Genotype", rbind("KO vs WT"=c(-1,1) ), conf=0.95 )
                   Estimate Std. Error   t value Pr(>|t|)   lower CI  upper CI
GenotypeKO vs WT 0.01683876  0.1095764 0.1536714   0.8779 -0.1981888 0.2318664
> 
> fit.contrast( model, "Time",
+             rbind("1 vs 2"=c(-1,1,0),
+                   "2 vs 3"=c(0,-1,1)
+                   ),
+             conf=0.95 )
               Estimate Std. Error     t value  Pr(>|t|)   lower CI   upper CI
Time1 vs 2 -0.003321336  0.1150931 -0.02885782 0.9769838 -0.2291746 0.22253196
Time2 vs 3 -0.131231895  0.1152710 -1.13846378 0.2552012 -0.3574344 0.09497061
> 
> cm.G <- rbind("KO vs WT"=c(-1,1) )
> cm.T <- rbind("1 vs 2"=c(-1,1,0),
+               "2 vs 3"=c(0,-1,1) )
> 
> # Compute contrasts and show SSQ decompositions
> 
> model <- model <- aov( y ~ Genotype + Time + Genotype:Time, data=data,
+                       contrasts=list(Genotype=make.contrasts(cm.G),
+                                      Time=make.contrasts(cm.T) )
+                       )
> 
> summary(model, split=list( Genotype=list( "KO vs WT"=1 ),
+                            Time = list( "1 vs 2" = 1,
+                                         "2 vs 3" = 2 ) ) )
                                  Df  Sum Sq Mean Sq F value Pr(>F)
Genotype                           1    1.17    1.17  1.1207 0.2900
  Genotype: KO vs WT               1    1.17    1.17  1.1207 0.2900
Time                               2    0.74    0.37  0.3526 0.7029
  Time: 1 vs 2                     1    0.18    0.18  0.1710 0.6793
  Time: 2 vs 3                     1    0.56    0.56  0.5342 0.4650
Genotype:Time                      2    1.15    0.58  0.5524 0.5758
  Genotype:Time: KO vs WT.1 vs 2   1    0.33    0.33  0.3131 0.5759
  Genotype:Time: KO vs WT.2 vs 3   1    0.83    0.83  0.7916 0.3738
Residuals                        994 1036.64    1.04               
> 
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  4.95 0.42 5.66 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')