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. > ### *
> ### > 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("multcomp-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('multcomp') Loading required package: mvtnorm > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "angina" > > ### * angina > > flush(stderr()); flush(stdout()) > > ### Name: angina > ### Title: Dose Response Data Set > ### Aliases: angina > ### Keywords: datasets > > ### ** Examples > > data(angina) > > # perform a dose-response analysis using simultaneous confidence > # intervals for Williams' contrasts > summary(ci <- simint(response~dose, data=angina, alternative="greater", + type="Williams")) Simultaneous 95% confidence intervals: Williams contrasts Call: simint.formula(formula = response ~ dose, data = angina, alternative = "greater", type = "Williams") Williams contrasts for factor dose Contrast matrix: dose0 dose1 dose2 dose3 dose4 C1 0 -1 0.00 0.0000000 0.0000000 1.0000000 C2 0 -1 0.00 0.0000000 0.5000000 0.5000000 C3 0 -1 0.00 0.3333333 0.3333333 0.3333333 C4 0 -1 0.25 0.2500000 0.2500000 0.2500000 Absolute Error Tolerance: 0.001 95 % quantile: 1.978 Coefficients: Estimate 5 % -- t value Std.Err. p raw p Bonf p adj C1 10.499 7.435 Inf 6.778 1.549 0 0 0 C2 7.747 5.093 Inf 5.775 1.341 0 0 0 C3 6.297 3.795 Inf 4.979 1.265 0 0 0 C4 5.247 2.824 Inf 4.284 1.225 0 0 0 > plot(ci, cex.axis=1.75, lwd=1.5) > > # consider the same data set as a changepoint problem > summary(simint(response~dose, data=angina, alternative="greater", + type="Changepoint")) Simultaneous 95% confidence intervals: Changepoint contrasts Call: simint.formula(formula = response ~ dose, data = angina, alternative = "greater", type = "Changepoint") Changepoint contrasts for factor dose Contrast matrix: dose0 dose1 dose2 dose3 dose4 C1 0 -1.0000000 0.2500000 0.2500000 0.2500000 0.2500000 C2 0 -0.5000000 -0.5000000 0.3333333 0.3333333 0.3333333 C3 0 -0.3333333 -0.3333333 -0.3333333 0.5000000 0.5000000 C4 0 -0.2500000 -0.2500000 -0.2500000 -0.2500000 1.0000000 Absolute Error Tolerance: 0.001 95 % quantile: 2.213 Coefficients: Estimate 5 % -- t value Std.Err. p raw p Bonf p adj C1 5.247 2.536 Inf 4.284 1.225 0 0 0 C2 5.250 3.037 Inf 5.250 1.000 0 0 0 C3 5.916 3.703 Inf 5.917 1.000 0 0 0 C4 7.877 5.167 Inf 6.433 1.225 0 0 0 > > # 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: McDermott contrasts Call: simtest.formula(formula = response ~ dose, data = angina, type = "McDermott", alternative = "greater", ttype = "logical") McDermott contrasts for factor dose Contrast matrix: dose0 dose1 dose2 dose3 dose4 C1 0 -1.0000000 1.0000000 0.0000000 0.00 0 C2 0 -0.5000000 -0.5000000 1.0000000 0.00 0 C3 0 -0.3333333 -0.3333333 -0.3333333 1.00 0 C4 0 -0.2500000 -0.2500000 -0.2500000 -0.25 1 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj C4 7.877 6.433 1.549 0.000 0.000 0.000 C3 3.164 2.502 1.341 0.008 0.024 0.024 C2 2.350 1.751 1.265 0.043 0.087 0.085 C1 2.095 1.353 1.225 0.091 0.091 0.091 > > > > cleanEx(); ..nameEx <- "cholesterol" > > ### * cholesterol > > flush(stderr()); flush(stdout()) > > ### 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: Tukey contrasts Call: simtest.formula(formula = response ~ trt, data = cholesterol, type = "Tukey", ttype = "logical") Contrast matrix: trt1time trt2times trt4times trtdrugD trtdrugE trt2times-trt1time 0 -1 1 0 0 0 trt4times-trt1time 0 -1 0 1 0 0 trtdrugD-trt1time 0 -1 0 0 1 0 trtdrugE-trt1time 0 -1 0 0 0 1 trt4times-trt2times 0 0 -1 1 0 0 trtdrugD-trt2times 0 0 -1 0 1 0 trtdrugE-trt2times 0 0 -1 0 0 1 trtdrugD-trt4times 0 0 0 -1 1 0 trtdrugE-trt4times 0 0 0 -1 0 1 trtdrugE-trtdrugD 0 0 0 0 -1 1 Adjusted P-Values p adj trtdrugE-trt1time 0.000 trtdrugE-trt2times 0.000 trtdrugD-trt1time 0.000 trtdrugE-trt4times 0.000 trt4times-trt1time 0.000 trtdrugD-trt2times 0.000 trtdrugE-trtdrugD 0.001 trt2times-trt1time 0.042 trt4times-trt2times 0.042 trtdrugD-trt4times 0.044 > > # 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: Tukey contrasts Call: simtest.formula(formula = response ~ trt, data = cholesterol, type = "Tukey", ttype = "free") Contrast matrix: trt1time trt2times trt4times trtdrugD trtdrugE trt2times-trt1time 0 -1 1 0 0 0 trt4times-trt1time 0 -1 0 1 0 0 trtdrugD-trt1time 0 -1 0 0 1 0 trtdrugE-trt1time 0 -1 0 0 0 1 trt4times-trt2times 0 0 -1 1 0 0 trtdrugD-trt2times 0 0 -1 0 1 0 trtdrugE-trt2times 0 0 -1 0 0 1 trtdrugD-trt4times 0 0 0 -1 1 0 trtdrugE-trt4times 0 0 0 -1 0 1 trtdrugE-trtdrugD 0 0 0 0 -1 1 Adjusted P-Values p adj trtdrugE-trt1time 0.000 trtdrugE-trt2times 0.000 trtdrugD-trt1time 0.000 trtdrugE-trt4times 0.000 trt4times-trt1time 0.000 trtdrugD-trt2times 0.001 trtdrugE-trtdrugD 0.002 trt2times-trt1time 0.057 trt4times-trt2times 0.063 trtdrugD-trt4times 0.063 > > # 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 <- mginv(t(x) %*% x) > 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: user-defined contrasts 95 % confidence intervals Estimate 2.5 % 97.5 % 2times-1time 3.443 -0.659 7.545 4times-1time 6.593 2.491 10.694 drugD-1time 9.579 5.478 13.681 drugE-1time 15.166 11.064 19.267 4times-2times 3.150 -0.952 7.251 drugD-2times 6.136 2.035 10.238 drugE-2times 11.723 7.621 15.824 drugD-4times 2.986 -1.115 7.088 drugE-4times 8.573 4.471 12.674 drugE-drugD 5.586 1.485 9.688 > csimtest(estpar=parm, df=df, covm=covm, cmatrix=cbind(0, contrast), + ttype="logical") Simultaneous tests: user-defined contrasts Contrast matrix: 1time 2times 4times drugD drugE 2times-1time 0 -1 1 0 0 0 4times-1time 0 -1 0 1 0 0 drugD-1time 0 -1 0 0 1 0 drugE-1time 0 -1 0 0 0 1 4times-2times 0 0 -1 1 0 0 drugD-2times 0 0 -1 0 1 0 drugE-2times 0 0 -1 0 0 1 drugD-4times 0 0 0 -1 1 0 drugE-4times 0 0 0 -1 0 1 drugE-drugD 0 0 0 0 -1 1 Adjusted P-Values p adj drugE-1time 0.000 drugE-2times 0.000 drugD-1time 0.000 drugE-4times 0.000 4times-1time 0.000 drugD-2times 0.000 drugE-drugD 0.001 2times-1time 0.042 4times-2times 0.042 drugD-4times 0.044 > > # only a subset of all pairwise hypotheses: > # > # * drug D versus all other formulations and > # * all pairwise comparisions for "1time", "2times" and "4times" > # > csubset = contrast[c(1,3,5,6,8,10),] > csubset 1time 2times 4times drugD drugE 2times-1time -1 1 0 0 0 drugD-1time -1 0 0 1 0 4times-2times 0 -1 1 0 0 drugD-2times 0 -1 0 1 0 drugD-4times 0 0 -1 1 0 drugE-drugD 0 0 0 -1 1 > simint(response ~ trt, data=cholesterol, cmatrix = csubset) Simultaneous confidence intervals: user-defined contrasts Call: simint.formula(formula = response ~ trt, data = cholesterol, cmatrix = csubset) 95 % confidence intervals Estimate 2.5 % 97.5 % 2times-1time 3.443 -0.422 7.308 drugD-1time 9.579 5.714 13.444 4times-2times 3.150 -0.715 7.015 drugD-2times 6.136 2.271 10.001 drugD-4times 2.986 -0.879 6.852 drugE-drugD 5.586 1.721 9.451 > > > > > cleanEx(); ..nameEx <- "contrMat" > > ### * contrMat > > flush(stderr()); flush(stdout()) > > ### Name: contrMat > ### Title: Contrast Matrices > ### Aliases: contr.Dunnett contrMat > ### Keywords: misc > > ### ** Examples > > n <- c(10,20,30,40) > a <- 2 > b <- 3 > names(n) <- paste("group", 1:4, sep="") > contrMat(n) # Dunnett is default group1 group2 group3 group4 group2-group1 -1 1 0 0 group3-group1 -1 0 1 0 group4-group1 -1 0 0 1 > contrMat(n, base=2) # use second level as baseline group1 group2 group3 group4 group1-group2 1 -1 0 0 group3-group2 0 -1 1 0 group4-group2 0 -1 0 1 > contr.Dunnett(names(n)) # it Moore-Penrose Inverse group2-group1 group3-group1 group4-group1 group1 -0.25 -0.25 -0.25 group2 0.75 -0.25 -0.25 group3 -0.25 0.75 -0.25 group4 -0.25 -0.25 0.75 > contrMat(n, type="Tukey") group1 group2 group3 group4 group2-group1 -1 1 0 0 group3-group1 -1 0 1 0 group4-group1 -1 0 0 1 group3-group2 0 -1 1 0 group4-group2 0 -1 0 1 group4-group3 0 0 -1 1 > contrMat(n, type="Sequen") group1 group2 group3 group4 group2-group1 -1 1 0 0 group3-group2 0 -1 1 0 group4-group3 0 0 -1 1 > contrMat(n, type="AVE") group1 group2 group3 group4 C 1 1.0000000 -0.2222222 -0.3333333 -0.4444444 C 2 -0.1250000 1.0000000 -0.3750000 -0.5000000 C 3 -0.1428571 -0.2857143 1.0000000 -0.5714286 C 4 -0.1666667 -0.3333333 -0.5000000 1.0000000 > contrMat(n, type="Changepoint") group1 group2 group3 group4 C1 -1.0000000 0.2222222 0.3333333 0.4444444 C2 -0.3333333 -0.6666667 0.4285714 0.5714286 C3 -0.1666667 -0.3333333 -0.5000000 1.0000000 > contrMat(n, type="Williams") group1 group2 group3 group4 C1 -1 0.0000000 0.0000000 1.0000000 C2 -1 0.0000000 0.4285714 0.5714286 C3 -1 0.2222222 0.3333333 0.4444444 > contrMat(n, type="Marcus") group1 group2 group3 group4 C1 -1.0000000 0.2222222 0.3333333 0.4444444 C2 -1.0000000 0.0000000 0.4285714 0.5714286 C3 -0.3333333 -0.6666667 0.4285714 0.5714286 C4 -1.0000000 0.0000000 0.0000000 1.0000000 C5 -0.3333333 -0.6666667 0.0000000 1.0000000 C6 -0.1666667 -0.3333333 -0.5000000 1.0000000 > contrMat(n, type="McDermott") group1 group2 group3 group4 C1 -1.0000000 1.0000000 0.0 0 C2 -0.3333333 -0.6666667 1.0 0 C3 -0.1666667 -0.3333333 -0.5 1 > contrMat(n, type="Tetrade", nlevel=c(a,b)) [,1] [,2] [,3] [,4] [,5] [,6] (11-12)-(21-22) 1 -1 0 -1 1 0 (11-13)-(21-23) 1 0 -1 -1 0 1 (12-13)-(22-23) 0 1 -1 0 -1 1 > > > > cleanEx(); ..nameEx <- "detergent" > > ### * detergent > > flush(stderr()); flush(stdout()) > > ### Name: detergent > ### Title: Detergent Durability Data Set > ### Aliases: detergent > ### Keywords: datasets > > ### ** Examples > > > 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") > > # adjusted p-values: either using the contrast matrix C > summary(simtest(plates ~ block+detergent, data=detergent, whichf="detergent", + ttype="logical", cmatrix=C)) Simultaneous tests: user-defined contrasts Call: simtest.formula(formula = plates ~ block + detergent, data = detergent, whichf = "detergent", ttype = "logical", cmatrix = C) user-defined contrasts for factor detergent, covariable: block Contrast matrix: 1 2 3 4 5 2-1 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 3-1 0 0 0 0 0 0 0 0 0 0 -1 0 1 0 0 4-1 0 0 0 0 0 0 0 0 0 0 -1 0 0 1 0 5-1 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 1 3-2 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 4-2 0 0 0 0 0 0 0 0 0 0 0 -1 0 1 0 5-2 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 1 4-3 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 5-3 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 1 5-4 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj 5-3 -7.933 -9.140 0.868 0.000 0.000 0.000 5-4 -6.533 -7.527 0.868 0.000 0.000 0.000 3-2 5.733 -6.606 0.868 0.000 0.000 0.000 5-1 -4.333 -4.993 0.868 0.000 0.001 0.000 4-2 4.333 -4.993 0.868 0.000 0.001 0.000 3-1 3.600 -4.148 0.868 0.001 0.003 0.003 5-2 -2.200 -2.535 0.868 0.022 0.044 0.043 4-1 2.200 -2.535 0.868 0.022 0.044 0.043 2-1 -2.133 -2.458 0.868 0.026 0.052 0.050 4-3 -1.400 -1.613 0.868 0.126 0.126 0.126 > > # or, more easily, type="Tukey" > summary(simtest(plates ~ detergent + block, data=detergent, + whichf="detergent", type="Tukey", ttype="logical")) Simultaneous tests: Tukey contrasts Call: simtest.formula(formula = plates ~ detergent + block, data = detergent, whichf = "detergent", type = "Tukey", ttype = "logical") Tukey contrasts for factor detergent, covariable: block Contrast matrix: detergent1 detergent2 detergent3 detergent4 detergent5 detergent2-detergent1 0 -1 1 0 0 0 detergent3-detergent1 0 -1 0 1 0 0 detergent4-detergent1 0 -1 0 0 1 0 detergent5-detergent1 0 -1 0 0 0 1 detergent3-detergent2 0 0 -1 1 0 0 detergent4-detergent2 0 0 -1 0 1 0 detergent5-detergent2 0 0 -1 0 0 1 detergent4-detergent3 0 0 0 -1 1 0 detergent5-detergent3 0 0 0 -1 0 1 detergent5-detergent4 0 0 0 0 -1 1 detergent2-detergent1 0 0 0 0 0 0 0 0 0 detergent3-detergent1 0 0 0 0 0 0 0 0 0 detergent4-detergent1 0 0 0 0 0 0 0 0 0 detergent5-detergent1 0 0 0 0 0 0 0 0 0 detergent3-detergent2 0 0 0 0 0 0 0 0 0 detergent4-detergent2 0 0 0 0 0 0 0 0 0 detergent5-detergent2 0 0 0 0 0 0 0 0 0 detergent4-detergent3 0 0 0 0 0 0 0 0 0 detergent5-detergent3 0 0 0 0 0 0 0 0 0 detergent5-detergent4 0 0 0 0 0 0 0 0 0 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj detergent5-detergent3 -7.933 -9.140 0.868 0.000 0.000 0.000 detergent5-detergent4 -6.533 -7.527 0.868 0.000 0.000 0.000 detergent3-detergent2 5.733 -6.606 0.868 0.000 0.000 0.000 detergent5-detergent1 -4.333 -4.993 0.868 0.000 0.001 0.001 detergent4-detergent2 4.333 -4.993 0.868 0.000 0.001 0.001 detergent3-detergent1 3.600 -4.148 0.868 0.001 0.003 0.003 detergent5-detergent2 -2.200 -2.535 0.868 0.022 0.044 0.043 detergent4-detergent1 2.200 -2.535 0.868 0.022 0.044 0.043 detergent2-detergent1 -2.133 -2.458 0.868 0.026 0.052 0.050 detergent4-detergent3 -1.400 -1.613 0.868 0.126 0.126 0.126 > > > > > cleanEx(); ..nameEx <- "litter" > > ### * litter > > flush(stderr()); flush(stdout()) > > ### Name: litter > ### Title: Litter Weights Data Set > ### Aliases: litter > ### Keywords: datasets > > ### ** Examples > > data(litter) > > # define the contrast matrix, either completly > C <- matrix(c(0,0,0,0,0,0,1,1,1,.75,.384,.887,-1,0,0,.25,.37,.113, + 0,-1,0,-.25,.246,-.339,0,0,-1,-.75,-1,-.661,0,0,0, + 0,0,0,0,0,0,0,0,0), ncol=7) > # numerate the contrasts > rownames(C) <- paste("C", 1:nrow(C), sep="") > > # simultaneous confidence intervals for above comparisons > # in an ANCOVA > summary(simint(weight ~ dose + gesttime + number, data=litter, + alternative="greater", cmatrix=C)) Simultaneous 95% confidence intervals: user-defined contrasts Call: simint.formula(formula = weight ~ dose + gesttime + number, data = litter, alternative = "greater", cmatrix = C) user-defined contrasts for factor dose, covariables: gesttime + number Contrast matrix: (Intercept) dose0 dose5 dose50 dose500 gesttime number C1 0 1.000 -1.000 0.000 0.000 0 0 C2 0 1.000 0.000 -1.000 0.000 0 0 C3 0 1.000 0.000 0.000 -1.000 0 0 C4 0 0.750 0.250 -0.250 -0.750 0 0 C5 0 0.384 0.370 0.246 -1.000 0 0 C6 0 0.887 0.113 -0.339 -0.661 0 0 Absolute Error Tolerance: 0.001 95 % quantile: 2.221 Coefficients: Estimate 5 % -- t value Std.Err. p raw p Bonf p adj C1 3.352 0.486 Inf 2.597 1.291 0.006 0.035 0.021 C2 2.291 -0.682 Inf 1.712 1.338 0.046 0.275 0.138 C3 2.675 -0.288 Inf 2.005 1.334 0.024 0.147 0.079 C4 1.741 -0.576 Inf 1.669 1.043 0.050 0.299 0.148 C5 0.871 -1.643 Inf 0.770 1.132 0.222 1.000 0.501 C6 2.166 -0.218 Inf 2.018 1.074 0.024 0.143 0.077 > > # or for dose only > D <- matrix(c(1,1,1,.75,.384,.887,-1,0,0,.25,.37,.113,0,-1,0,-.25, + .246,-.339,0,0,-1,-.75,-1,-.661),ncol = 4) > rownames(D) <- paste("D", 1:nrow(D), sep="") > > summary(simint(weight ~ dose + gesttime + number, data=litter, + whichf="dose", alternative="greater", cmatrix=D)) Simultaneous 95% confidence intervals: user-defined contrasts Call: simint.formula(formula = weight ~ dose + gesttime + number, data = litter, whichf = "dose", alternative = "greater", cmatrix = D) user-defined contrasts for factor dose, covariables: gesttime + number Contrast matrix: [,1] [,2] [,3] [,4] [,5] [,6] [,7] D1 0 1.000 -1.000 0.000 0.000 0 0 D2 0 1.000 0.000 -1.000 0.000 0 0 D3 0 1.000 0.000 0.000 -1.000 0 0 D4 0 0.750 0.250 -0.250 -0.750 0 0 D5 0 0.384 0.370 0.246 -1.000 0 0 D6 0 0.887 0.113 -0.339 -0.661 0 0 Absolute Error Tolerance: 0.001 95 % quantile: 2.221 Coefficients: Estimate 5 % -- t value Std.Err. p raw p Bonf p adj D1 3.352 0.485 Inf 2.597 1.291 0.006 0.035 0.021 D2 2.291 -0.682 Inf 1.712 1.338 0.046 0.275 0.138 D3 2.675 -0.288 Inf 2.005 1.334 0.024 0.147 0.079 D4 1.741 -0.576 Inf 1.669 1.043 0.050 0.299 0.149 D5 0.871 -1.643 Inf 0.770 1.132 0.222 1.000 0.501 D6 2.166 -0.218 Inf 2.018 1.074 0.024 0.143 0.077 > > > > > cleanEx(); ..nameEx <- "mginv" > > ### * mginv > > flush(stderr()); flush(stdout()) > > ### Name: mginv > ### Title: Generalized Inverse of a Matrix > ### Aliases: mginv > ### Keywords: algebra > > ### ** Examples > > ## Not run: > ##D # The function is currently defined as > ##D function(X, tol = sqrt(.Machine$double.eps)) > ##D { > ##D ## Generalized Inverse of a Matrix > ##D dnx <- dimnames(X) > ##D if(is.null(dnx)) dnx <- vector("list", 2) > ##D s <- svd(X) > ##D nz <- s$d > tol * s$d[1] > ##D structure( > ##D if(any(nz)) s$v[, nz] %*% (t(s$u[, nz])/s$d[nz]) else X, > ##D dimnames = dnx[2:1]) > ##D } > ## End(Not run) > > > cleanEx(); ..nameEx <- "recovery" > > ### * recovery > > flush(stderr()); flush(stdout()) > > ### Name: recovery > ### Title: Recovery Time Data Set > ### Aliases: recovery > ### Keywords: datasets > > ### ** Examples > > data(recovery) > > # one-sided simultaneous confidence intervals for Dunnett > # in the one-way layout > ci <- simint(minutes ~ blanket, data=recovery, conf.level=0.9, + alternative="less",eps=0.0001) > summary(ci) Simultaneous 90% confidence intervals: Dunnett contrasts Call: simint.formula(formula = minutes ~ blanket, data = recovery, conf.level = 0.9, alternative = "less", eps = 1e-04) Dunnett contrasts for factor blanket Contrast matrix: blanketb0 blanketb1 blanketb2 blanketb3 blanketb1-blanketb0 0 -1 1 0 0 blanketb2-blanketb0 0 -1 0 1 0 blanketb3-blanketb0 0 -1 0 0 1 Absolute Error Tolerance: 1e-04 90 % quantile: 1.8431 Coefficients: Estimate -- 90 % t value Std.Err. p raw p Bonf p adj blanketb1-blanketb0 -2.1333 -Inf 0.8226 -1.3302 1.6038 0.0958 0.2874 0.2412 blanketb2-blanketb0 -7.4667 -Inf -4.5107 -4.6556 1.6038 0.0000 0.0001 0.0001 blanketb3-blanketb0 -1.6667 -Inf -0.0359 -1.8837 0.8848 0.0337 0.1012 0.0924 > plot(ci,cex.axis=1.5, lwd=1.5) > > # 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) > print(test) Simultaneous confidence intervals: user-defined contrasts Call: simint.formula(formula = minutes ~ blanket, data = recovery, conf.level = 0.9, alternative = "less", eps = 1e-04, cmatrix = C) 90 % confidence intervals Estimate -- 90 % C1 -2.1333 -Inf 0.8225 C2 -7.4667 -Inf -4.5108 C3 -1.6667 -Inf -0.0360 > > # same results, but more detailed information using the summary method > summary(test) Simultaneous 90% confidence intervals: user-defined contrasts Call: simint.formula(formula = minutes ~ blanket, data = recovery, conf.level = 0.9, alternative = "less", eps = 1e-04, cmatrix = C) user-defined contrasts for factor blanket Contrast matrix: (Intercept) blanketb0 blanketb1 blanketb2 blanketb3 C1 0 -1 1 0 0 C2 0 -1 0 1 0 C3 0 -1 0 0 1 Absolute Error Tolerance: 1e-04 90 % quantile: 1.843 Coefficients: Estimate -- 90 % t value Std.Err. p raw p Bonf p adj C1 -2.1333 -Inf 0.8225 -1.3302 1.6038 0.0958 0.2874 0.2412 C2 -7.4667 -Inf -4.5108 -4.6556 1.6038 0.0000 0.0001 0.0001 C3 -1.6667 -Inf -0.0360 -1.8837 0.8848 0.0337 0.1012 0.0924 > > > > cleanEx(); ..nameEx <- "respiratory" > > ### * respiratory > > flush(stderr()); flush(stdout()) > > ### Name: respiratory > ### Title: Respiratory Health Data Set > ### Aliases: respiratory > ### Keywords: datasets > > ### ** Examples > > data(respiratory) > > # compute the contrast matrix in several steps > # overall active vs. placebo > CA <- c(0, 13, 0, 11, 0, 13, 0, 17, 0) > CP <- c(0, 0, 14, 0, 12, 0, 19, 0, 12) > CA <- CA/sum(CA) > CP <- CP/sum(CP) > C1 <- CP-CA > > # for older subgroup only > CAO <- c( 0, 13, 0, 0, 0, 13, 0, 0, 0 ) > CPO <- c( 0, 0, 14, 0, 0, 0, 19, 0, 0 ) > CAO <- CAO/sum(CAO) > CPO <- CPO/sum(CPO) > C2 <- CPO - CAO > > # for younger subgroup only > CAY <- c( 0, 0, 0, 11, 0, 0, 0, 17, 0 ) > CPY <- c( 0, 0, 0, 0, 12, 0, 0, 0, 12 ) > CAY <- CAY/sum(CAY) > CPY <- CPY/sum(CPY) > C3 <- CPY - CAY > > # subgroup with inital good health > CAG <- c( 0, 13, 0, 11, 0, 0, 0, 0, 0 ) > CPG <- c( 0, 0, 14, 0, 12, 0, 0, 0, 0 ) > CAG <- CAG/sum(CAG) > CPG <- CPG/sum(CPG) > C4 <- CPG - CAG > > # subgroup with inital poor health > CAP <- c( 0, 0, 0, 0, 0, 13, 0, 17, 0 ) > CPP <- c( 0, 0, 0, 0, 0, 0, 19, 0, 12 ) > CAP <- CAP/sum(CAP) > CPP <- CPP/sum(CPP) > C5 <- CPP - CAP > > # all 4 subgroup combinations of age and initial health condition > C6 <- c( 0, -1, 1, 0, 0, 0, 0, 0, 0 ) > C7 <- c( 0, 0, 0, 0, 0, -1, 1, 0, 0 ) > C8 <- c( 0, 0, 0, -1, 1, 0, 0, 0, 0 ) > C9 <- c( 0, 0, 0, 0, 0, 0, 0, -1, 1 ) > > # contrast matrix, note: first column is zero and corresponds to the > # intercept implicitly given in the formula > C <- rbind(C1, C2, C3, C4, C5, C6, C7, C8, C9) > > # numerate the contrasts > colnames(C) <- NULL > rownames(C) <- c("Overall", "Older", "Younger", "Good Init", "Poor Init", + "Old x Good", "Old x Poor", "Young x Good", "Young x Poor") > > # remove the intercept (not needed, simtest can deal with contrast matrices > # with and without a column of zeros for the intercept > C <- C[,-1] > > summary(simtest(Score ~ Treatment:AgeGroup:InitHealth, + data=respiratory, ttype="logical", + alternative="greater", cmatrix=C)) Simultaneous tests: user-defined contrasts Call: simtest.formula(formula = Score ~ Treatment:AgeGroup:InitHealth, data = respiratory, ttype = "logical", alternative = "greater", cmatrix = C) user-defined contrasts for factor Treatment:AgeGroup:InitHealth Contrast matrix: [,1] [,2] [,3] [,4] [,5] [,6] Overall 0 0.2407407 -0.2456140 0.2037037 -0.2105263 0.2407407 Older 0 0.5000000 -0.4242424 0.0000000 0.0000000 0.5000000 Younger 0 0.0000000 0.0000000 0.3928571 -0.5000000 0.0000000 Good Init 0 0.5416667 -0.5384615 0.4583333 -0.4615385 0.0000000 Poor Init 0 0.0000000 0.0000000 0.0000000 0.0000000 0.4333333 Old x Good 0 1.0000000 -1.0000000 0.0000000 0.0000000 0.0000000 Old x Poor 0 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 Young x Good 0 0.0000000 0.0000000 1.0000000 -1.0000000 0.0000000 Young x Poor 0 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 [,7] [,8] [,9] Overall -0.3333333 0.3148148 -0.2105263 Older -0.5757576 0.0000000 0.0000000 Younger 0.0000000 0.6071429 -0.5000000 Good Init 0.0000000 0.0000000 0.0000000 Poor Init -0.6129032 0.5666667 -0.3870968 Old x Good 0.0000000 0.0000000 0.0000000 Old x Poor -1.0000000 0.0000000 0.0000000 Young x Good 0.0000000 0.0000000 0.0000000 Young x Poor 0.0000000 1.0000000 -1.0000000 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj Older 1.073 4.074 0.191 0.000 0.000 0.000 Overall 0.736 3.857 0.263 0.000 0.001 0.001 Poor Init 0.861 3.346 0.279 0.001 0.003 0.003 Old x Good 1.219 3.149 0.284 0.001 0.006 0.006 Young x Poor 0.870 2.296 0.257 0.012 0.059 0.049 Old x Poor 0.815 2.255 0.387 0.013 0.059 0.049 Good Init 0.609 2.142 0.362 0.017 0.059 0.049 Younger 0.319 1.142 0.419 0.128 0.256 0.195 Young x Good -0.105 -0.249 0.379 0.598 0.598 0.598 > > > > cleanEx(); ..nameEx <- "simint" > > ### * simint > > flush(stderr()); flush(stdout()) > > ### Name: simint > ### Title: Simultaneous Confidence Intervals > ### Aliases: simint simint.default simint.formula simint.lm simint.glm > ### Keywords: htest > > ### ** 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 90% confidence intervals: Dunnett contrasts Call: simint.formula(formula = minutes ~ blanket, data = recovery, type = "Dunnett", conf.level = 0.9, alternative = "less", eps = 1e-04) Dunnett contrasts for factor blanket Contrast matrix: blanketb0 blanketb1 blanketb2 blanketb3 blanketb1-blanketb0 0 -1 1 0 0 blanketb2-blanketb0 0 -1 0 1 0 blanketb3-blanketb0 0 -1 0 0 1 Absolute Error Tolerance: 1e-04 90 % quantile: 1.8431 Coefficients: Estimate -- 90 % t value Std.Err. p raw p Bonf p adj blanketb1-blanketb0 -2.1333 -Inf 0.8226 -1.3302 1.6038 0.0958 0.2874 0.2412 blanketb2-blanketb0 -7.4667 -Inf -4.5107 -4.6556 1.6038 0.0000 0.0001 0.0001 blanketb3-blanketb0 -1.6667 -Inf -0.0359 -1.8837 0.8848 0.0337 0.1012 0.0924 > > # alternatively via a prespecified linear model > lmmod <- lm(minutes ~ blanket, data=recovery, contrasts=list(blanket = + "contr.Dunnett")) > summary(simint(lmmod, psubset=2:4, conf.level=0.9, + alternative="less",eps=0.0001)) Simultaneous 90% confidence intervals: model contrasts model contrasts for factor Contrast matrix: [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] 0 0 1 Absolute Error Tolerance: 1e-04 90 % quantile: 1.843 Coefficients: Estimate -- 90 % t value Std.Err. p raw p Bonf p adj blanketb1-b0 -2.1333 -Inf 0.8225 -1.3302 1.6038 0.0958 0.2874 0.2412 blanketb2-b0 -7.4667 -Inf -4.5108 -4.6556 1.6038 0.0000 0.0001 0.0001 blanketb3-b0 -1.6667 -Inf -0.0360 -1.8837 0.8848 0.0337 0.1012 0.0924 > > # Tukey confidence intervals, compare with TukeyHSD > > data(warpbreaks) > fm1 <- aov(breaks ~ wool + tension, data = warpbreaks) > tHSD <- TukeyHSD(fm1, "tension", ordered = FALSE) > print(tHSD) Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = breaks ~ wool + tension, data = warpbreaks) $tension diff lwr upr M-L -10.000000 -19.35342 -0.6465793 H-L -14.722222 -24.07564 -5.3688015 H-M -4.722222 -14.07564 4.6311985 > > mcHSD <- simint(breaks ~ wool + tension, data = warpbreaks, + whichf="tension", type="Tukey") > print(mcHSD) Simultaneous confidence intervals: Tukey contrasts Call: simint.formula(formula = breaks ~ wool + tension, data = warpbreaks, whichf = "tension", type = "Tukey") 95 % confidence intervals Estimate 2.5 % 97.5 % tensionM-tensionL -10.000 -19.352 -0.648 tensionH-tensionL -14.722 -24.074 -5.370 tensionH-tensionM -4.722 -14.074 4.630 > plot(mcHSD) > > ## Don't show: > tlow <- tHSD$tension[,2] > tupp <- tHSD$tension[,3] > mclow <- mcHSD$conf.int[,1] > mcupp <- mcHSD$conf.int[,2] > > a <- round(tlow, 1) > b <- round(mclow, 1) > names(a) <- NULL > names(b) <- NULL > stopifnot(all.equal(a, b)) > > a <- round(tupp, 1) > b <- round(mcupp, 1) > names(a) <- NULL > names(b) <- NULL > stopifnot(all.equal(a, b)) > ## End Don't show > > > > > cleanEx(); ..nameEx <- "simtest" > > ### * simtest > > flush(stderr()); flush(stdout()) > > ### Name: simtest > ### Title: Simultaneous Comparisons > ### Aliases: simtest.default simtest.formula simtest.lm simtest.glm simtest > ### Keywords: htest > > ### ** 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: Tukey contrasts Call: simtest.formula(formula = response ~ trt, data = cholesterol, type = "Tukey", ttype = "logical") Contrast matrix: trt1time trt2times trt4times trtdrugD trtdrugE trt2times-trt1time 0 -1 1 0 0 0 trt4times-trt1time 0 -1 0 1 0 0 trtdrugD-trt1time 0 -1 0 0 1 0 trtdrugE-trt1time 0 -1 0 0 0 1 trt4times-trt2times 0 0 -1 1 0 0 trtdrugD-trt2times 0 0 -1 0 1 0 trtdrugE-trt2times 0 0 -1 0 0 1 trtdrugD-trt4times 0 0 0 -1 1 0 trtdrugE-trt4times 0 0 0 -1 0 1 trtdrugE-trtdrugD 0 0 0 0 -1 1 Adjusted P-Values p adj trtdrugE-trt1time 0.000 trtdrugE-trt2times 0.000 trtdrugD-trt1time 0.000 trtdrugE-trt4times 0.000 trt4times-trt1time 0.000 trtdrugD-trt2times 0.000 trtdrugE-trtdrugD 0.001 trt2times-trt1time 0.042 trt4times-trt2times 0.042 trtdrugD-trt4times 0.044 > > # some examples for the formula interface, statistically non-sense! > > # response > y <- rnorm(21) > > # three factors > f1 <- factor(c(rep(c("A", "B", "C"), 7))) > f2 <- factor(c(rep("D", 10), rep("E", 11))) > > # and one continuous covariable > x <- rnorm(21) > testdata <- cbind(as.data.frame(y), f1, f2, x) > > # one factor only > summary(simtest(y ~ f1)) Simultaneous tests: Dunnett contrasts Call: simtest.formula(formula = y ~ f1) Dunnett contrasts for factor f1 Contrast matrix: f1A f1B f1C f1B-f1A 0 -1 1 0 f1C-f1A 0 -1 0 1 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj f1B-f1A 0.434 -0.769 0.565 0.452 0.904 0.666 f1C-f1A 0.036 -0.064 0.565 0.950 0.950 0.950 > > # one factor only, the same > summary(simtest(y ~ f1, data=testdata)) Simultaneous tests: Dunnett contrasts Call: simtest.formula(formula = y ~ f1, data = testdata) Dunnett contrasts for factor f1 Contrast matrix: f1A f1B f1C f1B-f1A 0 -1 1 0 f1C-f1A 0 -1 0 1 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj f1B-f1A 0.434 -0.769 0.565 0.452 0.904 0.666 f1C-f1A 0.036 -0.064 0.565 0.950 0.950 0.950 > > # and a continuous covariable > summary(simtest(y ~ f1 + x, data=testdata)) Simultaneous tests: Dunnett contrasts Call: simtest.formula(formula = y ~ f1 + x, data = testdata) Dunnett contrasts for factor f1, covariable: x Contrast matrix: f1A f1B f1C f1B-f1A 0 -1 1 0 0 f1C-f1A 0 -1 0 1 0 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj f1C-f1A -0.438 -0.817 0.521 0.425 0.851 0.627 f1B-f1A 0.078 -0.150 0.536 0.883 0.883 0.883 > > # without intercept > summary(simtest(y ~ f1 + x - 1, data=testdata)) Simultaneous tests: Dunnett contrasts Call: simtest.formula(formula = y ~ f1 + x - 1, data = testdata) Dunnett contrasts for factor f1, covariable: x Contrast matrix: f1A f1B f1C f1B-f1A -1 1 0 0 f1C-f1A -1 0 1 0 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj f1C-f1A -0.438 -0.817 0.521 0.425 0.851 0.627 f1B-f1A 0.078 -0.150 0.536 0.883 0.883 0.883 > > # with an additional factor as covariable > # use `whichf' to specify the term in the model to > # calculate p-values or confidence intervals for > summary(simtest(y ~ f1 + f2 + x - 1, data=testdata, whichf="f1")) Simultaneous tests: Dunnett contrasts Call: simtest.formula(formula = y ~ f1 + f2 + x - 1, data = testdata, whichf = "f1") Dunnett contrasts for factor f1, covariables: f2 + x Contrast matrix: f1A f1B f1C f1B-f1A -1 1 0 0 0 f1C-f1A -1 0 1 0 0 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj f1C-f1A -0.435 -0.791 0.535 0.440 0.881 0.645 f1B-f1A 0.086 -0.161 0.550 0.874 0.881 0.874 > > # with interaction terms > summary(simtest(y ~ f1*f2 + x - 1, data=testdata, whichf="f1")) Simultaneous tests: Dunnett contrasts Call: simtest.formula(formula = y ~ f1 * f2 + x - 1, data = testdata, whichf = "f1") Dunnett contrasts for factor f1, covariables: f2 +x + f1:f2 Contrast matrix: f1A f1B f1C f1B-f1A -1 1 0 0 0 0 0 0 f1C-f1A -1 0 1 0 0 0 0 0 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj f1C-f1A -0.422 -0.539 0.781 0.598 1 0.82 f1B-f1A 0.280 -0.358 0.783 0.725 1 0.82 > > # inference about the interactions term > # either Tetrade contrasts > summary(simtest(y ~ f1:f2, data=testdata, type="Tetrade")) Simultaneous tests: Tetrade contrasts Call: simtest.formula(formula = y ~ f1:f2, data = testdata, type = "Tetrade") Tetrade contrasts for factor f1:f2 Contrast matrix: [,1] [,2] [,3] [,4] [,5] [,6] [,7] (11-12)-(21-22) 0 1 -1 -1 1 0 0 (11-12)-(31-32) 0 1 -1 0 0 -1 1 (21-22)-(31-32) 0 0 0 1 -1 -1 1 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj (21-22)-(31-32) -1.658 -1.644 1.009 0.121 0.363 0.259 (11-12)-(21-22) 0.967 -0.959 1.009 0.353 0.706 0.543 (11-12)-(31-32) -0.691 -0.685 1.009 0.504 0.706 0.543 > > # or a user-defined contrast matrix > # note: this is a contrast matrix for the interactions only, > # the column for the intercept is added automatically > simtest(y ~ f1:f2, data=testdata, cmatrix=diag(6)) Warning in simtest.default(y = c(1.5778917949044, 0.596234109318454, -1.17357694087136, : at least one contrast is not estimable Simultaneous tests: user-defined contrasts Call: simtest.formula(formula = y ~ f1:f2, data = testdata, cmatrix = diag(6)) Contrast matrix: [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0 1 0 0 0 0 0 [2,] 0 0 1 0 0 0 0 [3,] 0 0 0 1 0 0 0 [4,] 0 0 0 0 1 0 0 [5,] 0 0 0 0 0 1 0 [6,] 0 0 0 0 0 0 1 Adjusted P-Values p adj [1,] 0.968 [2,] 0.991 [3,] 0.991 [4,] 0.991 [5,] 0.991 [6,] 0.991 > > # works too, if the column for the intercept is included > summary(simtest(y ~ f1:f2, data=testdata, cmatrix=cbind(0, diag(6)))) Warning in simtest.default(y = c(1.5778917949044, 0.596234109318454, -1.17357694087136, : at least one contrast is not estimable Simultaneous tests: user-defined contrasts Call: simtest.formula(formula = y ~ f1:f2, data = testdata, cmatrix = cbind(0, diag(6))) user-defined contrasts for factor f1:f2 Contrast matrix: (Intercept) f1A:f2D f1B:f2D f1C:f2D f1A:f2E f1B:f2E f1C:f2E [1,] 0 1 0 0 0 0 0 [2,] 0 0 1 0 0 0 0 [3,] 0 0 0 1 0 0 0 [4,] 0 0 0 0 1 0 0 [5,] 0 0 0 0 0 1 0 [6,] 0 0 0 0 0 0 1 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj [1,] -0.411 -0.685 0.531 0.504 1 0.968 [2,] 0.243 -0.457 0.600 0.654 1 0.992 [3,] 0.254 -0.424 0.600 0.678 1 0.992 [4,] -0.245 -0.409 0.600 0.688 1 0.992 [5,] -0.143 -0.269 0.531 0.792 1 0.992 [6,] 0.045 -0.084 0.531 0.934 1 0.992 > > # additional covariable > summary(simtest(y ~ f1:f2 + x, data=testdata, cmatrix=diag(6))) Warning in simtest.default(y = c(1.5778917949044, 0.596234109318454, -1.17357694087136, : at least one contrast is not estimable Simultaneous tests: user-defined contrasts Call: simtest.formula(formula = y ~ f1:f2 + x, data = testdata, cmatrix = diag(6)) user-defined contrasts for factor f1:f2, covariable: x Contrast matrix: [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [1,] 0 0 1 0 0 0 0 0 [2,] 0 0 0 1 0 0 0 0 [3,] 0 0 0 0 1 0 0 0 [4,] 0 0 0 0 0 1 0 0 [5,] 0 0 0 0 0 0 1 0 [6,] 0 0 0 0 0 0 0 1 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj [1,] -0.420 -0.812 0.485 0.430 1 0.935 [2,] 0.384 -0.715 0.537 0.486 1 0.945 [3,] -0.318 -0.593 0.536 0.562 1 0.951 [4,] 0.104 -0.215 0.554 0.833 1 0.994 [5,] 0.080 -0.145 0.488 0.887 1 0.994 [6,] -0.029 -0.060 0.517 0.953 1 0.994 > > # again with intercept and covariables in included in cmatrix > summary(simtest(y ~ f1:f2 + x, data=testdata, + cmatrix=cbind(0, diag(6), 0))) Warning in simtest.default(y = c(1.5778917949044, 0.596234109318454, -1.17357694087136, : at least one contrast is not estimable Simultaneous tests: user-defined contrasts Call: simtest.formula(formula = y ~ f1:f2 + x, data = testdata, cmatrix = cbind(0, diag(6), 0)) user-defined contrasts for factor f1:f2, covariable: x Contrast matrix: (Intercept) x f1A:f2D f1B:f2D f1C:f2D f1A:f2E f1B:f2E f1C:f2E [1,] 0 1 0 0 0 0 0 0 [2,] 0 0 1 0 0 0 0 0 [3,] 0 0 0 1 0 0 0 0 [4,] 0 0 0 0 1 0 0 0 [5,] 0 0 0 0 0 1 0 0 [6,] 0 0 0 0 0 0 1 0 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj [1,] -0.585 -2.223 0.263 0.043 0.259 0.203 [2,] 0.384 -0.715 0.485 0.486 1.000 0.945 [3,] -0.318 -0.593 0.537 0.562 1.000 0.952 [4,] 0.104 -0.215 0.536 0.833 1.000 0.994 [5,] 0.080 -0.145 0.554 0.887 1.000 0.994 [6,] -0.029 -0.060 0.488 0.953 1.000 0.994 > > > > > cleanEx(); ..nameEx <- "tire" > > ### * tire > > flush(stderr()); flush(stdout()) > > ### Name: tire > ### Title: Tire Wear Data Set > ### Aliases: tire > ### Keywords: datasets > > ### ** Examples > > data(tire) > 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 for difference of regression functions > summary(simint(cost ~ make*mph, data=tire, + cmatrix=C, eps=0.001)) Simultaneous 95% confidence intervals: user-defined contrasts Call: simint.formula(formula = cost ~ make * mph, data = tire, cmatrix = C, eps = 0.001) user-defined contrasts for factor make, covariables: mph + make:mph Contrast matrix: (Intercept) makeA makeB mph makeA:mph makeB:mph C1 0 1 -1 0 10 -10 C2 0 1 -1 0 15 -15 C3 0 1 -1 0 20 -20 C4 0 1 -1 0 25 -25 C5 0 1 -1 0 30 -30 C6 0 1 -1 0 35 -35 C7 0 1 -1 0 40 -40 C8 0 1 -1 0 45 -45 C9 0 1 -1 0 50 -50 C10 0 1 -1 0 55 -55 C11 0 1 -1 0 60 -60 C12 0 1 -1 0 65 -65 C13 0 1 -1 0 70 -70 Absolute Error Tolerance: 0.001 95 % quantile: 2.647 Coefficients: Estimate 2.5 % 97.5 % t value Std.Err. p raw p Bonf p adj C1 -4.107 -6.527 -1.686 -4.492 0.914 0.000 0.005 0.001 C2 -3.454 -5.594 -1.314 -4.272 0.808 0.001 0.008 0.002 C3 -2.801 -4.681 -0.921 -3.945 0.710 0.001 0.015 0.004 C4 -2.148 -3.798 -0.499 -3.448 0.623 0.003 0.043 0.011 C5 -1.496 -2.958 -0.033 -2.707 0.552 0.016 0.202 0.045 C6 -0.843 -2.181 0.495 -1.668 0.505 0.115 1.000 0.264 C7 -0.190 -1.484 1.104 -0.389 0.489 0.703 1.000 0.924 C8 0.463 -0.875 1.801 0.916 0.505 0.373 1.000 0.653 C9 1.116 -0.347 2.578 2.019 0.552 0.061 0.787 0.152 C10 1.768 0.119 3.418 2.838 0.623 0.012 0.154 0.035 C11 2.421 0.541 4.301 3.410 0.710 0.004 0.047 0.011 C12 3.074 0.934 5.214 3.802 0.808 0.002 0.020 0.005 C13 3.727 1.306 6.147 4.076 0.914 0.001 0.011 0.003 > > > > cleanEx(); ..nameEx <- "waste" > > ### * waste > > flush(stderr()); flush(stdout()) > > ### Name: waste > ### Title: Industrial Waste Data Set > ### Aliases: waste > ### Keywords: datasets > > ### ** 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.1)) Simultaneous 95% confidence intervals: Tetrade contrasts Call: simint.formula(formula = waste ~ envir:temp, data = waste, type = "Tetrade", eps = 0.1) Tetrade contrasts for factor envir:temp Contrast matrix: [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] (11-12)-(21-22) 0 1 -1 0 -1 1 0 0 0 0 0 0 (11-13)-(21-23) 0 1 0 -1 -1 0 1 0 0 0 0 0 (12-13)-(22-23) 0 0 1 -1 0 -1 1 0 0 0 0 0 (11-12)-(31-32) 0 1 -1 0 0 0 0 -1 1 0 0 0 (11-13)-(31-33) 0 1 0 -1 0 0 0 -1 0 1 0 0 (12-13)-(32-33) 0 0 1 -1 0 0 0 0 -1 1 0 0 (11-12)-(41-42) 0 1 -1 0 0 0 0 0 0 0 -1 1 (11-13)-(41-43) 0 1 0 -1 0 0 0 0 0 0 -1 0 (12-13)-(42-43) 0 0 1 -1 0 0 0 0 0 0 0 -1 (11-12)-(51-52) 0 1 -1 0 0 0 0 0 0 0 0 0 (11-13)-(51-53) 0 1 0 -1 0 0 0 0 0 0 0 0 (12-13)-(52-53) 0 0 1 -1 0 0 0 0 0 0 0 0 (21-22)-(31-32) 0 0 0 0 1 -1 0 -1 1 0 0 0 (21-23)-(31-33) 0 0 0 0 1 0 -1 -1 0 1 0 0 (22-23)-(32-33) 0 0 0 0 0 1 -1 0 -1 1 0 0 (21-22)-(41-42) 0 0 0 0 1 -1 0 0 0 0 -1 1 (21-23)-(41-43) 0 0 0 0 1 0 -1 0 0 0 -1 0 (22-23)-(42-43) 0 0 0 0 0 1 -1 0 0 0 0 -1 (21-22)-(51-52) 0 0 0 0 1 -1 0 0 0 0 0 0 (21-23)-(51-53) 0 0 0 0 1 0 -1 0 0 0 0 0 (22-23)-(52-53) 0 0 0 0 0 1 -1 0 0 0 0 0 (31-32)-(41-42) 0 0 0 0 0 0 0 1 -1 0 -1 1 (31-33)-(41-43) 0 0 0 0 0 0 0 1 0 -1 -1 0 (32-33)-(42-43) 0 0 0 0 0 0 0 0 1 -1 0 -1 (31-32)-(51-52) 0 0 0 0 0 0 0 1 -1 0 0 0 (31-33)-(51-53) 0 0 0 0 0 0 0 1 0 -1 0 0 (32-33)-(52-53) 0 0 0 0 0 0 0 0 1 -1 0 0 (41-42)-(51-52) 0 0 0 0 0 0 0 0 0 0 1 -1 (41-43)-(51-53) 0 0 0 0 0 0 0 0 0 0 1 0 (42-43)-(52-53) 0 0 0 0 0 0 0 0 0 0 0 1 [,13] [,14] [,15] [,16] (11-12)-(21-22) 0 0 0 0 (11-13)-(21-23) 0 0 0 0 (12-13)-(22-23) 0 0 0 0 (11-12)-(31-32) 0 0 0 0 (11-13)-(31-33) 0 0 0 0 (12-13)-(32-33) 0 0 0 0 (11-12)-(41-42) 0 0 0 0 (11-13)-(41-43) 1 0 0 0 (12-13)-(42-43) 1 0 0 0 (11-12)-(51-52) 0 -1 1 0 (11-13)-(51-53) 0 -1 0 1 (12-13)-(52-53) 0 0 -1 1 (21-22)-(31-32) 0 0 0 0 (21-23)-(31-33) 0 0 0 0 (22-23)-(32-33) 0 0 0 0 (21-22)-(41-42) 0 0 0 0 (21-23)-(41-43) 1 0 0 0 (22-23)-(42-43) 1 0 0 0 (21-22)-(51-52) 0 -1 1 0 (21-23)-(51-53) 0 -1 0 1 (22-23)-(52-53) 0 0 -1 1 (31-32)-(41-42) 0 0 0 0 (31-33)-(41-43) 1 0 0 0 (32-33)-(42-43) 1 0 0 0 (31-32)-(51-52) 0 -1 1 0 (31-33)-(51-53) 0 -1 0 1 (32-33)-(52-53) 0 0 -1 1 (41-42)-(51-52) 0 -1 1 0 (41-43)-(51-53) -1 -1 0 1 (42-43)-(52-53) -1 0 -1 1 Absolute Error Tolerance: 0.1 95 % quantile: 3.6 Coefficients: Estimate 2.5 % 97.5 % t value Std.Err. p raw p Bonf p adj (11-12)-(21-22) -1.8 -7.3 3.7 -1.2 1.5 0.3 1.0 1.0 (11-13)-(21-23) -0.2 -5.7 5.3 -0.1 1.5 0.9 1.0 1.0 (12-13)-(22-23) 1.6 -3.9 7.0 1.0 1.5 0.3 1.0 1.0 (11-12)-(31-32) -2.4 -7.8 3.1 -1.5 1.5 0.1 1.0 0.8 (11-13)-(31-33) -1.7 -7.2 3.8 -1.1 1.5 0.3 1.0 0.9 (12-13)-(32-33) 0.7 -4.8 6.1 0.4 1.5 0.7 1.0 1.0 (11-12)-(41-42) 2.1 -3.4 7.6 1.4 1.5 0.2 1.0 0.9 (11-13)-(41-43) 3.7 -1.8 9.2 2.4 1.5 0.0 0.8 0.3 (12-13)-(42-43) 1.6 -3.8 7.1 1.1 1.5 0.3 1.0 1.0 (11-12)-(51-52) 1.3 -4.2 6.7 0.8 1.5 0.4 1.0 1.0 (11-13)-(51-53) 2.0 -3.5 7.5 1.3 1.5 0.2 1.0 0.9 (12-13)-(52-53) 0.7 -4.7 6.2 0.5 1.5 0.6 1.0 1.0 (21-22)-(31-32) -0.6 -6.1 4.9 -0.4 1.5 0.7 1.0 1.0 (21-23)-(31-33) -1.5 -7.0 4.0 -1.0 1.5 0.3 1.0 1.0 (22-23)-(32-33) -0.9 -6.4 4.6 -0.6 1.5 0.6 1.0 1.0 (21-22)-(41-42) 3.9 -1.6 9.3 2.5 1.5 0.0 0.7 0.3 (21-23)-(41-43) 3.9 -1.5 9.4 2.6 1.5 0.0 0.6 0.3 (22-23)-(42-43) 0.1 -5.4 5.6 0.1 1.5 1.0 1.0 1.0 (21-22)-(51-52) 3.0 -2.4 8.5 2.0 1.5 0.1 1.0 0.6 (21-23)-(51-53) 2.2 -3.3 7.7 1.4 1.5 0.2 1.0 0.9 (22-23)-(52-53) -0.8 -6.3 4.6 -0.5 1.5 0.6 1.0 1.0 (31-32)-(41-42) 4.4 -1.0 9.9 2.9 1.5 0.0 0.3 0.1 (31-33)-(41-43) 5.4 0.0 10.9 3.5 1.5 0.0 0.1 0.0 (32-33)-(42-43) 1.0 -4.5 6.5 0.6 1.5 0.5 1.0 1.0 (31-32)-(51-52) 3.6 -1.8 9.1 2.4 1.5 0.0 1.0 0.3 (31-33)-(51-53) 3.7 -1.8 9.2 2.4 1.5 0.0 0.9 0.3 (32-33)-(52-53) 0.1 -5.4 5.6 0.1 1.5 1.0 1.0 1.0 (41-42)-(51-52) -0.8 -6.3 4.7 -0.5 1.5 0.6 1.0 1.0 (41-43)-(51-53) -1.7 -7.2 3.7 -1.1 1.5 0.3 1.0 1.0 (42-43)-(52-53) -0.9 -6.4 4.6 -0.6 1.5 0.6 1.0 1.0 > summary(simtest(waste ~ envir:temp, data=waste, + type="Tetrade", eps=0.1)) Simultaneous tests: Tetrade contrasts Call: simtest.formula(formula = waste ~ envir:temp, data = waste, type = "Tetrade", eps = 0.1) Tetrade contrasts for factor envir:temp Contrast matrix: [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] (11-12)-(21-22) 0 1 -1 0 -1 1 0 0 0 0 0 0 (11-13)-(21-23) 0 1 0 -1 -1 0 1 0 0 0 0 0 (12-13)-(22-23) 0 0 1 -1 0 -1 1 0 0 0 0 0 (11-12)-(31-32) 0 1 -1 0 0 0 0 -1 1 0 0 0 (11-13)-(31-33) 0 1 0 -1 0 0 0 -1 0 1 0 0 (12-13)-(32-33) 0 0 1 -1 0 0 0 0 -1 1 0 0 (11-12)-(41-42) 0 1 -1 0 0 0 0 0 0 0 -1 1 (11-13)-(41-43) 0 1 0 -1 0 0 0 0 0 0 -1 0 (12-13)-(42-43) 0 0 1 -1 0 0 0 0 0 0 0 -1 (11-12)-(51-52) 0 1 -1 0 0 0 0 0 0 0 0 0 (11-13)-(51-53) 0 1 0 -1 0 0 0 0 0 0 0 0 (12-13)-(52-53) 0 0 1 -1 0 0 0 0 0 0 0 0 (21-22)-(31-32) 0 0 0 0 1 -1 0 -1 1 0 0 0 (21-23)-(31-33) 0 0 0 0 1 0 -1 -1 0 1 0 0 (22-23)-(32-33) 0 0 0 0 0 1 -1 0 -1 1 0 0 (21-22)-(41-42) 0 0 0 0 1 -1 0 0 0 0 -1 1 (21-23)-(41-43) 0 0 0 0 1 0 -1 0 0 0 -1 0 (22-23)-(42-43) 0 0 0 0 0 1 -1 0 0 0 0 -1 (21-22)-(51-52) 0 0 0 0 1 -1 0 0 0 0 0 0 (21-23)-(51-53) 0 0 0 0 1 0 -1 0 0 0 0 0 (22-23)-(52-53) 0 0 0 0 0 1 -1 0 0 0 0 0 (31-32)-(41-42) 0 0 0 0 0 0 0 1 -1 0 -1 1 (31-33)-(41-43) 0 0 0 0 0 0 0 1 0 -1 -1 0 (32-33)-(42-43) 0 0 0 0 0 0 0 0 1 -1 0 -1 (31-32)-(51-52) 0 0 0 0 0 0 0 1 -1 0 0 0 (31-33)-(51-53) 0 0 0 0 0 0 0 1 0 -1 0 0 (32-33)-(52-53) 0 0 0 0 0 0 0 0 1 -1 0 0 (41-42)-(51-52) 0 0 0 0 0 0 0 0 0 0 1 -1 (41-43)-(51-53) 0 0 0 0 0 0 0 0 0 0 1 0 (42-43)-(52-53) 0 0 0 0 0 0 0 0 0 0 0 1 [,13] [,14] [,15] [,16] (11-12)-(21-22) 0 0 0 0 (11-13)-(21-23) 0 0 0 0 (12-13)-(22-23) 0 0 0 0 (11-12)-(31-32) 0 0 0 0 (11-13)-(31-33) 0 0 0 0 (12-13)-(32-33) 0 0 0 0 (11-12)-(41-42) 0 0 0 0 (11-13)-(41-43) 1 0 0 0 (12-13)-(42-43) 1 0 0 0 (11-12)-(51-52) 0 -1 1 0 (11-13)-(51-53) 0 -1 0 1 (12-13)-(52-53) 0 0 -1 1 (21-22)-(31-32) 0 0 0 0 (21-23)-(31-33) 0 0 0 0 (22-23)-(32-33) 0 0 0 0 (21-22)-(41-42) 0 0 0 0 (21-23)-(41-43) 1 0 0 0 (22-23)-(42-43) 1 0 0 0 (21-22)-(51-52) 0 -1 1 0 (21-23)-(51-53) 0 -1 0 1 (22-23)-(52-53) 0 0 -1 1 (31-32)-(41-42) 0 0 0 0 (31-33)-(41-43) 1 0 0 0 (32-33)-(42-43) 1 0 0 0 (31-32)-(51-52) 0 -1 1 0 (31-33)-(51-53) 0 -1 0 1 (32-33)-(52-53) 0 0 -1 1 (41-42)-(51-52) 0 -1 1 0 (41-43)-(51-53) -1 -1 0 1 (42-43)-(52-53) -1 0 -1 1 Absolute Error Tolerance: 0.1 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj (31-33)-(41-43) 5.4 -3.5 1.5 0.0 0.1 0.0 (31-32)-(41-42) 4.4 -2.9 1.5 0.0 0.3 0.2 (21-23)-(41-43) 3.9 -2.6 1.5 0.0 0.6 0.3 (21-22)-(41-42) 3.9 -2.5 1.5 0.0 0.6 0.3 (11-13)-(41-43) 3.7 -2.4 1.5 0.0 0.7 0.3 (31-33)-(51-53) 3.7 -2.4 1.5 0.0 0.7 0.3 (31-32)-(51-52) 3.6 -2.4 1.5 0.0 0.8 0.3 (21-22)-(51-52) 3.0 -2.0 1.5 0.1 1.0 0.5 (11-12)-(31-32) -2.4 -1.5 1.5 0.1 1.0 0.8 (21-23)-(51-53) 2.2 -1.4 1.5 0.2 1.0 0.8 (11-12)-(41-42) 2.1 -1.4 1.5 0.2 1.0 0.9 (11-13)-(51-53) 2.0 -1.3 1.5 0.2 1.0 0.9 (11-12)-(21-22) -1.8 -1.2 1.5 0.3 1.0 0.9 (41-43)-(51-53) -1.7 -1.1 1.5 0.3 1.0 0.9 (11-13)-(31-33) -1.7 -1.1 1.5 0.3 1.0 0.9 (12-13)-(42-43) 1.6 -1.1 1.5 0.3 1.0 0.9 (12-13)-(22-23) 1.6 -1.0 1.5 0.3 1.0 0.9 (21-23)-(31-33) -1.5 -1.0 1.5 0.3 1.0 0.9 (11-12)-(51-52) 1.3 -0.8 1.5 0.4 1.0 1.0 (32-33)-(42-43) 1.0 -0.6 1.5 0.5 1.0 1.0 (22-23)-(32-33) -0.9 -0.6 1.5 0.6 1.0 1.0 (42-43)-(52-53) -0.9 -0.6 1.5 0.6 1.0 1.0 (22-23)-(52-53) -0.8 -0.5 1.5 0.6 1.0 1.0 (41-42)-(51-52) -0.8 -0.5 1.5 0.6 1.0 1.0 (12-13)-(52-53) 0.7 -0.5 1.5 0.6 1.0 1.0 (12-13)-(32-33) 0.7 -0.4 1.5 0.7 1.0 1.0 (21-22)-(31-32) -0.6 -0.4 1.5 0.7 1.0 1.0 (11-13)-(21-23) -0.2 -0.1 1.5 0.9 1.0 1.0 (22-23)-(42-43) 0.1 -0.1 1.5 1.0 1.0 1.0 (32-33)-(52-53) 0.1 -0.1 1.5 1.0 1.0 1.0 > > > > > ### *