## Code preceded by R> in the paper "The cg Package for Comparison of Groups" ## NOTE that the some of the code in Section 6.2 is actually run inside ## S-PLUS. See commented block. ## print() wrappers are added in case this is run as a single script. ## cat() calls are also added to identify sections of code execution. ## Section 2.1 cat("\nSection 2.1\n") library("cg") ## successful load required for any of the following code to run print(canine) ## Section 2.2 cat("\nSection 2.2\n") print(gmcsfcens, row.names=FALSE) ## Section 3.2 cat("\nSection 3.2\n") print(canine.listfmt) ## Section 4.1 cat("\nSection 4.1\n") canine.data <- prepareCGOneFactorData(canine, format = "groupcolumns", analysisname = "Canine", endptname = "Prostate Volume", endptunits = expression(plain(cm)^3), digits = 1, logscale = TRUE) descriptiveTable(canine.data) canine.fit <- fit(canine.data) errorBarGraph(canine.fit) comparisonsTable(canine.fit, model = "olsonly") varianceGraph(canine.fit, model = "olsonly") qqGraph(canine.fit, model = "extended") canine.samplesize <- samplesizeTable(canine.fit, direction = "increasing", mmdvec = c(10, 25, 50, 75, 100)) samplesizeGraph(canine.samplesize) ## Section 4.2 cat("\nSection 4.2\n") gmcsfcens.data <- prepareCGOneFactorData(gmcsfcens, format = "groupcolumns", analysisname = "cytokine", endptname = "GM-CSF (pg/ml)", logscale = TRUE, digits = 1) boxplot(gmcsfcens.data) descriptiveTable(gmcsfcens.data) kmGraph(gmcsfcens.data, distfcn = "cumulative") gmcsfcens.fit <- fit(gmcsfcens.data) gmcsfcens.comps <- comparisonsTable(gmcsfcens.fit, type = "allgroupstocontrol", refgrp = "PBS/Tg 197") comparisonsGraph(gmcsfcens.comps) comparisonsTable(gmcsfcens.fit, type = "allgroupstocontrol", refgrp = "PBS/Tg 197", mcadjust = TRUE) ## Section 6.1 cat("\nSection 6.1\n") ## Abbreviated paper listing head(gmcsfcens.listfmt1); cat("\n...\n"); tail(gmcsfcens.listfmt1) ## Full listing print(gmcsfcens.listfmt1) gmcsfcens.listfmt1.data <- prepareCGOneFactorData(gmcsfcens.listfmt1, format = "listed", analysisname = "cytokine", endptname = "GM-CSF (pg/ml)", logscale = TRUE, digits = 1) ## Abbreviated paper listing head(gmcsfcens.listfmt2); cat("\n...\n"); tail(gmcsfcens.listfmt2) ## Full listing print(gmcsfcens.listfmt2) gmcsfcens.listfmt2.data <- prepareCGOneFactorData(gmcsfcens.listfmt2, format = "listed", analysisname = "cytokine", endptname = "GM-CSF (pg/ml)", logscale = TRUE, digits = 1, leftcensor = TRUE) ## Abbreviated paper listing head(gmcsfcens.listfmt3); cat("\n...\n"); tail(gmcsfcens.listfmt3) ## Full listing print(gmcsfcens.listfmt3) gmcsfcens.listfmt3.data <- prepareCGOneFactorData(gmcsfcens.listfmt3, format = "listed", analysisname = "cytokine", endptname = "GM-CSF (pg/ml)", logscale = TRUE, digits = 1) ## Section 6.2 cat("\nSection 6.2\n") group.s <- c(3, 5, 10) n.s <- c(3, 5, 10) balanced.configs <- vector("list", length = length(group.s) * length(n.s)) i <- 0 for(g in group.s) { for(n in n.s) { i <- i + 1 balanced.configs[[i]] <- rep(n, g) } } data.frame(configs = sapply(balanced.configs, function(x) paste(x, collapse = ","))) unbalance.theconfigs <- function(x, delta) { lapply(x[unlist(lapply(balanced.configs, function(x) { all(x > 3) } ))], function(x) { add <- sample((-delta):delta, size = length(x), replace = TRUE) while(length(unique(add)) == 1) { add <- sample((-delta):delta, size = length(x), replace = TRUE) } return(x + add) }) } set.seed(47) mild.unbalanced.configs <- unbalance.theconfigs(balanced.configs, 1) unbalanced.configs <- unbalance.theconfigs(balanced.configs, 2) ## Output objects for use in S-PLUS dump(c("balanced.configs","mild.unbalanced.configs", "unbalanced.configs"), file = "p1Rdump.ssc") ########################################################################## ## Code for use in S-PLUS ########################################################################## if(Hmisc:::.SV4.) { source("p1Rdump.ssc") crit.dunnett <- function(x) { lapply(x, function(x) { qdunnett(p = .95, k = length(x), df = sum(x) - length(x), nvec = x, control = 1) }) } balanced.dunnett <- crit.dunnett(balanced.configs) mild.unbalanced.dunnett <- crit.dunnett(mild.unbalanced.configs) unbalanced.dunnett <- crit.dunnett(unbalanced.configs) ## Output objects for use back in R dump(c("balanced.dunnett","mild.unbalanced.dunnett", "unbalanced.dunnett"), file = "p1sscdump.R", oldStyle = TRUE) } ########################################################################## ## End of code for use in S-PLUS ########################################################################## source("p1sscdump.R") crit.tukey <- function(x) { lapply(x, function(x) { qtukey(.95, length(x), sum(x) - length(x)) / sqrt(2) }) } balanced.tukey <- crit.tukey(balanced.configs) mild.unbalanced.tukey <- crit.tukey(mild.unbalanced.configs) unbalanced.tukey <- crit.tukey(unbalanced.configs) crit.multcomp <- function(x, type) { lapply(x, function(x) { d <- data.frame(grp = factor(rep(1:length(x), x), 1:length(x)), endpt = rnorm(sum(x))) f <- aov(endpt ~ grp, data = d) m <- glht(f, linfct = mcp(grp = type)) ca <- adjusted_calpha() return(as.vector(ca(m, 0.95))) }) } cat("\nPlease wait while critical points are calculated from multcomp::glht...\n") flush.console() balanced.tukey.multcomp <- crit.multcomp(balanced.configs, type = "Tukey") mild.unbalanced.tukey.multcomp <- crit.multcomp(mild.unbalanced.configs, type = "Tukey") unbalanced.tukey.multcomp <- crit.multcomp(unbalanced.configs, type = "Tukey") balanced.dunnett.multcomp <- crit.multcomp(balanced.configs, type = "Dunnett") mild.unbalanced.dunnett.multcomp <- crit.multcomp(mild.unbalanced.configs, type = "Dunnett") unbalanced.dunnett.multcomp <- crit.multcomp(unbalanced.configs, type = "Dunnett") cat("\nBalanced Tukey\n") print( data.frame(numgrps = rep(c(3, 5, 10), each = 3), samplesizes = sapply(balanced.configs, function(x) { paste(x, collapse = ",") } ), multcomp = unlist(balanced.tukey.multcomp), crit.tukey = unlist(balanced.tukey)) ) cat("\nBalanced Dunnett\n") print( data.frame(numgrps = rep(c(3, 5, 10), each = 3), samplesizes = sapply(balanced.configs, function(x) { paste(x, collapse = ",") } ), multcomp = unlist(balanced.dunnett.multcomp), crit.dunnett = unlist(balanced.dunnett)) ) cat("\nMildly Unbalanced Tukey\n") print( data.frame(numgrps = rep(c(3, 5, 10), each = 2), samplesizes = sapply(mild.unbalanced.configs, function(x) { paste(x, collapse = ",") } ), multcomp = unlist(mild.unbalanced.tukey.multcomp), crit.tukey = unlist(mild.unbalanced.tukey)) ) cat("\nMildly Unbalanced Dunnett\n") print( data.frame(numgrps = rep(c(3, 5, 10), each = 2), samplesizes = sapply(mild.unbalanced.configs, function(x) { paste(x, collapse = ",") } ), multcomp = unlist(mild.unbalanced.dunnett.multcomp), crit.dunnett = unlist(mild.unbalanced.dunnett)) ) cat("\nUnbalanced Tukey\n") print( data.frame(numgrps = rep(c(3, 5, 10), each = 2), samplesizes = sapply(unbalanced.configs, function(x) { paste(x, collapse = ",") } ), multcomp = unlist(unbalanced.tukey.multcomp), crit.tukey = unlist(unbalanced.tukey)) ) cat("\nUnbalanced Dunnett\n") print( data.frame(numgrps = rep(c(3, 5, 10), each = 2), samplesizes = sapply(unbalanced.configs, function(x) { paste(x, collapse = ",") } ), multcomp = unlist(unbalanced.dunnett.multcomp), crit.dunnett = unlist(unbalanced.dunnett)) )