Nothing
R/p03FormalOneFactorData.R
In cg: Compare Groups, Analytically and Graphically
Defines functions
comparisons
pairwisecompsmatrix
errorbargraph
grpsummary
comparisonsgraph
validErrorDf
validComparisonType
validEstimates
validGrpNames
validN
validCnames
validAddPct
Documented in
comparisons
comparisonsgraph
errorbargraph
grpsummary
pairwisecompsmatrix
validComparisonType
validErrorDf
validEstimates
validGrpNames
validN
## $Id: p03FormalOneFactorData.R 6053 2015-02-22 20:23:45Z bpikouni $
## One-Factor Unpaired Groups Case
## Formal Analysis methods for One-Factor Unpaired Groups Data
setClass("cgOneFactorGlobalTest",
representation(ols.gpval="numericOrNULL", rr.gpval="numericOrNULL",
aft.gpval="numericOrNULL", uv.gpval="numericOrNULL",
settings="list"),
prototype(ols.gpval=NULL, rr.gpval=NULL, aft.gpval=NULL,
uv.gpval=NULL,
settings=list()))
setMethod("globalTest", "cgOneFactorFit",
globalTest.cgOneFactorFit <-
function(fit, display="print", ...) {
##
## PURPOSE: Derive global F-test p-values
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
## initializations
aft <- ols <- rr <- uv <- FALSE
rr.gpval <- ols.gpval <- aft.gpval <- uv.gpval <- NULL
##
rrfit <- fit@rrfit
olsfit <- fit@olsfit
aftfit <- fit@aftfit
uvfit <- fit@uvfit
settings <- fit@settings
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
if(class(aftfit)[1]=="survreg") {
aft <- TRUE
validArgModel(...)
}
else if(class(uvfit)[1]=="gls") {
uv <- TRUE
validArgModel(...)
}
if(class(rrfit)[1]=="rlm" && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(class(olsfit)[1]=="lm" && model!="rronly" && !aft && !uv) {
ols <- TRUE
if(!rr) model <- "olsonly"
}
if(rr || ols) {
df.residual <- olsfit$df.residual
df.grpf <- length(settings$grpnames) - 1
upform <- update(formula(olsfit), . ~ grpf)
}
dfru <- olsfit$dfru
if(ols) {
## calculations mimic print.summary.lm code snippet
olsfit <- lm(upform, data=dfru)
ols.gpval <- with(summary(olsfit), pf(fstatistic[1],
fstatistic[2],
fstatistic[3], lower.tail
= FALSE))
}
if(rr) {
## if only 2 groups, we rely on the comparisonstable method for computing
## the pvalue
if(df.grpf==1) {
rr.gpval <- pairwisecompsmatrix(comparisonsTable(fit,
type="pairwisereflect")@rr.comprs,
grpnames=settings$grpnames)[3, 3]
}
else {
## As of current writing,
## The MASS summary.rlm does not compute or facilitate construction of a
## global F, or equivalently, an R^2. For now we use a placeholder method
## based on ad-hoc (witchcraft) development, which is essentially to
## re-fit a linear model with lm() and weights from the RR fit.
olsw.rrfit <- lm(upform, data=dfru, weights=rrfit$w)
rr.Fobs <- { (1 / (1 / summary.lm(olsw.rrfit)$r.squared - 1))*
(df.residual/df.grpf) }
rr.gpval <- 1 - pf(rr.Fobs, df.grpf, df.residual)
}
}
else if(aft) {
upform <- update(formula(aftfit), . ~ grpf)
aftfit <- survreg(upform, data = dfru,
dist = "gaussian", maxiter = aftfit$maxIter)
## calculations mimic print.survreg code snippet
chisq <- 2 * diff(aftfit$loglik)
df <- sum(aftfit$df) - aftfit$idf
aft.gpval <- 1 - pchisq(chisq, df)
}
else if(uv) {
upform <- update(formula(uvfit), . ~ grpf)
uvfit <- gls(model = upform, data = dfru,
weights = varIdent(form = ~1 | grpf))
uv.gpval <- anova(uvfit)$"p-value"[2]
}
x <- new("cgOneFactorGlobalTest",
ols.gpval=ols.gpval, rr.gpval=rr.gpval,
aft.gpval=aft.gpval, uv.gpval=uv.gpval,
settings=settings)
if(display=="print") {
print(x)
}
else if (display=="show"){
showDefault(x)
}
## else show nothing
invisible(x)
})
setMethod("print", "cgOneFactorGlobalTest",
print.cgOneFactorGlobalTest <-
function(x, title=NULL, endptname=NULL, ...) {
##
## PURPOSE: Semi-formatted print version of Global test p-value
## Three digits are used, and anything smaller than or equal to
## 0.0005 gets "< 0.001"
##
## NOTE: Had to use x as an argument because of the system defined
## generic. I would have preferred to use object; hence the first
## statement below.
##
object <- x
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
olsgrps <- object@ols.gpval
rrgrps <- object@rr.gpval
aftgrps <- object@aft.gpval
uvgrps <- object@uv.gpval
ols <- rr <- uv <- aft <- FALSE ## Initializations
if(!is.null(aftgrps)) {
aft <- TRUE
validArgModel(...)
}
else if(!is.null(uvgrps)) {
uv <- TRUE
validArgModel(...)
}
if(!is.null(rrgrps) && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(!is.null(olsgrps) && (model!="rronly") && !aft && !uv) {
ols <- TRUE
}
settings <- object@settings
alpha <- settings$alpha
if(is.null(title)) {
title <- paste("Global Test P-value of", settings$analysisname)
}
else {
validCharacter(title)
}
if(is.null(endptname)) {
endptname <- settings$endptname
if(!is.character(endptname)) {
endptname <- ""
}
}
else {
validCharacter(endptname)
}
cat(title,"\n")
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
if(ols) {
cat("\nClassical Least Squares Model Fit:", fmtPvalue(olsgrps), "\n")
}
if(rr) {
cat("\nResistant & Robust Model Fit:", fmtPvalue(rrgrps), "\n")
}
if(aft) {
cat("\nAccelerated Failure Time Model Fit:", fmtPvalue(aftgrps), "\n")
}
if(uv) {
cat("\nUnequal Variances Model Fit:", fmtPvalue(uvgrps), "\n")
}
invisible()
})
setMethod("show", "cgOneFactorGlobalTest",
show.cgOneFactorGlobalTest <- function(object) showDefault(object))
## Comparisons Tables
comparisons <- function(estimates,
varcovmatrix, errordf=Inf,
endptscale, mcadjust=FALSE,
alpha=0.05,
type="pairwisereflect",
contrastmatrix=NULL, n,
offset=NULL, cnames="derive",
analysisname="", endptname="",
digits=NULL,
addpct=FALSE,
display="print") {
##
## PURPOSE: Function for computing comparisons results
## and placing in a listed form.
##
## Input argument checking
if(!hasArg(estimates)) reportInvalidArg("estimates")
if(!hasArg(varcovmatrix)) reportInvalidArg("varcovmatrix")
validAlpha(alpha)
validBoolean(mcadjust)
validBoolean(addpct)
validErrorDf(errordf, varcovmatrix, n, mcadjust)
validComparisonType(type, errordf)
validCnames(type, cnames, nrow(contrastmatrix))
estimates <- validEstimates(estimates)
endptscale <- validArgMatch(endptscale, c("log","original"))
validAddPct(addpct, endptscale)
display <- validArgMatch(display, c("print","none","show"))
ngrps <- length(estimates)
grpnames <- names(estimates)
if(is.null(contrastmatrix)) {
L <- contrastMatrix(list(grpnames), type=type)
## if type="allgroupstocontrol", then the first grpname is set to be
## the reference group level.
}
else if(type=="custom") {
L <- contrastmatrix
}
## If only two groups were analyzed, then multiplicity adjustment has
## no effect. So we ensure mcadjust is set to FALSE
if(ngrps==2) { mcadjust <- FALSE }
if(mcadjust) {
set.seed(17)
multcompInform()
glhtobj <- list(object=NULL, linfct=L, coef=estimates,
vcov=varcovmatrix,
type="user-defined", alternative="two.sided",
df=errordf,
rhs=rep(0, nrow(L)))
class(glhtobj) <- "glht"
mcpci.fit <- confint(glhtobj, level=1-alpha)
mcptest.fit <- summary(glhtobj)
estdiff <- mcpci.fit$confint[,"Estimate"]
sediff <- mcptest.fit$test$sigma
lowerci <- mcpci.fit$confint[,"lwr"]
upperci <- mcpci.fit$confint[,"upr"]
pval <- mcptest.fit$test$pvalues
}
else if(!mcadjust) {
estdiff <- L %*% estimates
sediff <- sqrt(diag(L %*% varcovmatrix %*% t(L)))
if(errordf=="approx") {
## Apply Satterthwaite approximation for special case of pairwise
## or all groups to control comparisons (-1 vs 1 contrasts)
incidence <- t(apply(L, 1, function(x) { x!=0 }))
thestnderrs <- sqrt(diag(varcovmatrix))
errordf <- apply(incidence, 1,
function(x, thestnderrs, n) {
sepair <- thestnderrs[x]
npair <- n[x]
sese <- sqrt(sepair[1]^2 + sepair[2]^2)
return(sese^4/(sepair[1]^4/(npair[1] - 1) +
sepair[2]^4/(npair[2] - 1)))
}, thestnderrs=thestnderrs, n=n)
}
tcrit <- qt(1 - alpha/2, errordf)
lowerci <- estdiff - tcrit * sediff
upperci <- estdiff + tcrit * sediff
pval <- 2 * (1 - pt(abs(estdiff/sediff), errordf))
}
mcp <- data.frame(estimate=estdiff, se=sediff,
lowerci=lowerci, upperci=upperci,
pval=pval)
## Log Case:
## Transformations to Percent Difference Expression Scale
## And we need to specially handle the asymmetry of percent change,
## since decreases cannot be greater than 100%
## also note that if offset is a valid number,
## some calculations will be replaced.
if(endptscale=="log") {
## Focus only on the percent differences for now
logscalest <- mcp[,1]
if(is.null(offset)) {
mcp[, 1] <- 100 * ( exp(logscalest) - 1 ) # Point Estimate
mcp[, 2] <- 100 * ( exp(logscalest) * mcp[, 2] ) # Std Err
mcp[, 3:4] <- 100 * ( exp(mcp[, 3:4]) - 1) # Confidence Limits
}
else { ## offset
validNumeric(offset, positive=TRUE)
mA <- exp(mcp$meanA) - offset
mB <- exp(mcp$meanB) - offset
mAmBratio <- mA/mB
correctionFactor <- (mA * (mB + offset))/(mB * (mA + offset))
mcp[, 1] <- 100 * ( mAmBratio - 1 ) # Point Estimate
mcp[, 2] <- 100 * ( mAmBratio * mcp[, 2] ) # Std Err
mcp[, 3:4] <- 100 * ( correctionFactor*exp(mcp[, 3:4]) - 1) # Confidence Limits
}
}
## Fetch, arrange, and add individual component estimates
if(type!="custom") {
complist <- strsplit(row.names(L), split=" vs. ", fixed=TRUE)
compA <- sapply(complist, function(x) x[1])
compB <- sapply(complist, function(x) x[2])
compAlist <- as.list(compA)
compBlist <- as.list(compB)
compAindx <- unlist(sapply(compAlist,
function(x) which(grpnames==x)))
compBindx <- unlist(sapply(compBlist,
function(x) which(grpnames==x)))
sevec <- sqrt(diag(varcovmatrix))
mcp$meanA <- estimates[compAindx]
mcp$seA <- sevec[compAindx]
mcp$meanB <- estimates[compBindx]
mcp$seB <- sevec[compBindx]
}
else if(type=="custom") {
L.A <- t(apply(L, 1, function(x) ifelse(x > 0, x, 0)))
## For the B side of the contrast, we need to remove the minus sign
L.B <- t(apply(L, 1, function(x) ifelse(x < 0, -x, 0)))
mcp$meanA <- L.A %*% estimates
mcp$seA <- sqrt(diag(L.A %*% varcovmatrix %*% t(L.A)))
mcp$meanB <- L.B %*% estimates
mcp$seB <- sqrt(diag(L.B %*% varcovmatrix %*% t(L.B)))
if(length(cnames)==1 && cnames=="derive") {
names.ests <- names(estimates)
names.A <- apply(L.A, 1,
function(x) {
indx <- (x > 0)
coef.ests <- x[indx]
coef.ests <- ifelse(coef.ests%%1!=0,
signif(coef.ests, 3), coef.ests)
label <- paste(paste(coef.ests, names.ests[indx], sep="*"),
collapse="+")
if(length(coef.ests) > 1) {
label <- paste("(", label, ")", sep="")
}
gsub("1\\*","", label)
})
names.B <- apply(L.B, 1,
function(x) {
indx <- (x > 0)
coef.ests <- x[indx]
coef.ests <- ifelse(coef.ests%%1!=0,
signif(coef.ests, 3), coef.ests)
label <- paste(paste(coef.ests, names.ests[indx], sep="*"),
collapse="+")
if(length(coef.ests) > 1) {
label <- paste("(", label, ")", sep="")
}
gsub("1\\*","", label)
})
rownames(mcp) <- paste(names.A, names.B, sep=" vs. ")
}
}
if(endptscale=="log") {
## Now focus on the individual estimates
if(is.null(offset)) {
mcp$meanA <- exp(mcp$meanA)
mcp$meanB <- exp(mcp$meanB)
}
else { ## offset
mcp$meanA <- exp(mcp$meanA) - offset
mcp$meanB <- exp(mcp$meanB) - offset
}
mcp$seA <- mcp$meanA * mcp$seA
mcp$seB <- mcp$meanB * mcp$seB
}
if(cnames[1]!="derive") { row.names(mcp) <- cnames }
if(endptscale=="original" && addpct==TRUE) {
mcp$pctdiff <- 100*(mcp$meanA / mcp$meanB - 1)
}
if(display=="print") {
fmt.mcp <- mcp
cat(paste("Comparisons Table",
if(analysisname!="") paste("for", analysisname), "\n"))
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
diffmetric <- "Differences (A vs. B)"
if(endptscale=="log") {
diffmetric <- paste("Percent", diffmetric)
}
cat(diffmetric, "\n")
pvalfmt <- fmtPvalue(fmt.mcp$pval)
if(regexpr("Percent", diffmetric) > 0) {
fmt.mcp$estimate <- fmtPercent(fmt.mcp$estimate)
fmt.mcp$se <- fmtPercent(fmt.mcp$se)
fmt.mcp$lowerci <- fmtPercent(fmt.mcp$lowerci)
fmt.mcp$upperci <- fmtPercent(fmt.mcp$upperci)
fmt.mcp$pval <- pvalfmt
if(!is.null(digits)) {
fmt.mcp$meanA <- fround(fmt.mcp$meanA, digits=digits)
fmt.mcp$seA <- fround(fmt.mcp$seA, digits=digits)
fmt.mcp$meanB <- fround(fmt.mcp$meanB, digits=digits)
fmt.mcp$seB <- fround(fmt.mcp$seB, digits=digits)
}
{ names(fmt.mcp)[is.element(names(fmt.mcp),
c("meanA","seA",
"meanB","seB"))] <- c("geomeanA","seA",
"geomeanB","seB") }
{ names(mcp)[is.element(names(mcp), c("meanA","seA",
"meanB","seB"))] <-
c("geomeanA","seA",
"geomeanB","seB") }
}
else { ## Simple Differences
pctdiff <- mcp$pctdiff
if(!is.null(digits)) fmt.mcp <- fround(fmt.mcp, digits)
fmt.mcp$pval <- pvalfmt
if(endptscale=="original" && addpct==TRUE) {
fmt.mcp$pctdiff <- fmtPercent(pctdiff)
}
}
cat(paste(round(100*(1-alpha), 0), "% Confidence ",
"(alpha of ", alpha, ")",
if(mcadjust) {", Multiplicity Adjusted"},
"\n",
sep=""))
curwidth <- getOption("width")
on.exit(options(width=curwidth), add=TRUE)
if(curwidth < 500) { options(width=500) }
print(fmt.mcp, quote=FALSE)
}
else if (display=="show"){
showDefault(mcp)
}
## else show nothing
invisible(mcp)
}
setClass("cgOneFactorComparisonsTable",
representation(ols.comprs="dataframeMatrixOrNULL",
rr.comprs="dataframeMatrixOrNULL",
aft.comprs="dataframeMatrixOrNULL",
uv.comprs="dataframeMatrixOrNULL",
settings="list"),
prototype(ols.comprs=NULL,
rr.comprs=NULL,
aft.comprs=NULL,
uv.comprs=NULL,
settings=list()))
setMethod("comparisonsTable", "cgOneFactorFit",
comparisonsTable.cgOneFactorFit <-
function(fit,
## mcadjust=FALSE,
type="pairwisereflect",
## contrastmatrix=NULL,
## refgrp=NULL,
alpha=0.05, addpct=FALSE,
display="print", ...) {
##
## PURPOSE: Produce multiple comparisons of groups based
## on a one-factor unpaired samples fit.
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names=c("mcadjust", "contrastmatrix", "refgrp","model"))
type <- validArgMatch(type, c("pairwisereflect","pairwise",
"allgroupstocontrol", "custom"))
validAlpha(alpha)
validBoolean(addpct)
display <- validArgMatch(display, c("print","none","show"))
mcadjustarg <- getDotsArgName(dots, "mcadjustarg")
if(!is.na(mcadjustarg)) {
mcadjust <- eval(parse(text=paste("dots$", mcadjustarg, sep="")))
validBoolean(mcadjust)
}
else {
mcadjust <- FALSE
}
contrastmatrixarg <- getDotsArgName(dots, "contrastmatrix")
if(is.na(contrastmatrixarg)) {
contrastmatrix <- 1
contrastmatrix <- as.null(contrastmatrix)
}
refgrparg <- getDotsArgName(dots, "refgrp")
settings <- fit@settings
if(!is.na(refgrparg)) {
refgrp <- eval(parse(text=paste("dots$", refgrparg, sep="")))
if(!is.null(refgrp)) {
refgrp <- validArgMatch(refgrp, settings$grpnames)
}
if(type!="allgroupstocontrol") {
stop(cgMessage("If the refgrp= argument is specified",
"then the type= argument needs to be set to",
"\"allgroupstocontrol\"."))
}
}
else refgrp <- settings$refgrp
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
## initializations
aft <- ols <- rr <- uv <- FALSE
aft.comprs <- ols.comprs <- rr.comprs <- uv.comprs <- NULL
##
## settings <- fit@settings
endptscale <- settings$endptscale
rrfit <- fit@rrfit
olsfit <- fit@olsfit
aftfit <- fit@aftfit
uvfit <- fit@uvfit
grpnames <- settings$grpnames
offset <- settings$addconstant
validAddPct(addpct, endptscale)
if(class(aftfit)[1]=="survreg") {
aft <- TRUE
validArgModel(...)
}
else if(class(uvfit)[1]=="gls") {
uv <- TRUE
validArgModel(...)
thens <- with(uvfit$dfru, sapply(split(endpt, grpf), length))
}
if(class(rrfit)[1]=="rlm" && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(class(olsfit)[1]=="lm" && model!="rronly" && !aft && !uv) {
ols <- TRUE
if(!rr) model <- "olsonly"
}
numberofgrps <- length(grpnames)
grpnamesindex <- 1:numberofgrps
if(!is.null(refgrp) && type=="allgroupstocontrol") {
grpnamesindex <- c(which(grpnames==refgrp), which(grpnames!=refgrp))
}
if(rr || ols) {
df.residual <- olsfit$df.residual
}
if(ols) {
olsestimates <- olsfit$coef
varcov.olsfit <- vcov(olsfit)
ols.comprs <- comparisons(olsestimates[grpnamesindex],
varcov.olsfit[grpnamesindex, grpnamesindex],
errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
type=type,
contrastmatrix=contrastmatrix,
offset=offset, addpct=addpct, display="none")
if(mcadjust) { multcompDone("Classical Least Squares") }
}
if(rr) {
rrestimates <- rrfit$coef
summ.rrfit <- summary(rrfit, method="XtWX", ...)
stddev <- summ.rrfit$stddev
varcov.rrfit <- (stddev^2) * summ.rrfit$cov.unscaled
rr.comprs <- comparisons(rrestimates[grpnamesindex],
varcov.rrfit[grpnamesindex, grpnamesindex],
errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
type=type,
contrastmatrix=contrastmatrix,
offset=offset,
addpct=addpct, display="none")
if(mcadjust) { multcompDone("Resistant & Robust") }
}
if(aft) {
aftestimates <- aftfit$coef
## estscaleindex <- length(aftestimates) + 1
df.residual <- aftfit$df.residual
varcov.aftfit <- vcov(aftfit)
aft.comprs <- comparisons(aftestimates[grpnamesindex],
varcov.aftfit[grpnamesindex, grpnamesindex],
errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
type=type,
contrastmatrix=contrastmatrix,
offset=offset, display="none")
if(mcadjust) { multcompDone("Accelerated Failure Time") }
}
if(uv) {
uvestimates <- uvfit$coef
varcov.uvfit <- vcov(uvfit)
uv.comprs <- comparisons(uvestimates[grpnamesindex],
varcov.uvfit[grpnamesindex, grpnamesindex],
errordf="approx",
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
type=type,
contrastmatrix=contrastmatrix,
n=with(uvfit$dfru,
sapply(split(endpt, grpf),
length)[grpnamesindex]),
offset=offset, display="none")
}
## slightly different from the original fit object settings slot
settings <- list(endptscale=settings$endptscale,
analysisname=settings$analysisname,
endptname=settings$endptname,
endptlabel=makeEndptLabel(settings$endptname,
settings$endptunits),
alpha=alpha,
mcadjust=mcadjust,
digits=settings$digits,
stamps=settings$stamps,
type=switch(type,
pairwise="Pairwise Comparisons (halfset)",
pairwisereflect="All Pairwise Comparisons",
allgroupstocontrol="All Groups versus Control",
custom=""),
sandaft=if(aft) TRUE else FALSE,
addpct=addpct)
returnObj <- new("cgOneFactorComparisonsTable",
ols.comprs=ols.comprs,
rr.comprs=rr.comprs,
aft.comprs=aft.comprs,
uv.comprs=uv.comprs,
settings=settings)
if(display=="print") {
print(returnObj)
}
else if(display=="show") {
show(returnObj)
}
## else display=="none"
invisible(returnObj)
}
)
pairwisecompsmatrix <- function(comps, grpnames) {
##
## PURPOSE: Function for rearranging pairwise comparisons results
## and placing in a matrix. The output matrix is designed to be
## used as input to formatting or further processing
##
## comps is assumed to be a data frame such as
## created by the comparisons function
grpnames <- validGrpNames(grpnames)
ngrps <- length(grpnames)
if(nrow(comps)!=ngrps*(ngrps-1)) {
stop(cgMessage("The number of rows in the comps dataframe",
"needs to have entries for all possible",
"pairwise comparisons, i.e. the number of groups",
"times the number of groups - 1"))
}
##ngrps <- round(uniroot(function(x) x^2 - x - nrow(comps),
## c(2, 1000))$root, 0)
indx <- matrix(1:(ngrps^2) - as.vector(col(diag(ngrps))), ncol=ngrps)
mcp.indx <- indx[lower.tri(indx)]
mcp.opp.indx <- setdiff(1:nrow(comps), mcp.indx)
mcp <- comps[mcp.indx,]
mcp.opp <- comps[mcp.opp.indx,]
## The combined rows of mcp and mcp.opp
## pertain to all possible unique pairwise comparisons
## amongst the groups. So we create a matrix for each statistical
## quantity that is of dimension ngrps by ngrps, putting the two
## triangular parts that come from mcp and mcp.opp together
## Construct list of each individual quantity
mcp.matrices <- lapply(seq(along=mcp),
function(i, x, y, ngrps) {
z <- matrix(NA, nrow=ngrps, ncol=ngrps)
z[lower.tri(z)] <- x[[i]]
z[upper.tri(z)] <- y[[i]]
return(z)
}, x=mcp, y=mcp.opp, ngrps=ngrps)
names(mcp.matrices) <- c("diff","se","lcl","ucl","pval")
## Now we have to combine the matrices into one big table
## which is probably intended for import into something like
## EXCEL for formatting eventually. This will involve
## interleaving rows and columns and adding some buffer sections to
## columns. A tricky part comes with the
## "lower CI and upper CI" part, because we want these to now be on the
## same row for a given comparison. The dimension of the final table
## matrix will be 4*ngrps rows by 2*ngrps columns.
mcpout.nrows <- 4*ngrps
mcpout.ncols <- 2*ngrps
mcp.out <- matrix(NA, nrow=mcpout.nrows, ncol=mcpout.ncols)
## fill-in the ucl's first
mcp.out[rep( c(rep(FALSE, mcpout.nrows),
rep(c(rep(FALSE, 3), TRUE), ngrps)), ngrps ) ] <- mcp.matrices$ucl
## second, the other quantities
mcp.out[rep( c(rep(TRUE, mcpout.nrows),
rep(FALSE, mcpout.nrows) ),
ngrps) ] <- do.call("rbind",
lapply(mcp.matrices[c("diff","se",
"pval",
"lcl")],
as.vector))
## lastly, matrix labels
therownames <- character(length(grpnames)*4)
thecolnames <- character(length(grpnames)*2)
therownames[seq(1, length(therownames), by=4)] <- grpnames
thecolnames[seq(1, length(thecolnames), by=2)] <- grpnames
grpnameheaders <- list(therownames, thecolnames)
dimnames(mcp.out) <- grpnameheaders
# return
mcp.out
}
setMethod("print", "cgOneFactorComparisonsTable",
print.cgOneFactorComparisonsTable <-
function(x, digits=NULL, title=NULL, endptname=NULL, ...) {
##
## PURPOSE: Semi-formatted print version of Comparisons Table
##
## NOTE: Had to use x as an argument because of the system defined
## generic. I would have preferred to use object; hence the first
## statement below.
##
object <- x
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
olscomprs <- object@ols.comprs
rrcomprs <- object@rr.comprs
aftcomprs <- object@aft.comprs
uvcomprs <- object@uv.comprs
ols <- rr <- uv <- aft <- FALSE ## Initializations
if(!is.null(aftcomprs)) {
aft <- TRUE
validArgModel(...)
}
else if(!is.null(uvcomprs)) {
uv <- TRUE
validArgModel(...)
}
if(!is.null(rrcomprs) && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(!is.null(olscomprs) && (model!="rronly") && !aft && !uv) {
ols <- TRUE
}
settings <- object@settings
alpha <- settings$alpha
mcadjust <- settings$mcadjust
addpct <- settings$addpct
if(is.null(digits)) {
digits <- settings$digits
}
else {
digits <- validArgDigits(digits)
}
curscipen <- getOption("scipen")
on.exit(options(scipen=curscipen), add=TRUE)
options(scipen=9)
if(is.null(title)) {
title <- paste("Comparisons Table of", settings$analysisname)
}
else {
validCharacter(title)
}
if(is.null(endptname)) {
endptname <- settings$endptname
if(!is.character(endptname)) {
endptname <- ""
}
}
else {
validCharacter(endptname)
}
cat(title,"\n")
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
diffmetric <- "Differences (A vs. B)"
if(settings$endptscale=="log") {
diffmetric <- paste("Percent", diffmetric)
}
cat(diffmetric, "\n")
fmtdig <- function(x, diffmetric, digits) {
pvalfmt <- fmtPvalue(x$pval)
fmt.x <- x
if(regexpr("Percent", diffmetric) > 0) {
fmt.x$estimate <- fmtPercent(x$estimate)
fmt.x$se <- fmtPercent(x$se)
fmt.x$lowerci <- fmtPercent(x$lowerci)
fmt.x$upperci <- fmtPercent(x$upperci)
fmt.x$pval <- pvalfmt
fmt.x$meanA <- fround(x$meanA, digits=digits)
fmt.x$seA <- fround(x$seA, digits=digits)
fmt.x$meanB <- fround(x$meanB, digits=digits)
fmt.x$seB <- fround(x$seB, digits=digits)
{ names(fmt.x)[is.element(names(fmt.x), c("meanA","seA",
"meanB","seB"))] <-
c("geomeanA","seA",
"geomeanB","seB") }
}
else { ## Simple Differences
fmt.x <- fround(x, digits)
if(addpct==TRUE) {
pctdiff.fmt <- fmtPercent(x$pctdiff)
fmt.x$pctdiff <- pctdiff.fmt
}
fmt.x$pval <- pvalfmt
}
return(fmt.x)
}
informConfidence <- function() {
cat(paste(round(100*(1-alpha), 0), "% Confidence ",
"(alpha of ", alpha, ")",
if(mcadjust) {", Multiplicity Adjusted"},
"\n",
sep=""))
}
curwidth <- getOption("width")
on.exit(options(width=curwidth), add=TRUE)
if(curwidth < 500) { options(width=500) }
if(ols) {
cat("\nClassical Least Squares Model Fit\n")
informConfidence()
print(fmtdig(olscomprs, diffmetric, digits), quote=FALSE)
}
if(rr) {
cat("\nResistant & Robust Model Fit\n")
informConfidence()
print(fmtdig(rrcomprs, diffmetric, digits), quote=FALSE)
}
if(aft) {
afttitle <- "\nAccelerated Failure Time Model Fit\n"
if(settings$sandaft) {
cat(paste(afttitle,
"with Sandwich Variance-Covariance Estimate\n",
sep=""))
}
else {
cat(afttitle)
}
informConfidence()
print(fmtdig(aftcomprs, diffmetric, digits), quote=FALSE)
}
if(uv) {
cat("\nUnequal Variances Model Fit\n")
informConfidence()
print(fmtdig(uvcomprs, diffmetric, digits), quote=FALSE)
}
invisible()
})
setMethod("show", "cgOneFactorComparisonsTable",
show.cgOneFactorComparisonsTable <- function(object) showDefault(object))
## Error Bar Graph
errorbargraph <- function(estimates, centralvar,
critpoint,
endptscale="log",
analysisname="",
endptname="",
alpha=0.05,
digits=NULL,
approxstamp=FALSE,
titlestamp=TRUE,
offset=NULL,
ticklabels=NULL, ...) {
##
## PURPOSE: Construct an error bar graph based on pairwise
## multiple comparisons (method of Andrews, Sarner, and Snee, 1980)
##
## argument handling
estimates <- validEstimates(estimates)
endptscale <- validArgMatch(endptscale, c("log","original"))
validAlpha(alpha)
validBoolean(approxstamp)
validBoolean(titlestamp)
options(warn=-1)
curpar <- par(new=FALSE, mgp=c(3,0.25,0), tck=-0.010)
options(warn=0)
on.exit(par(curpar))
numberofgrps <- length(estimates)
grpnames <- names(estimates)
errorbarlength <- critpoint*sqrt(2*centralvar)/2
bardata <- list(lower=estimates - errorbarlength,
upper=estimates + errorbarlength)
logscale <- ifelse(endptscale=="log", TRUE, FALSE)
rangedata <- if(logscale) exp(unlist(bardata)) else unlist(bardata)
if(logscale) {
parmar <- par(mar=c(5, 4, 4, 3) + 0.1)
curpar$mar <- parmar$mar
errbar(1:numberofgrps, estimates/log(10),
yplus=bardata$upper/log(10),
yminus=bardata$lower/log(10),
ylab=endptname, xlab="",
pch=16,
xlim=c(0.5, numberofgrps + 0.5), axes=FALSE)
mtext("log-spaced", side=2, line=2.25, cex=0.7)
tickmarks <- setupAxisTicks(rangedata,
logscale=logscale, digits=digits,
offset=offset)
if(!is.null(ticklabels)) {
tickmarks <- makeTickMarks(ticklabels, tickmarks,
offset=offset)
}
axis(2, at=log10(tickmarks), labels=names(tickmarks),
cex.axis=0.8, adj=1, las=1)
log10endpt <- unlist(bardata)/log(10)
axis(4, at=pretty(log10endpt), pretty(log10endpt), cex.axis=0.7,
adj=0, las=1, tck=-0.0075)
mtext(side=4, text=if(is.expression(endptname)) {
catCharExpr("Log10 scale of", endptname)
} else { paste("Log10 scale of", endptname, sep=" ") },
line=2, cex=0.7, adj=0)
}
else {
errbar(1:numberofgrps, estimates,
yplus=bardata$upper,
yminus=bardata$lower,
ylab=endptname, xlab="",
pch=16,
xlim=c(0.5, numberofgrps + 0.5), axes=FALSE)
tickmarks <- setupAxisTicks(rangedata,
logscale=logscale, digits=digits)
if(!is.null(ticklabels)) {
tickmarks <- makeTickMarks(ticklabels, tickmarks,
offset=offset)
}
axis(2, at=tickmarks, labels=names(tickmarks),
cex.axis=0.8, adj=1, las=1)
}
## Axes Customization
grpnameticksettings <- setupGrpNameTicks(grpnames, 1:numberofgrps)
plotGrpNameTicks((grpnames), settings=grpnameticksettings)
minmaxTicks(if(logscale) exp(unlist(bardata)) else
unlist(bardata),
theaxis="y", logscale=logscale, digits=digits,
offset=offset)
## Annotations
if(titlestamp) {
title(main=paste("Error Bar Graph\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(alphapercent=100*alpha)
}
if(approxstamp) errorBarGraphApproximateStamp()
box()
invisible()
}
setMethod("errorBarGraph", "cgOneFactorFit",
errorBarGraph.cgOneFactorFit <-
function(fit, mcadjust=FALSE,
alpha=0.05, cgtheme=TRUE, device="single", ...) {
##
## PURPOSE: Construct an error bar graph based on pairwise
## multiple comparisons (method of Andrews, Sarner, and Snee, 1980)
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names=c("model", "ticklabels"))
if(class(fit@uvfit)[1]=="gls" || class(fit@aftfit)[1]=="survreg") {
stop(cgMessage("There is no errorBarGraph method",
"defined for a fitted model that allowed",
"unequal variances or censored observations.",
"You may wish to look",
"at the errorbargraph function",
"to construct a graph if your standard errors",
"are not too different."))
}
validBoolean(cgtheme)
validAlpha(alpha)
validBoolean(mcadjust)
ticklabelsarg <- getDotsArgName(dots, "ticklabels")
if(!is.na(ticklabelsarg)) {
ticklabels <- eval(parse(text=paste("dots$", ticklabelsarg, sep="")))
validList(ticklabels, names=c("mod","marks"),
argname="ticklabels")
}
else {
ticklabels <- NULL
}
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly", "rronly"))
}
else {
model <- "both"
}
device <- validArgMatch(device, c("single", "multiple", "ask"))
settings <- fit@settings
offset <- settings$addconstant
rrfit <- fit@rrfit
olsfit <- fit@olsfit
digits <- settings$digits
analysisname <- settings$analysisname
endptlabel <- makeEndptLabel(settings$endptname, settings$endptunits)
endptscale <- settings$endptscale
stamps <- settings$stamps
alphapercent <- round(100*alpha, 0)
confquantile <- 1 - alpha/2
grpnames <- settings$grpnames
numberofgrps <- length(grpnames)
if(numberofgrps==2) { mcadjust <- FALSE }
L <- contrastMatrix(list(grpnames), "pairwise")
ols <- rr <- FALSE ## Initializations
if(class(rrfit)[1]=="rlm" && model!="olsonly") {
rr <- TRUE
}
if(class(olsfit)[1]=="lm" && model!="rronly") {
ols <- TRUE
}
df.residual <- olsfit$df.residual
thetitle <- "Error Bar Graph"
dfru <- olsfit$dfru
if(ols) {
n <- with(dfru, tapply(endpt, grpf, length))
## Calculate the harmonic mean to account for unequal sample sizes
nharmonic <- n ## default all n equal
if(length(unique(n)) > 1) {
nharmonic <- with(olsfit$dfru,
length(n) / sum(tapply(endpt, grpf,
function(x) {
1/length(x)} )
))
}
olsestimates <- olsfit$coef
olscentralstderr <- summary(olsfit)$sigma/sqrt(nharmonic)
olscentralvar <- olscentralstderr^2
olscritpoint <-
if(!mcadjust) { qt(confquantile, df.residual) }
else {
set.seed(17)
multcompInform()
glhtobj <- list(object=NULL, linfct=L, coef=olsestimates,
vcov=vcov(olsfit),
type="user-defined", alternative="two.sided",
df=df.residual,
rhs=rep(0, nrow(L)))
class(glhtobj) <- "glht"
## assign this value:
attr(confint(glhtobj, level=1-alpha)$confint, "calpha")
}
if(mcadjust) { multcompDone("Classical Least Squares") }
}
if(rr) {
rrestimates <- rrfit$coef
summ.rrfit <- summary(rrfit, method="XtWX", correlation=FALSE)
varcov.rrfit <- (summ.rrfit$stddev^2) * summ.rrfit$cov.unscaled
## a robust / resistant estimate of center of the std errors
rrcentralvar <- rlm(diag(varcov.rrfit) ~ 1)[[1]]
rrcentralstderr <- sqrt(rrcentralvar)
rrcritpoint <-
if(!mcadjust) { qt(confquantile, df.residual) }
else {
set.seed(17)
multcompInform()
glhtobj <- list(object=NULL, linfct=L, coef=olsestimates,
vcov=varcov.rrfit,
type="user-defined", alternative="two.sided",
df=df.residual,
rhs=rep(0, nrow(L)))
class(glhtobj) <- "glht"
## assign this value:
attr(confint(glhtobj, level=1-alpha)$confint, "calpha")
}
if(mcadjust) { multcompDone("Resistant & Robust") }
}
if(rr && ols && is.element(model, "both") && device=="single") {
ols.dfr <- data.frame(type=rep("Classical", numberofgrps),
grpf=factorInSeq(grpnames), estimate=olsestimates,
lower=olsestimates -
olscritpoint*olscentralstderr/sqrt(2),
upper=olsestimates +
olscritpoint*olscentralstderr/sqrt(2))
rr.dfr <- data.frame(type=rep("Resistant & Robust", numberofgrps),
grpf=factorInSeq(grpnames), estimate=rrestimates,
lower=rrestimates -
rrcritpoint*rrcentralstderr/sqrt(2),
upper=rrestimates +
rrcritpoint*rrcentralstderr/sqrt(2))
all.dfr <- rbind(ols.dfr, rr.dfr)
thetitle <- paste(thetitle, "s", sep="")
theapproxmsg <-
if(n[1]!=nharmonic[1]) {
"Error bars are approximate."
}
else {
"Error bars in Resistant & Robust panel are approximate."
}
cgDevice(cgtheme=cgtheme)
trellispanelstg <- trellis.par.get("clip")$panel
trellis.par.set("clip", list(panel="off"))
on.exit(trellis.par.set("clip", list(panel=trellispanelstg)),
add=TRUE)
trellisparstg2 <- trellis.par.get("layout.widths")$ylab
trellis.par.set("layout.widths", list(ylab=3))
on.exit(trellis.par.set("layout.widths", list(panel=trellisparstg2)),
add=TRUE)
trellisparstg3 <- trellis.par.get("axis.components")
trellis.par.set("axis.components",
list(left=list(pad1=0.2, pad2=2),
bottom=list(pad1=0.2, pad2=1),
top=list(pad1=1, pad2=1),
right=list(pad1=0.5, pad2=1)
))
on.exit(trellis.par.set("axis.components", trellisparstg3),
add=TRUE)
ally <- unlist(all.dfr[, c("estimate","lower","upper")])
if(endptscale=="log") {
thegraph <- xyplot(Cbind(log10(exp(estimate)),
log10(exp(lower)),
log10(exp(upper))) ~ as.numeric(grpf)
| type,
data=all.dfr,
digits=digits,
panel = function(x, y, ...) {
grpnameticksettings <- setupGrpNameTicks(grpnames,
grplocation=
1:length(grpnames),
cexinit=0.8,
cexthreshold=0.5,
grid=TRUE)
plotGrpNameTicks(grpnames, grpnameticksettings, grid=TRUE)
if(is.element(panel.number(),
c(1, nlevels(all.dfr$type)))) {
panelside <- ifelse(panel.number()==1,
"left", "right")
ycex <- 0.6
tickmarks <- setupAxisTicks(exp(ally),
logscale=TRUE,
grid=TRUE,
digits=digits,
offset=offset,
ycex=ycex)
if(!is.null(ticklabels)) {
tickmarks <- makeTickMarks(ticklabels, tickmarks,
offset=offset)
}
panel.axis(side=panelside,
at=log10(tickmarks),
labels=names(tickmarks),
tck=0.20, text.cex=ycex,
rot=0,
outside=TRUE)
}
panel.xYplot(x, y, subscripts=FALSE,
method="bars",
label.curves=FALSE, ...)
yrange <- range(unwind(attr(y, "other")))
panel.text(x=rep(0, 2), y=yrange,
labels=paste("",
c(fround(min(10^yrange), digits),
fround(max(10^yrange), digits))),
col="blue", adj=0, cex=0.6)
},
layout=c(2,1), aspect=1, pch=16,
col="black",
as.table=TRUE,
xlim=c(0, numberofgrps + 0.5),
ylim=rangeExtend(log10(exp(ally))),
ylab=list(cex=0.7,
label=endptlabel),
xlab="",
scales=list(x=list(at=1:numberofgrps,
labels=NULL, tck=c(0.15, 0), axs="i"),
y=list(labels=NULL, tck=0)),
main=list(label=paste(thetitle, "\n",
settings$analysisname, sep=""), cex=1.1),
par.strip.text=list(cex=0.7)
)
}
else {
thegraph <- xyplot(Cbind(estimate,
lower,
upper) ~ as.numeric(grpf) | type,
data=all.dfr,
digits=digits,
panel = function(x, y, ...) {
grpnameticksettings <- setupGrpNameTicks(grpnames,
grplocation=
1:length(grpnames),
cexinit=0.8,
cexthreshold=0.5,
grid=TRUE)
plotGrpNameTicks(grpnames, grpnameticksettings, grid=TRUE)
if(is.element(panel.number(),
c(1, nlevels(all.dfr$type)))) {
panelside <- ifelse(panel.number()==1,
"left", "right")
ycex <- 0.6
tickmarks <- setupAxisTicks(ally,
logscale=FALSE,
grid=TRUE,
digits=digits,
offset=offset,
ycex=ycex)
if(!is.null(ticklabels)) {
tickmarks <- makeTickMarks(ticklabels, tickmarks,
offset=offset)
}
panel.axis(side=panelside,
at=tickmarks,
labels=names(tickmarks),
tck=0.20, text.cex=ycex,
rot=0,
outside=TRUE)
}
panel.xYplot(x, y, subscripts=FALSE,
method="bars",
label.curves=FALSE, ...)
yrange <- range(unwind(attr(y, "other")))
panel.text(x=rep(0, 2), y=yrange,
labels=paste("",
c(fround(min(yrange), digits),
fround(max(yrange), digits))),
col="blue", adj=0, cex=0.6)
},
layout=c(2,1), aspect=1, pch=16,
col="black",
as.table=TRUE,
xlim=c(0, numberofgrps + 0.5),
ylim=rangeExtend(ally),
ylab=list(cex=0.7,
label=endptlabel),
xlab="",
scales=list(x=list(at=1:numberofgrps,
labels=NULL, tck=c(0.15, 0), axs="i"),
y=list(labels=NULL, tck=0)),
main=list(label=paste(thetitle, "\n",
settings$analysisname, sep=""), cex=1.1),
par.strip.text=list(cex=0.7)
)
}
print(thegraph) ## , position=c(0,0,0.95,1))
if(stamps) graphStampCG()
errorBarGraphStamp(mcadjust, alphapercent, grid=TRUE)
errorBarGraphApproximateStamp(grid=TRUE, msg=theapproxmsg)
}
else if((model=="olsonly" || (ols && !rr && model=="both")) &&
device=="single") {
errorbargraph(olsestimates, olscentralvar,
olscritpoint,
endptscale,
analysisname,
endptname=endptlabel,
alpha=alpha,
digits=digits,
approxstamp={nharmonic[1]!=n[1]},
titlestamp=FALSE, offset=offset,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Error Bar Graph, Classical analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(mcadjust, alphapercent)
}
else if(model=="rronly" && rr && device=="single") {
errorbargraph(rrestimates, rrcentralvar,
rrcritpoint,
endptscale,
analysisname,
endptname=endptlabel,
alpha=alpha,
digits=digits,
approxstamp=TRUE,
titlestamp=FALSE,
offset=offset,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Error Bar Graph, Resistant & Robust analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(mcadjust, alphapercent)
}
else if(rr && ols &&
is.element(device, c("ask","multiple")) &&
is.element(model, "both")) {
device <- validArgMatch(device, c("multiple", "ask"))
if(device=="ask") {
op <- par(ask = TRUE)
on.exit(par(op), add=TRUE)
}
errorbargraph(olsestimates, olscentralvar,
olscritpoint,
endptscale,
analysisname,
endptname=endptlabel,
alpha=alpha,
digits=digits,
approxstamp={nharmonic[1]!=n[1]},
titlestamp=FALSE, offset=offset,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Error Bar Graph, Classical analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(mcadjust, alphapercent)
if(device=="multiple") {
if(!is.null(dots$model)) dots$model <- NULL ## since we only
## want dots arguments for trellis.device in next call
do.call("cgDevice", c(list(new=TRUE), dots))
cat(cgMessage("A new graphics device has been generated",
"to hold the Resistant & Robust",
"errorBarGraph version.",
"The Classical Least Squares version is on the previous",
"device.\n",
warning=TRUE))
}
errorbargraph(rrestimates, rrcentralvar,
rrcritpoint,
endptscale,
analysisname,
endptname=endptlabel,
alpha=alpha,
digits=digits,
approxstamp=TRUE,
titlestamp=FALSE,
offset=offset,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Error Bar Graph, Resistant & Robust analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(mcadjust, alphapercent)
}
else {
stop(cgMessage("The chosen device and model arguments",
"are not compatible either with each other",
"or with the fitted model(s) in the fit object.",
seeHelpFile("errorBarGraph")))
}
invisible()
}
)
grpsummary <- function(estimates,
varcovmatrix,
n, errordf=Inf,
endptscale="log", mcadjust=FALSE,
alpha=0.05,
offset=NULL,
analysisname="", endptname="",
digits=NULL, display="print", ...) {
##
## PURPOSE: Function for computing group summary results.
## The output matrix is designed to be
## used as input to formatting or further processing
validBoolean(mcadjust)
n <- validN(n, estimates)
estimates <- validEstimates(estimates)
endptscale <- validArgMatch(endptscale, c("log","original"))
validAlpha(alpha)
display <- validArgMatch(display, c("print","none","show"))
if(!is.numeric(errordf) || any(errordf < 1)) {
stop(cgMessage("The errordf argument needs to be numeric",
"and (each value) set to be 1 or greater."))
}
ngrps <- length(estimates)
grpnames <- names(estimates)
## If only two groups were analyzed, then multiplicity adjustment has
## no effect. So we ensure mcadjust is set to FALSE
if(ngrps==2) { mcadjust <- FALSE }
if(mcadjust) {
if(length(errordf) > 1) {
stop(cgMessage("The errordf can only be a single value",
"when mcadjust=FALSE.",
seeHelpFile("grpsummary")))
}
set.seed(17)
multcompInform()
glhtobj <- list(object=NULL, linfct=diag(ngrps),
coef=estimates,
vcov=varcovmatrix,
type="user-defined", alternative="two.sided",
df=errordf,
rhs=rep(0, nrow(diag(ngrps))))
class(glhtobj) <- "glht"
mcpci.fit <- confint(glhtobj, level=1-alpha)
mcptest.fit <- summary(glhtobj)
se <- as.vector(mcptest.fit)$test$sigma
lowerci <- mcpci.fit$confint[,"lwr"]
upperci <- mcpci.fit$confint[,"upr"]
}
else if(!mcadjust) {
se <- sqrt(diag(varcovmatrix))
tcrit <- qt(1 - alpha/2, errordf)
lowerci <- estimates - tcrit * se
upperci <- estimates + tcrit * se
}
thetable <- data.frame(n=n,
estimate=estimates, se=se,
lowerci=lowerci, upperci=upperci)
row.names(thetable) <- grpnames
## Handle log case if needed
## Transformations to Percent Difference Expression Scale
if(endptscale=="log") {
logscale.est <- thetable[, 2]
logscale.sd <- thetable[, 3]
thetable[, 2] <- exp(logscale.est) # Point Estimate
thetable[, 3] <- exp(logscale.est) * logscale.sd # Std Err
if(any(logscale.sd > 0.50)) {
warning(cgMessage("There is at least one",
"group standard error in the log scale",
"that exceeds 0.50, so the",
"estimated geometric mean standard",
"errors may be nonsensical and should be",
"cautiously regarded.")
)
}
thetable[, 4:5] <- exp(thetable[, 4:5]) # Confidence Limits
}
if(!is.null(offset)) {
validNumeric(offset, positive=TRUE)
thetable[, c(2,4,5)] <- (thetable[, c(2,4,5)] - offset)
}
if(display=="print") {
fmt.table <- thetable
cat(paste("Group Summary Table",
if(analysisname!="") paste("for", analysisname), "\n"))
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
metric <- "Means"
if(endptscale=="log") {
metric <- paste("Geometric", metric)
}
else {
metric <- paste("Arithmetic", metric)
}
cat(metric, "\n")
if(!is.null(digits)) {
fmt.table$estimate <- fround(fmt.table$estimate, digits)
fmt.table$se <- fround(fmt.table$se, digits)
fmt.table$lowerci <- fround(fmt.table$lowerci, digits)
fmt.table$upperci <- fround(fmt.table$upperci, digits)
}
fmt.table$n <- fround(fmt.table$n, 0)
cat(paste(round(100*(1-alpha), 0), "% Confidence ",
"(alpha of ", alpha, ")",
if(mcadjust) {", Multiplicity Adjusted"},
"\n",
sep=""))
curwidth <- getOption("width")
on.exit(options(width=curwidth), add=TRUE)
if(curwidth < 500) { options(width=500) }
print(fmt.table, quote=FALSE)
return(invisible(thetable))
}
else if (display=="show") {
showDefault(thetable)
}
## else show nothing
invisible(thetable)
}
setClass("cgOneFactorGrpSummaryTable",
representation(ols.grps="dataframeOrNULL",
rr.grps="dataframeOrNULL",
aft.grps="dataframeOrNULL",
uv.grps="dataframeOrNULL",
settings="list"),
prototype(ols.grps=NULL, rr.grps=NULL, aft.grps=NULL, uv.grps=NULL,
settings=list()))
setMethod("grpSummaryTable", "cgOneFactorFit",
function(fit,
mcadjust=FALSE,
alpha=0.05, display="print", ...) {
##
## PURPOSE: Summarize each group distribution based
## on the model fit.
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
settings <- fit@settings
validAlpha(alpha)
validBoolean(mcadjust)
display <- validArgMatch(display, c("print","none","show"))
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly", "rronly"))
}
else {
model <- "both"
}
##
## initializations
aft <- ols <- rr <- uv <- FALSE
aft.grps <- ols.grps <- rr.grps <- uv.grps <- NULL
##
endptscale <- settings$endptscale
rrfit <- fit@rrfit
olsfit <- fit@olsfit
aftfit <- fit@aftfit
uvfit <- fit@uvfit
grpnames <- settings$grpnames
offset <- settings$addconstant
if(class(aftfit)[1]=="survreg") {
aft <- TRUE
validArgModel(...)
}
else if(class(uvfit)[1]=="gls") {
uv <- TRUE
validArgModel(...)
}
if(class(rrfit)[1]=="rlm" && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(class(olsfit)[1]=="lm" && model!="rronly" && !aft && !uv) {
ols <- TRUE
if(!rr) model <- "olsonly"
}
n <- with(olsfit$dfru, sapply(split(endpt, grpf), length))
if(rr || ols) {
df.residual <- olsfit$df.residual
}
if(ols) {
olsestimates <- olsfit$coef
varcov.olsfit <- vcov(olsfit)
ols.grps <- grpsummary(olsestimates,
varcov.olsfit,
n=n, errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
offset=offset, display="none", ...)
if(mcadjust) { multcompDone("Classical Least Squares") }
}
if(rr) {
rrestimates <- rrfit$coef
summ.rrfit <- summary(rrfit, method="XtWX")
stddev <- summ.rrfit$stddev
varcov.rrfit <- (stddev^2) * summ.rrfit$cov.unscaled
rr.grps <- grpsummary(rrestimates,
varcov.rrfit,
n=n, errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
offset=offset, display="none", ...)
if(mcadjust) { multcompDone("Resistant & Robust") }
}
if(aft) {
aftestimates <- aftfit$coef
df.residual <- aftfit$df.residual
varcov.aftfit <- chop.matrix(vcov(aftfit))
aft.grps <- grpsummary(aftestimates,
varcov.aftfit,
n=n, errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
offset=offset, display="none",...)
}
if(uv) {
if(mcadjust) {
stop(cgMessage("The mcadjust argument cannot be set",
"to TRUE when the fitted model allowed",
"for unequal variances."))
}
uv.grps <- grpsummary(uvfit$coef,
vcov(uvfit),
n=n, errordf=n-1,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
offset=offset, display="none",...)
}
settings <- list(endptscale=settings$endptscale,
analysisname=settings$analysisname,
endptname=settings$endptname,
alpha=alpha,
mcadjust=mcadjust,
digits=settings$digits,
sandaft=if(aft) TRUE else FALSE)
returnObj <- new("cgOneFactorGrpSummaryTable",
ols.grps=ols.grps, rr.grps=rr.grps,
aft.grps=aft.grps, uv.grps=uv.grps,
settings=settings)
if(display=="print") {
print(returnObj)
}
else if(display=="show") {
show(returnObj)
}
## else display=="none"
invisible(returnObj)
})
setMethod("print", "cgOneFactorGrpSummaryTable",
print.cgOneFactorGrpSummaryTable <-
function(x, digits=NULL, title=NULL, endptname=NULL, ...) {
##
## PURPOSE: Semi-formatted print version of Comparisons Table
##
## NOTE: Had to use x as an argument because of the system defined
## generic. I would have preferred to use object; hence the first
## statement below.
##
object <- x
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
olsgrps <- object@ols.grps
rrgrps <- object@rr.grps
aftgrps <- object@aft.grps
uvgrps <- object@uv.grps
ols <- rr <- uv <- aft <- FALSE ## Initializations
if(!is.null(aftgrps)) {
aft <- TRUE
validArgModel(...)
}
else if(!is.null(uvgrps)) {
uv <- TRUE
validArgModel(...)
}
if(!is.null(rrgrps) && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(!is.null(olsgrps) && (model!="rronly") && !aft && !uv) {
ols <- TRUE
}
settings <- object@settings
alpha <- settings$alpha
mcadjust <- settings$mcadjust
if(is.null(digits)) {
digits <- settings$digits
}
else {
digits <- validArgDigits(digits)
}
curscipen <- getOption("scipen")
on.exit(options(scipen=curscipen), add=TRUE)
options(scipen=9)
if(is.null(title)) {
title <- paste("Group Summary Table of", settings$analysisname)
}
else {
validCharacter(title)
}
if(is.null(endptname)) {
endptname <- settings$endptname
if(!is.character(endptname)) {
endptname <- ""
}
}
else {
validCharacter(endptname)
}
cat(title,"\n")
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
metric <- "Means"
if(settings$endptscale=="log") {
metric <- paste("Geometric", metric)
}
else {
metric <- paste("Arithmetic", metric)
}
cat(metric, "\n")
fmtdig <- function(x, digits) {
x$estimate <- fround(x$estimate, digits)
x$se <- fround(x$se, digits)
x$lowerci <- fround(x$lowerci, digits)
x$upperci <- fround(x$upperci, digits)
x$n <- fround(x$n, 0)
return(x)
}
informConfidence <- function() {
cat(paste(round(100*(1-alpha), 0), "% Confidence ",
"(alpha of ", alpha, ")",
if(mcadjust) {", Multiplicity Adjusted"},
"\n",
sep=""))
}
curwidth <- getOption("width")
on.exit(options(width=curwidth), add=TRUE)
if(curwidth < 500) { options(width=500) }
if(ols) {
cat("\nClassical Least Squares Model Fit\n")
informConfidence()
print(fmtdig(olsgrps, digits), quote=FALSE)
}
if(rr) {
cat("\nResistant & Robust Model Fit\n")
informConfidence()
print(fmtdig(rrgrps, digits), quote=FALSE)
}
if(aft) {
afttitle <- "\nAccelerated Failure Time Model Fit\n"
if(settings$sandaft) {
cat(paste(afttitle,
"with Sandwich Variance-Covariance Estimate\n",
sep=""))
}
else {
cat(afttitle)
}
informConfidence()
print(fmtdig(aftgrps, digits), quote=FALSE)
}
if(uv) {
cat("\nUnequal Variances Model Fit\n")
informConfidence()
print(fmtdig(uvgrps, digits), quote=FALSE)
}
invisible()
})
setMethod("show", "cgOneFactorGrpSummaryTable",
show.cgOneFactorGrpSummaryTable <- function(object) showDefault(object))
comparisonsgraph <- function(compstable,
difftype,
analysisname="",
endptname="",
alpha=0.05,
digits=NULL,
titlestamp=TRUE,
explanation=TRUE,
wraplength=20,
cex.comps=0.7,
ticklabels=NULL,
...) {
##
## PURPOSE: Plot Confidence Intervals of Comparisons;
## use horizontal error bars
##
## compstable needs to be a dataframe of proper format
difftype <- validArgMatch(difftype, c("percent","amount","simple"))
validAlpha(alpha)
validBoolean(titlestamp)
validBoolean(explanation)
digits <- validArgDigits(digits)
options(warn=-1)
curpar <- par(new=FALSE, mgp=c(3,0.25,0), tck=-0.010, mar=c(5,7,4,2)+0.1)
options(warn=0)
on.exit(par(curpar), add=TRUE)
numberofcomprs <- nrow(compstable)
comprnames <- row.names(compstable)
est <- compstable[, "estimate"]
lower <- compstable[, "lowerci"]
upper <- compstable[, "upperci"]
compys <- numberofcomprs:1
if(difftype=="percent") {
est <- log(pctToRatio(est))
lower <- log(pctToRatio(lower))
upper <- log(pctToRatio(upper))
}
allx <- c(lower, upper, est)
## We include zero because we always want the
## "no-difference" reference line to be visible
xmin <- min(c(0, allx))
xmax <- max(c(0, allx))
if(difftype=="percent") {
## snippet from Hmisc's errbar function
plot(est/log(10), compys, ylim = c(0.5, numberofcomprs + 0.5),
xlim = c(xmin, xmax)/log(10),
ylab="",
xlab = "",
pch = 16, axes=FALSE, ...)
mtext(side=1, text="log-spaced", line=4, cex=0.7)
ycoord <- par()$usr[3:4]
segments(lower/log(10), compys, upper/log(10), compys)
smidge <- 0.015 * (ycoord[2] - ycoord[1])/2
segments(lower/log(10), compys - smidge, lower/log(10), compys + smidge)
segments(upper/log(10), compys - smidge, upper/log(10), compys + smidge)
## zero reference line
abline(v=0, lty=3)
## Axes Customization
log10endpt <- log10(exp(allx))
xratioticks <- setupAxisTicks(exp(allx), ratio=TRUE, percent=TRUE,
axis="x", digits=digits)
if(!is.null(ticklabels)) {
xratioticks <- makeTickMarks(ticklabels, xratioticks,
percent=TRUE)
}
axis(1,at=log10(xratioticks),
labels=names(xratioticks),
tck=-.010, cex.axis=0.8)
xlabchar <- "Percent Difference"
}
else { ## difftype = "amount" or "simple" (Simple Difference)
plot(est, compys, ylim = c(0.5, numberofcomprs + 0.5),
xlim = c(xmin, xmax),
ylab="",
xlab="",## paste("Difference in\n ", endptname, sep=""),
pch=16, axes=FALSE, ...)
ycoord <- par()$usr[3:4]
segments(lower, compys, upper, compys)
smidge <- 0.015 * (ycoord[2] - ycoord[1])/2
segments(lower, compys - smidge, lower, compys + smidge)
segments(upper, compys - smidge, upper, compys + smidge)
## zero reference line
abline(v=0, lty=3)
## Axes Customization
xdiffticks <- setupAxisTicks((allx), ratio=FALSE, difference=TRUE,
logscale=FALSE,
axis="x", digits=digits)
if(!is.null(ticklabels)) {
xdiffticks <- makeTickMarks(ticklabels, xdiffticks,
percent=FALSE)
}
axis(1, at=xdiffticks,
labels=names(xdiffticks),
tck=-.010, cex=0.8)
xlabchar <- "Difference"
}
mtext(side=1, text=xlabchar, line=2)
if(is.expression(endptname)) {
mtext(side=1, text=catCharExpr("in ", endptname), line=3)
}
else {
mtext(side=1, text=paste("in", endptname), line=3)
}
## Wrap long comparison labels
for(i in seq(along=comprnames)) {
##if(nchar(comprnames[i]) > wraplength) {
## comprnames[i] <- gsub("vs.","vs.\n", comprnames[i])
##}
if(nchar(comprnames[i]) > wraplength) {
comprnames[i] <- paste(strwrap(comprnames[i], width=wraplength,
exdent=1), collapse="\n")
}
}
axis(2, at=compys, labels=comprnames,
adj=1, cex.axis=cex.comps, las=1)
box()
## Titles
par(curpar)
if(explanation) {
comparisonsGraphStamp(alphapercent=100*alpha)
}
if(titlestamp) {
title(main=paste("Comparisons Graph\n",
analysisname, sep=""), line=2, cex.main=1.1)
}
## min-max annotation
bardata <- c(est, lower, upper)
minmaxTicks(if(difftype=="percent") exp(bardata) else bardata,
theaxis="x", percent=if(difftype=="percent") TRUE else FALSE,
ratio=if(difftype=="percent") TRUE else FALSE,
logscale=if(difftype=="percent") TRUE else FALSE,
digits=digits)
## Add 0
text(y=par("usr")[3], x=0,
labels="0\n", col="blue", adj=0, cex=0.7)
invisible()
}
setMethod("comparisonsGraph", "cgOneFactorComparisonsTable",
comparisonsGraph.cgOneFactorComparisonsTable <-
function(compstable, cgtheme=TRUE, device="single",
wraplength=20, cex.comps=0.7, ...) {
##
## PURPOSE: create a graph of CI's on differences
## to judge comparisons
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names=c("model", "ticklabels"))
validBoolean(cgtheme)
settings <- compstable@settings
difftype <- if(settings$endptscale=="log") "percent" else "simple"
alpha <- settings$alpha
alphapercent <- round(100*alpha, 0)
mcadjust <- settings$mcadjust
analysisname <- settings$analysisname
endptlabel <- settings$endptlabel
stamps <- settings$stamps
digits <- settings$digits
rrcomprs <- compstable@rr.comprs
olscomprs <- compstable@ols.comprs
aftcomprs <- compstable@aft.comprs
uvcomprs <- compstable@uv.comprs
ticklabelsarg <- getDotsArgName(dots, "ticklabels")
if(!is.na(ticklabelsarg)) {
ticklabels <- eval(parse(text=paste("dots$", ticklabelsarg, sep="")))
validList(ticklabels, names=c("mod","marks"),
argname="ticklabels")
}
else {
ticklabels <- NULL
}
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
dots$model <- "both"
model <- "both"
}
aft <- ols <- rr <- uv <- FALSE ## Initializations
if(!is.null(aftcomprs)) {
aft <- TRUE
validArgModel(...)
}
else if(!is.null(uvcomprs)) {
uv <- TRUE
validArgModel(...)
}
if(!is.null(rrcomprs) && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(!is.null(olscomprs) && (model!="rronly") && !aft && !uv) {
ols <- TRUE
}
device <- validArgMatch(device, c("single","multiple", "ask"))
thetitle <- "Comparisons Graph"
if(rr && ols && is.element(model, "both") && device=="single") {
all.dfr <- rbind(olscomprs, rrcomprs)
all.dfr$typef <- factorInSeq(c(rep("Classical",
nrow(olscomprs)),
rep("Resistant & Robust",
nrow(rrcomprs))))
thetitle <- paste(thetitle, "s", sep="")
numberofcomprs <- nrow(olscomprs)
comprnames <- row.names(olscomprs)
cgDevice(cgtheme=cgtheme)
trellispanelstg <- trellis.par.get("clip")$panel
trellis.par.set("clip", list(panel="off"))
on.exit(trellis.par.set("clip", list(panel=trellispanelstg)),
add=TRUE)
trellisparstg2 <- trellis.par.get("layout.widths")
ymargin <- max(nchar(comprnames))/2
trellis.par.set("layout.widths", list(left.padding=ymargin,
ylab=ymargin))
on.exit(trellis.par.set("layout.widths", trellisparstg2),
add=TRUE)
trellisparstg3 <- trellis.par.get("axis.components")
trellis.par.set("axis.components",
list(left=list(pad1=0.2, pad2=2),
bottom=list(pad1=0.2, pad2=1),
top=list(pad1=1, pad2=1),
right=list(pad1=0.5, pad2=1)
))
on.exit(trellis.par.set("axis.components", trellisparstg3),
add=TRUE)
allx <- unlist(all.dfr[, c("estimate","lowerci","upperci")])
est <- all.dfr[, "estimate"]
lower <- all.dfr[, "lowerci"]
upper <- all.dfr[, "upperci"]
compys <- numberofcomprs:1
if(difftype=="percent") {
est <- log(pctToRatio(est))
lower <- log(pctToRatio(lower))
upper <- log(pctToRatio(upper))
allx <- log(pctToRatio(allx))
xmin <- min(c(0, allx))
xmax <- max(c(0, allx))
## ideas in the next call adapted from errbar function
thegraph <- xyplot(compys ~ Cbind(est,
lower,
upper)
| typef,
data=all.dfr,
digits=digits,
panel = function(x, y, ...) {
xlower <- attr(x, "other")[, 1]
xupper <- attr(x, "other")[, 2]
xcex <- 0.6
xratioticks <-
setupAxisTicks(exp(c(xlower, x, xupper)),
ratio=TRUE, percent=TRUE,
axis="x",
grid=TRUE, xcex=xcex,
digits=digits)
if(!is.null(ticklabels)) {
xratioticks <-
makeTickMarks(ticklabels,
xratioticks,
percent=TRUE)
}
panel.axis(side="bottom",
at=log(xratioticks),
labels=names(xratioticks),
text.cex=xcex,
tck=0.20, outside=TRUE, rot=0)
panel.points(x, compys, pch=16, col="black")
lsegments(xlower, compys, xupper, compys)
##ycoord <- par()$usr[3:4]
##smidge <- 0.015 * (ycoord[2] - ycoord[1])
smidge <- 0.015 * max(compys)
lsegments(xlower, compys - smidge,
xlower, compys + smidge)
lsegments(xupper, compys - smidge,
xupper, compys + smidge)
## zero reference line
panel.abline(v=0, lty=3)
## Add 0 if not present
if(all(xratioticks!=1)) {
panel.axis(side="bottom",
at=0, labels="0", text.cex=0.6,
tck=0.20, outside=TRUE)
}
if(panel.number()==1) {
## Wrap long comparison labels
for(i in seq(along=comprnames)) {
if(nchar(comprnames[i]) > wraplength) {
comprnames[i] <-
paste(strwrap(comprnames[i],
width=wraplength,
exdent=1),
collapse="\n")
}
}
panel.axis(side="left",
at=compys,
labels=comprnames,
tck=0.20, text.cex=cex.comps,
rot=0,
outside=TRUE)
}
xrange <- range(c(xlower, xupper))
panel.text(y=rep(0.5, 2), x=xrange,
labels=paste(
c(fmtRatioToPercent(exp(min(xrange))),
fmtRatioToPercent(exp(max(xrange)))), "\n"),
col="blue", adj=0, cex=0.6)
},
layout=c(2,1), aspect=1, pch=16,
col="black",
as.table=TRUE,
ylim = c(0.5, numberofcomprs + 0.5),
xlim=rangeExtend(c(xmin, xmax), pct=10),
xlab="Percent Difference",
ylab="",
scales=list(y=list(at=1:numberofcomprs,
labels=NULL, tck=c(0.15, 0), axs="i"),
x=list(labels=NULL, tck=0)),
main=list(label=paste(thetitle, "\n",
analysisname, sep=""), cex=1.1),
par.strip.text=list(cex=0.7)
)
}
else {
xmin <- min(c(0, allx))
xmax <- max(c(0, allx))
## ideas in the next call adapted from Hmisc::errbar function
thegraph <- xyplot(compys ~ Cbind(est,
lower,
upper)
| typef,
data=all.dfr,
digits=digits,
panel = function(x, y, ...) {
xlower <- attr(x, "other")[, 1]
xupper <- attr(x, "other")[, 2]
xcex <- 0.6
xdiffticks <-
setupAxisTicks(c(xlower, x, xupper),
difference=TRUE,
axis="x",
logscale=FALSE,
grid=TRUE,
xcex=0.6,
digits=digits)
if(!is.null(ticklabels)) {
xdiffticks <-
makeTickMarks(ticklabels,
xdiffticks,
percent=FALSE)
}
panel.axis(side="bottom",
at=xdiffticks,
labels=names(xdiffticks),
text.cex=xcex,
tck=0.20, outside=TRUE, rot=0)
panel.points(x, compys, pch=16, col="black")
lsegments(xlower, compys, xupper, compys)
## ycoord <- par()$usr[3:4]
## smidge <- 0.015 * (ycoord[2] - ycoord[1])
smidge <- 0.015 * max(compys)
lsegments(xlower, compys - smidge,
xlower, compys + smidge)
lsegments(xupper, compys - smidge,
xupper, compys + smidge)
## zero reference line
panel.abline(v=0, lty=3)
## Add 0 if not present
if(all(xdiffticks!=0)) {
panel.axis(side="bottom",
at=0, labels="0", text.cex=0.6,
tck=0.20, outside=TRUE)
}
if(panel.number()==1) {
## Wrap long comparison labels
for(i in seq(along=comprnames)) {
if(nchar(comprnames[i]) > wraplength) {
comprnames[i] <-
paste(strwrap(comprnames[i],
width=wraplength,
exdent=1),
collapse="\n")
}
}
panel.axis(side="left",
at=compys,
labels=comprnames,
tck=0.20, text.cex=cex.comps,
rot=0,
outside=TRUE)
}
xrange <- range(c(xlower, xupper))
panel.text(y=rep(0.5, 2), x=xrange,
labels=paste(
c(fmtDifference(min(xrange),
digits=digits),
fmtDifference(max(xrange),
digits=digits)), "\n"),
col="blue", adj=0, cex=0.6)
},
layout=c(2,1), aspect=1, pch=16,
col="black",
as.table=TRUE,
ylim=c(0.5, numberofcomprs + 0.5),
xlim=rangeExtend(c(xmin, xmax), pct=10),
xlab="Difference",
ylab="",
scales=list(y=list(at=1:numberofcomprs,
labels=NULL, tck=c(0.15, 0), axs="i"),
x=list(labels=NULL, tck=0)),
main=list(label=paste(thetitle, "\n",
analysisname, sep=""), cex=1.1),
par.strip.text=list(cex=0.7)
)
}
print(thegraph)
seekViewport(trellis.vpname("xlab"))
grid.text(catCharExpr("in ", endptlabel), y = unit(-1, "lines"))
upViewport(0)
comparisonsGraphStamp(mcadjust, alphapercent, grid=TRUE,
desc=settings$type)
if(stamps) graphStampCG()
}
else if((model=="olsonly" || (ols && !rr && model=="both")) &&
device=="single") {
comparisonsgraph(olscomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Classical analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else if((model=="rronly" || (!ols && rr && model=="both")) &&
device=="single") {
## else if(model=="rronly" & rr & device=="single") {
comparisonsgraph(rrcomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Resistant & Robust analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else if(rr && ols &&
is.element(device, c("ask","multiple")) &&
is.element(model, "both")) {
device <- validArgMatch(device, c("multiple", "ask"))
if(device=="ask") {
op <- par(ask = TRUE)
on.exit(par(op), add=TRUE)
}
comparisonsgraph(olscomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Classical analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
if(device=="multiple") {
## dots$model <- NULL ## avoid unused argument error
eval(parse(text=paste("dots$", modelarg, " <- NULL", sep="")))
do.call("cgDevice", c(list(new=TRUE), dots))
cat(cgMessage("A new graphics device has been generated",
"to hold the Resistant & Robust",
"Comparisons Graph version.",
"The Classical Least Squares version is on the previous",
"device.\n",
warning=TRUE))
}
comparisonsgraph(rrcomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Resistant & Robust analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else if(aft && device=="single") {
comparisonsgraph(aftcomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
explanation=FALSE, ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, AFT analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else if(uv && device=="single") {
comparisonsgraph(uvcomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Unequal Variances analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else {
stop(cgMessage("The chosen device and model arguments",
"are not compatible either with each other",
"or with the fitted model(s) in the fit object.",
seeHelpFile("comparisonsGraph")))
}
invisible()
}
)
validErrorDf <- function(errordf, varcovmatrix,
n, mcadjust) {
##
if(errordf=="approx") {
## First check that estimates are uncorrelated, i.e
## the varcovmatrix is diagonal.
if(sum(varcovmatrix[lower.tri(varcovmatrix)],
varcovmatrix[upper.tri(varcovmatrix)])!=0) {
stop(cgMessage("The varcovmatrix argument appears to have",
"nonzero off-diagonal elements. Currently,",
"the \"approx\" setting for the errordf",
"argument will only work for uncorrelated estimates."))
}
if(!isAllEqual(c(length(diag(varcovmatrix)),
length(n)))) {
stop(cgMessage("The lengths of the n argument and",
"the number of diagonal elements in the varcovmatrix",
"argument appear to differ. They need to be the same."))
}
if(mcadjust) {
stop(cgMessage("The mcadjust argument must be FALSE when the",
"errordf argument is \"approx\"."))
}
}
else if(length(errordf)!=1 || !is.numeric(errordf) || errordf < 1) {
stop(cgMessage("The errordf argument needs to be length 1",
"and set to be 1 or greater or the value \"approx\"."))
}
else return(TRUE)
}
validComparisonType <- function(type, errordf) {
type <- try(match.arg(type, c("pairwisereflect","pairwise",
"allgroupstocontrol", "custom")))
if(class(type)=="try-error") {
stop(cgMessage("The type argument must be one of \"pairwisereflect\",",
"\"pairwise\", \"allgroupstocontrol\",",
"or \"custom\"."))
}
if(errordf=="approx" && !any(is.element(type,
c("pairwisereflect","pairwise",
"allgroupstocontrol")))) {
stop(cgMessage("The type argument must be one of \"pairwisereflect\",",
"\"pairwise\", or \"allgroupstocontrol\"",
"when the errordf argument is \"approx\"."))
}
else if(!any(is.element(type, c("pairwisereflect","pairwise",
"allgroupstocontrol", "custom")))) {
stop(cgMessage("The type argument must be one of \"pairwisereflect\",",
"\"pairwise\", \"allgroupstocontrol\",",
"or \"custom\"."))
}
else return(TRUE)
}
validEstimates <- function(x) {
if(!is.numeric(x) || length(x) < 2 ) {
stop(cgMessage("The estimates argument needs to be numeric",
"and be of length two or more."
))
}
if(is.null(thenames <- names(x))) {
stop(cgMessage("The estimates arguments needs to have",
"a names attribute to use for",
"thes group labels."))
}
if(any(is.na(thenames)) || any(thenames=="")) {
stop(cgMessage("One or more of the groups names in the",
"names attribute of estimates",
"appears to be missing."))
}
else return(x)
}
validGrpNames <- function(x) {
if(!is.character(x) || length(x) < 2 ) {
stop(cgMessage("The grpnames argument needs to be character",
"and be of length two or more.",
))
}
if(any(is.na(x)) || any(x=="")) {
stop(cgMessage("One or more of the groups names in the",
"names attribute of estimates",
"appears to be missing.",
))
}
else return(x)
}
validN <- function(n, estimates) {
if(!is.numeric(n) || any(n < 1)) {
stop(cgMessage("The n argument vector needs to be numeric",
"and all elements must be one or greater."))
}
else if(length(n)!=length(estimates)) {
stop(cgMessage("The n argument vector needs to be of the",
"same length as the estimates vector."))
}
else return(as.integer(round(n, 0)))
}
validCnames <- function(type, x, reqlength) {
if(length(x)==1 && x=="derive") return(TRUE)
if(type!="custom") {
stop(cgMessage("The cnames argument must be left at its",
"default value of \"derive\" when type is",
"not \"custom\"."))
}
if(length(x) != reqlength) {
stop(cgMessage("The cnames argument must be either the",
"value \"derive\", ",
"or must be of length", reqlength,
", the number of rows of the",
"contrastmatrix argument."))
}
return(TRUE)
}
validAddPct <- function(addpct, endptscale) {
if(addpct && endptscale=="log") {
stop(cgMessage("The addpct argument cannot be TRUE,",
"and must be set to FALSE,",
"when the \"endptscale\" argument is",
"\"log\"."))
}
return(TRUE)
}
Try the cg package in your browser
Any scripts or data that you put into this service are public.
cg documentation built on May 2, 2019, 9:26 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions comparisons pairwisecompsmatrix errorbargraph grpsummary comparisonsgraph validErrorDf validComparisonType validEstimates validGrpNames validN validCnames validAddPct
Documented in comparisons comparisonsgraph errorbargraph grpsummary pairwisecompsmatrix validComparisonType validErrorDf validEstimates validGrpNames validN
## $Id: p03FormalOneFactorData.R 6053 2015-02-22 20:23:45Z bpikouni $
## One-Factor Unpaired Groups Case
## Formal Analysis methods for One-Factor Unpaired Groups Data
setClass("cgOneFactorGlobalTest",
representation(ols.gpval="numericOrNULL", rr.gpval="numericOrNULL",
aft.gpval="numericOrNULL", uv.gpval="numericOrNULL",
settings="list"),
prototype(ols.gpval=NULL, rr.gpval=NULL, aft.gpval=NULL,
uv.gpval=NULL,
settings=list()))
setMethod("globalTest", "cgOneFactorFit",
globalTest.cgOneFactorFit <-
function(fit, display="print", ...) {
##
## PURPOSE: Derive global F-test p-values
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
## initializations
aft <- ols <- rr <- uv <- FALSE
rr.gpval <- ols.gpval <- aft.gpval <- uv.gpval <- NULL
##
rrfit <- fit@rrfit
olsfit <- fit@olsfit
aftfit <- fit@aftfit
uvfit <- fit@uvfit
settings <- fit@settings
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
if(class(aftfit)[1]=="survreg") {
aft <- TRUE
validArgModel(...)
}
else if(class(uvfit)[1]=="gls") {
uv <- TRUE
validArgModel(...)
}
if(class(rrfit)[1]=="rlm" && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(class(olsfit)[1]=="lm" && model!="rronly" && !aft && !uv) {
ols <- TRUE
if(!rr) model <- "olsonly"
}
if(rr || ols) {
df.residual <- olsfit$df.residual
df.grpf <- length(settings$grpnames) - 1
upform <- update(formula(olsfit), . ~ grpf)
}
dfru <- olsfit$dfru
if(ols) {
## calculations mimic print.summary.lm code snippet
olsfit <- lm(upform, data=dfru)
ols.gpval <- with(summary(olsfit), pf(fstatistic[1],
fstatistic[2],
fstatistic[3], lower.tail
= FALSE))
}
if(rr) {
## if only 2 groups, we rely on the comparisonstable method for computing
## the pvalue
if(df.grpf==1) {
rr.gpval <- pairwisecompsmatrix(comparisonsTable(fit,
type="pairwisereflect")@rr.comprs,
grpnames=settings$grpnames)[3, 3]
}
else {
## As of current writing,
## The MASS summary.rlm does not compute or facilitate construction of a
## global F, or equivalently, an R^2. For now we use a placeholder method
## based on ad-hoc (witchcraft) development, which is essentially to
## re-fit a linear model with lm() and weights from the RR fit.
olsw.rrfit <- lm(upform, data=dfru, weights=rrfit$w)
rr.Fobs <- { (1 / (1 / summary.lm(olsw.rrfit)$r.squared - 1))*
(df.residual/df.grpf) }
rr.gpval <- 1 - pf(rr.Fobs, df.grpf, df.residual)
}
}
else if(aft) {
upform <- update(formula(aftfit), . ~ grpf)
aftfit <- survreg(upform, data = dfru,
dist = "gaussian", maxiter = aftfit$maxIter)
## calculations mimic print.survreg code snippet
chisq <- 2 * diff(aftfit$loglik)
df <- sum(aftfit$df) - aftfit$idf
aft.gpval <- 1 - pchisq(chisq, df)
}
else if(uv) {
upform <- update(formula(uvfit), . ~ grpf)
uvfit <- gls(model = upform, data = dfru,
weights = varIdent(form = ~1 | grpf))
uv.gpval <- anova(uvfit)$"p-value"[2]
}
x <- new("cgOneFactorGlobalTest",
ols.gpval=ols.gpval, rr.gpval=rr.gpval,
aft.gpval=aft.gpval, uv.gpval=uv.gpval,
settings=settings)
if(display=="print") {
print(x)
}
else if (display=="show"){
showDefault(x)
}
## else show nothing
invisible(x)
})
setMethod("print", "cgOneFactorGlobalTest",
print.cgOneFactorGlobalTest <-
function(x, title=NULL, endptname=NULL, ...) {
##
## PURPOSE: Semi-formatted print version of Global test p-value
## Three digits are used, and anything smaller than or equal to
## 0.0005 gets "< 0.001"
##
## NOTE: Had to use x as an argument because of the system defined
## generic. I would have preferred to use object; hence the first
## statement below.
##
object <- x
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
olsgrps <- object@ols.gpval
rrgrps <- object@rr.gpval
aftgrps <- object@aft.gpval
uvgrps <- object@uv.gpval
ols <- rr <- uv <- aft <- FALSE ## Initializations
if(!is.null(aftgrps)) {
aft <- TRUE
validArgModel(...)
}
else if(!is.null(uvgrps)) {
uv <- TRUE
validArgModel(...)
}
if(!is.null(rrgrps) && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(!is.null(olsgrps) && (model!="rronly") && !aft && !uv) {
ols <- TRUE
}
settings <- object@settings
alpha <- settings$alpha
if(is.null(title)) {
title <- paste("Global Test P-value of", settings$analysisname)
}
else {
validCharacter(title)
}
if(is.null(endptname)) {
endptname <- settings$endptname
if(!is.character(endptname)) {
endptname <- ""
}
}
else {
validCharacter(endptname)
}
cat(title,"\n")
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
if(ols) {
cat("\nClassical Least Squares Model Fit:", fmtPvalue(olsgrps), "\n")
}
if(rr) {
cat("\nResistant & Robust Model Fit:", fmtPvalue(rrgrps), "\n")
}
if(aft) {
cat("\nAccelerated Failure Time Model Fit:", fmtPvalue(aftgrps), "\n")
}
if(uv) {
cat("\nUnequal Variances Model Fit:", fmtPvalue(uvgrps), "\n")
}
invisible()
})
setMethod("show", "cgOneFactorGlobalTest",
show.cgOneFactorGlobalTest <- function(object) showDefault(object))
## Comparisons Tables
comparisons <- function(estimates,
varcovmatrix, errordf=Inf,
endptscale, mcadjust=FALSE,
alpha=0.05,
type="pairwisereflect",
contrastmatrix=NULL, n,
offset=NULL, cnames="derive",
analysisname="", endptname="",
digits=NULL,
addpct=FALSE,
display="print") {
##
## PURPOSE: Function for computing comparisons results
## and placing in a listed form.
##
## Input argument checking
if(!hasArg(estimates)) reportInvalidArg("estimates")
if(!hasArg(varcovmatrix)) reportInvalidArg("varcovmatrix")
validAlpha(alpha)
validBoolean(mcadjust)
validBoolean(addpct)
validErrorDf(errordf, varcovmatrix, n, mcadjust)
validComparisonType(type, errordf)
validCnames(type, cnames, nrow(contrastmatrix))
estimates <- validEstimates(estimates)
endptscale <- validArgMatch(endptscale, c("log","original"))
validAddPct(addpct, endptscale)
display <- validArgMatch(display, c("print","none","show"))
ngrps <- length(estimates)
grpnames <- names(estimates)
if(is.null(contrastmatrix)) {
L <- contrastMatrix(list(grpnames), type=type)
## if type="allgroupstocontrol", then the first grpname is set to be
## the reference group level.
}
else if(type=="custom") {
L <- contrastmatrix
}
## If only two groups were analyzed, then multiplicity adjustment has
## no effect. So we ensure mcadjust is set to FALSE
if(ngrps==2) { mcadjust <- FALSE }
if(mcadjust) {
set.seed(17)
multcompInform()
glhtobj <- list(object=NULL, linfct=L, coef=estimates,
vcov=varcovmatrix,
type="user-defined", alternative="two.sided",
df=errordf,
rhs=rep(0, nrow(L)))
class(glhtobj) <- "glht"
mcpci.fit <- confint(glhtobj, level=1-alpha)
mcptest.fit <- summary(glhtobj)
estdiff <- mcpci.fit$confint[,"Estimate"]
sediff <- mcptest.fit$test$sigma
lowerci <- mcpci.fit$confint[,"lwr"]
upperci <- mcpci.fit$confint[,"upr"]
pval <- mcptest.fit$test$pvalues
}
else if(!mcadjust) {
estdiff <- L %*% estimates
sediff <- sqrt(diag(L %*% varcovmatrix %*% t(L)))
if(errordf=="approx") {
## Apply Satterthwaite approximation for special case of pairwise
## or all groups to control comparisons (-1 vs 1 contrasts)
incidence <- t(apply(L, 1, function(x) { x!=0 }))
thestnderrs <- sqrt(diag(varcovmatrix))
errordf <- apply(incidence, 1,
function(x, thestnderrs, n) {
sepair <- thestnderrs[x]
npair <- n[x]
sese <- sqrt(sepair[1]^2 + sepair[2]^2)
return(sese^4/(sepair[1]^4/(npair[1] - 1) +
sepair[2]^4/(npair[2] - 1)))
}, thestnderrs=thestnderrs, n=n)
}
tcrit <- qt(1 - alpha/2, errordf)
lowerci <- estdiff - tcrit * sediff
upperci <- estdiff + tcrit * sediff
pval <- 2 * (1 - pt(abs(estdiff/sediff), errordf))
}
mcp <- data.frame(estimate=estdiff, se=sediff,
lowerci=lowerci, upperci=upperci,
pval=pval)
## Log Case:
## Transformations to Percent Difference Expression Scale
## And we need to specially handle the asymmetry of percent change,
## since decreases cannot be greater than 100%
## also note that if offset is a valid number,
## some calculations will be replaced.
if(endptscale=="log") {
## Focus only on the percent differences for now
logscalest <- mcp[,1]
if(is.null(offset)) {
mcp[, 1] <- 100 * ( exp(logscalest) - 1 ) # Point Estimate
mcp[, 2] <- 100 * ( exp(logscalest) * mcp[, 2] ) # Std Err
mcp[, 3:4] <- 100 * ( exp(mcp[, 3:4]) - 1) # Confidence Limits
}
else { ## offset
validNumeric(offset, positive=TRUE)
mA <- exp(mcp$meanA) - offset
mB <- exp(mcp$meanB) - offset
mAmBratio <- mA/mB
correctionFactor <- (mA * (mB + offset))/(mB * (mA + offset))
mcp[, 1] <- 100 * ( mAmBratio - 1 ) # Point Estimate
mcp[, 2] <- 100 * ( mAmBratio * mcp[, 2] ) # Std Err
mcp[, 3:4] <- 100 * ( correctionFactor*exp(mcp[, 3:4]) - 1) # Confidence Limits
}
}
## Fetch, arrange, and add individual component estimates
if(type!="custom") {
complist <- strsplit(row.names(L), split=" vs. ", fixed=TRUE)
compA <- sapply(complist, function(x) x[1])
compB <- sapply(complist, function(x) x[2])
compAlist <- as.list(compA)
compBlist <- as.list(compB)
compAindx <- unlist(sapply(compAlist,
function(x) which(grpnames==x)))
compBindx <- unlist(sapply(compBlist,
function(x) which(grpnames==x)))
sevec <- sqrt(diag(varcovmatrix))
mcp$meanA <- estimates[compAindx]
mcp$seA <- sevec[compAindx]
mcp$meanB <- estimates[compBindx]
mcp$seB <- sevec[compBindx]
}
else if(type=="custom") {
L.A <- t(apply(L, 1, function(x) ifelse(x > 0, x, 0)))
## For the B side of the contrast, we need to remove the minus sign
L.B <- t(apply(L, 1, function(x) ifelse(x < 0, -x, 0)))
mcp$meanA <- L.A %*% estimates
mcp$seA <- sqrt(diag(L.A %*% varcovmatrix %*% t(L.A)))
mcp$meanB <- L.B %*% estimates
mcp$seB <- sqrt(diag(L.B %*% varcovmatrix %*% t(L.B)))
if(length(cnames)==1 && cnames=="derive") {
names.ests <- names(estimates)
names.A <- apply(L.A, 1,
function(x) {
indx <- (x > 0)
coef.ests <- x[indx]
coef.ests <- ifelse(coef.ests%%1!=0,
signif(coef.ests, 3), coef.ests)
label <- paste(paste(coef.ests, names.ests[indx], sep="*"),
collapse="+")
if(length(coef.ests) > 1) {
label <- paste("(", label, ")", sep="")
}
gsub("1\\*","", label)
})
names.B <- apply(L.B, 1,
function(x) {
indx <- (x > 0)
coef.ests <- x[indx]
coef.ests <- ifelse(coef.ests%%1!=0,
signif(coef.ests, 3), coef.ests)
label <- paste(paste(coef.ests, names.ests[indx], sep="*"),
collapse="+")
if(length(coef.ests) > 1) {
label <- paste("(", label, ")", sep="")
}
gsub("1\\*","", label)
})
rownames(mcp) <- paste(names.A, names.B, sep=" vs. ")
}
}
if(endptscale=="log") {
## Now focus on the individual estimates
if(is.null(offset)) {
mcp$meanA <- exp(mcp$meanA)
mcp$meanB <- exp(mcp$meanB)
}
else { ## offset
mcp$meanA <- exp(mcp$meanA) - offset
mcp$meanB <- exp(mcp$meanB) - offset
}
mcp$seA <- mcp$meanA * mcp$seA
mcp$seB <- mcp$meanB * mcp$seB
}
if(cnames[1]!="derive") { row.names(mcp) <- cnames }
if(endptscale=="original" && addpct==TRUE) {
mcp$pctdiff <- 100*(mcp$meanA / mcp$meanB - 1)
}
if(display=="print") {
fmt.mcp <- mcp
cat(paste("Comparisons Table",
if(analysisname!="") paste("for", analysisname), "\n"))
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
diffmetric <- "Differences (A vs. B)"
if(endptscale=="log") {
diffmetric <- paste("Percent", diffmetric)
}
cat(diffmetric, "\n")
pvalfmt <- fmtPvalue(fmt.mcp$pval)
if(regexpr("Percent", diffmetric) > 0) {
fmt.mcp$estimate <- fmtPercent(fmt.mcp$estimate)
fmt.mcp$se <- fmtPercent(fmt.mcp$se)
fmt.mcp$lowerci <- fmtPercent(fmt.mcp$lowerci)
fmt.mcp$upperci <- fmtPercent(fmt.mcp$upperci)
fmt.mcp$pval <- pvalfmt
if(!is.null(digits)) {
fmt.mcp$meanA <- fround(fmt.mcp$meanA, digits=digits)
fmt.mcp$seA <- fround(fmt.mcp$seA, digits=digits)
fmt.mcp$meanB <- fround(fmt.mcp$meanB, digits=digits)
fmt.mcp$seB <- fround(fmt.mcp$seB, digits=digits)
}
{ names(fmt.mcp)[is.element(names(fmt.mcp),
c("meanA","seA",
"meanB","seB"))] <- c("geomeanA","seA",
"geomeanB","seB") }
{ names(mcp)[is.element(names(mcp), c("meanA","seA",
"meanB","seB"))] <-
c("geomeanA","seA",
"geomeanB","seB") }
}
else { ## Simple Differences
pctdiff <- mcp$pctdiff
if(!is.null(digits)) fmt.mcp <- fround(fmt.mcp, digits)
fmt.mcp$pval <- pvalfmt
if(endptscale=="original" && addpct==TRUE) {
fmt.mcp$pctdiff <- fmtPercent(pctdiff)
}
}
cat(paste(round(100*(1-alpha), 0), "% Confidence ",
"(alpha of ", alpha, ")",
if(mcadjust) {", Multiplicity Adjusted"},
"\n",
sep=""))
curwidth <- getOption("width")
on.exit(options(width=curwidth), add=TRUE)
if(curwidth < 500) { options(width=500) }
print(fmt.mcp, quote=FALSE)
}
else if (display=="show"){
showDefault(mcp)
}
## else show nothing
invisible(mcp)
}
setClass("cgOneFactorComparisonsTable",
representation(ols.comprs="dataframeMatrixOrNULL",
rr.comprs="dataframeMatrixOrNULL",
aft.comprs="dataframeMatrixOrNULL",
uv.comprs="dataframeMatrixOrNULL",
settings="list"),
prototype(ols.comprs=NULL,
rr.comprs=NULL,
aft.comprs=NULL,
uv.comprs=NULL,
settings=list()))
setMethod("comparisonsTable", "cgOneFactorFit",
comparisonsTable.cgOneFactorFit <-
function(fit,
## mcadjust=FALSE,
type="pairwisereflect",
## contrastmatrix=NULL,
## refgrp=NULL,
alpha=0.05, addpct=FALSE,
display="print", ...) {
##
## PURPOSE: Produce multiple comparisons of groups based
## on a one-factor unpaired samples fit.
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names=c("mcadjust", "contrastmatrix", "refgrp","model"))
type <- validArgMatch(type, c("pairwisereflect","pairwise",
"allgroupstocontrol", "custom"))
validAlpha(alpha)
validBoolean(addpct)
display <- validArgMatch(display, c("print","none","show"))
mcadjustarg <- getDotsArgName(dots, "mcadjustarg")
if(!is.na(mcadjustarg)) {
mcadjust <- eval(parse(text=paste("dots$", mcadjustarg, sep="")))
validBoolean(mcadjust)
}
else {
mcadjust <- FALSE
}
contrastmatrixarg <- getDotsArgName(dots, "contrastmatrix")
if(is.na(contrastmatrixarg)) {
contrastmatrix <- 1
contrastmatrix <- as.null(contrastmatrix)
}
refgrparg <- getDotsArgName(dots, "refgrp")
settings <- fit@settings
if(!is.na(refgrparg)) {
refgrp <- eval(parse(text=paste("dots$", refgrparg, sep="")))
if(!is.null(refgrp)) {
refgrp <- validArgMatch(refgrp, settings$grpnames)
}
if(type!="allgroupstocontrol") {
stop(cgMessage("If the refgrp= argument is specified",
"then the type= argument needs to be set to",
"\"allgroupstocontrol\"."))
}
}
else refgrp <- settings$refgrp
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
## initializations
aft <- ols <- rr <- uv <- FALSE
aft.comprs <- ols.comprs <- rr.comprs <- uv.comprs <- NULL
##
## settings <- fit@settings
endptscale <- settings$endptscale
rrfit <- fit@rrfit
olsfit <- fit@olsfit
aftfit <- fit@aftfit
uvfit <- fit@uvfit
grpnames <- settings$grpnames
offset <- settings$addconstant
validAddPct(addpct, endptscale)
if(class(aftfit)[1]=="survreg") {
aft <- TRUE
validArgModel(...)
}
else if(class(uvfit)[1]=="gls") {
uv <- TRUE
validArgModel(...)
thens <- with(uvfit$dfru, sapply(split(endpt, grpf), length))
}
if(class(rrfit)[1]=="rlm" && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(class(olsfit)[1]=="lm" && model!="rronly" && !aft && !uv) {
ols <- TRUE
if(!rr) model <- "olsonly"
}
numberofgrps <- length(grpnames)
grpnamesindex <- 1:numberofgrps
if(!is.null(refgrp) && type=="allgroupstocontrol") {
grpnamesindex <- c(which(grpnames==refgrp), which(grpnames!=refgrp))
}
if(rr || ols) {
df.residual <- olsfit$df.residual
}
if(ols) {
olsestimates <- olsfit$coef
varcov.olsfit <- vcov(olsfit)
ols.comprs <- comparisons(olsestimates[grpnamesindex],
varcov.olsfit[grpnamesindex, grpnamesindex],
errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
type=type,
contrastmatrix=contrastmatrix,
offset=offset, addpct=addpct, display="none")
if(mcadjust) { multcompDone("Classical Least Squares") }
}
if(rr) {
rrestimates <- rrfit$coef
summ.rrfit <- summary(rrfit, method="XtWX", ...)
stddev <- summ.rrfit$stddev
varcov.rrfit <- (stddev^2) * summ.rrfit$cov.unscaled
rr.comprs <- comparisons(rrestimates[grpnamesindex],
varcov.rrfit[grpnamesindex, grpnamesindex],
errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
type=type,
contrastmatrix=contrastmatrix,
offset=offset,
addpct=addpct, display="none")
if(mcadjust) { multcompDone("Resistant & Robust") }
}
if(aft) {
aftestimates <- aftfit$coef
## estscaleindex <- length(aftestimates) + 1
df.residual <- aftfit$df.residual
varcov.aftfit <- vcov(aftfit)
aft.comprs <- comparisons(aftestimates[grpnamesindex],
varcov.aftfit[grpnamesindex, grpnamesindex],
errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
type=type,
contrastmatrix=contrastmatrix,
offset=offset, display="none")
if(mcadjust) { multcompDone("Accelerated Failure Time") }
}
if(uv) {
uvestimates <- uvfit$coef
varcov.uvfit <- vcov(uvfit)
uv.comprs <- comparisons(uvestimates[grpnamesindex],
varcov.uvfit[grpnamesindex, grpnamesindex],
errordf="approx",
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
type=type,
contrastmatrix=contrastmatrix,
n=with(uvfit$dfru,
sapply(split(endpt, grpf),
length)[grpnamesindex]),
offset=offset, display="none")
}
## slightly different from the original fit object settings slot
settings <- list(endptscale=settings$endptscale,
analysisname=settings$analysisname,
endptname=settings$endptname,
endptlabel=makeEndptLabel(settings$endptname,
settings$endptunits),
alpha=alpha,
mcadjust=mcadjust,
digits=settings$digits,
stamps=settings$stamps,
type=switch(type,
pairwise="Pairwise Comparisons (halfset)",
pairwisereflect="All Pairwise Comparisons",
allgroupstocontrol="All Groups versus Control",
custom=""),
sandaft=if(aft) TRUE else FALSE,
addpct=addpct)
returnObj <- new("cgOneFactorComparisonsTable",
ols.comprs=ols.comprs,
rr.comprs=rr.comprs,
aft.comprs=aft.comprs,
uv.comprs=uv.comprs,
settings=settings)
if(display=="print") {
print(returnObj)
}
else if(display=="show") {
show(returnObj)
}
## else display=="none"
invisible(returnObj)
}
)
pairwisecompsmatrix <- function(comps, grpnames) {
##
## PURPOSE: Function for rearranging pairwise comparisons results
## and placing in a matrix. The output matrix is designed to be
## used as input to formatting or further processing
##
## comps is assumed to be a data frame such as
## created by the comparisons function
grpnames <- validGrpNames(grpnames)
ngrps <- length(grpnames)
if(nrow(comps)!=ngrps*(ngrps-1)) {
stop(cgMessage("The number of rows in the comps dataframe",
"needs to have entries for all possible",
"pairwise comparisons, i.e. the number of groups",
"times the number of groups - 1"))
}
##ngrps <- round(uniroot(function(x) x^2 - x - nrow(comps),
## c(2, 1000))$root, 0)
indx <- matrix(1:(ngrps^2) - as.vector(col(diag(ngrps))), ncol=ngrps)
mcp.indx <- indx[lower.tri(indx)]
mcp.opp.indx <- setdiff(1:nrow(comps), mcp.indx)
mcp <- comps[mcp.indx,]
mcp.opp <- comps[mcp.opp.indx,]
## The combined rows of mcp and mcp.opp
## pertain to all possible unique pairwise comparisons
## amongst the groups. So we create a matrix for each statistical
## quantity that is of dimension ngrps by ngrps, putting the two
## triangular parts that come from mcp and mcp.opp together
## Construct list of each individual quantity
mcp.matrices <- lapply(seq(along=mcp),
function(i, x, y, ngrps) {
z <- matrix(NA, nrow=ngrps, ncol=ngrps)
z[lower.tri(z)] <- x[[i]]
z[upper.tri(z)] <- y[[i]]
return(z)
}, x=mcp, y=mcp.opp, ngrps=ngrps)
names(mcp.matrices) <- c("diff","se","lcl","ucl","pval")
## Now we have to combine the matrices into one big table
## which is probably intended for import into something like
## EXCEL for formatting eventually. This will involve
## interleaving rows and columns and adding some buffer sections to
## columns. A tricky part comes with the
## "lower CI and upper CI" part, because we want these to now be on the
## same row for a given comparison. The dimension of the final table
## matrix will be 4*ngrps rows by 2*ngrps columns.
mcpout.nrows <- 4*ngrps
mcpout.ncols <- 2*ngrps
mcp.out <- matrix(NA, nrow=mcpout.nrows, ncol=mcpout.ncols)
## fill-in the ucl's first
mcp.out[rep( c(rep(FALSE, mcpout.nrows),
rep(c(rep(FALSE, 3), TRUE), ngrps)), ngrps ) ] <- mcp.matrices$ucl
## second, the other quantities
mcp.out[rep( c(rep(TRUE, mcpout.nrows),
rep(FALSE, mcpout.nrows) ),
ngrps) ] <- do.call("rbind",
lapply(mcp.matrices[c("diff","se",
"pval",
"lcl")],
as.vector))
## lastly, matrix labels
therownames <- character(length(grpnames)*4)
thecolnames <- character(length(grpnames)*2)
therownames[seq(1, length(therownames), by=4)] <- grpnames
thecolnames[seq(1, length(thecolnames), by=2)] <- grpnames
grpnameheaders <- list(therownames, thecolnames)
dimnames(mcp.out) <- grpnameheaders
# return
mcp.out
}
setMethod("print", "cgOneFactorComparisonsTable",
print.cgOneFactorComparisonsTable <-
function(x, digits=NULL, title=NULL, endptname=NULL, ...) {
##
## PURPOSE: Semi-formatted print version of Comparisons Table
##
## NOTE: Had to use x as an argument because of the system defined
## generic. I would have preferred to use object; hence the first
## statement below.
##
object <- x
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
olscomprs <- object@ols.comprs
rrcomprs <- object@rr.comprs
aftcomprs <- object@aft.comprs
uvcomprs <- object@uv.comprs
ols <- rr <- uv <- aft <- FALSE ## Initializations
if(!is.null(aftcomprs)) {
aft <- TRUE
validArgModel(...)
}
else if(!is.null(uvcomprs)) {
uv <- TRUE
validArgModel(...)
}
if(!is.null(rrcomprs) && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(!is.null(olscomprs) && (model!="rronly") && !aft && !uv) {
ols <- TRUE
}
settings <- object@settings
alpha <- settings$alpha
mcadjust <- settings$mcadjust
addpct <- settings$addpct
if(is.null(digits)) {
digits <- settings$digits
}
else {
digits <- validArgDigits(digits)
}
curscipen <- getOption("scipen")
on.exit(options(scipen=curscipen), add=TRUE)
options(scipen=9)
if(is.null(title)) {
title <- paste("Comparisons Table of", settings$analysisname)
}
else {
validCharacter(title)
}
if(is.null(endptname)) {
endptname <- settings$endptname
if(!is.character(endptname)) {
endptname <- ""
}
}
else {
validCharacter(endptname)
}
cat(title,"\n")
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
diffmetric <- "Differences (A vs. B)"
if(settings$endptscale=="log") {
diffmetric <- paste("Percent", diffmetric)
}
cat(diffmetric, "\n")
fmtdig <- function(x, diffmetric, digits) {
pvalfmt <- fmtPvalue(x$pval)
fmt.x <- x
if(regexpr("Percent", diffmetric) > 0) {
fmt.x$estimate <- fmtPercent(x$estimate)
fmt.x$se <- fmtPercent(x$se)
fmt.x$lowerci <- fmtPercent(x$lowerci)
fmt.x$upperci <- fmtPercent(x$upperci)
fmt.x$pval <- pvalfmt
fmt.x$meanA <- fround(x$meanA, digits=digits)
fmt.x$seA <- fround(x$seA, digits=digits)
fmt.x$meanB <- fround(x$meanB, digits=digits)
fmt.x$seB <- fround(x$seB, digits=digits)
{ names(fmt.x)[is.element(names(fmt.x), c("meanA","seA",
"meanB","seB"))] <-
c("geomeanA","seA",
"geomeanB","seB") }
}
else { ## Simple Differences
fmt.x <- fround(x, digits)
if(addpct==TRUE) {
pctdiff.fmt <- fmtPercent(x$pctdiff)
fmt.x$pctdiff <- pctdiff.fmt
}
fmt.x$pval <- pvalfmt
}
return(fmt.x)
}
informConfidence <- function() {
cat(paste(round(100*(1-alpha), 0), "% Confidence ",
"(alpha of ", alpha, ")",
if(mcadjust) {", Multiplicity Adjusted"},
"\n",
sep=""))
}
curwidth <- getOption("width")
on.exit(options(width=curwidth), add=TRUE)
if(curwidth < 500) { options(width=500) }
if(ols) {
cat("\nClassical Least Squares Model Fit\n")
informConfidence()
print(fmtdig(olscomprs, diffmetric, digits), quote=FALSE)
}
if(rr) {
cat("\nResistant & Robust Model Fit\n")
informConfidence()
print(fmtdig(rrcomprs, diffmetric, digits), quote=FALSE)
}
if(aft) {
afttitle <- "\nAccelerated Failure Time Model Fit\n"
if(settings$sandaft) {
cat(paste(afttitle,
"with Sandwich Variance-Covariance Estimate\n",
sep=""))
}
else {
cat(afttitle)
}
informConfidence()
print(fmtdig(aftcomprs, diffmetric, digits), quote=FALSE)
}
if(uv) {
cat("\nUnequal Variances Model Fit\n")
informConfidence()
print(fmtdig(uvcomprs, diffmetric, digits), quote=FALSE)
}
invisible()
})
setMethod("show", "cgOneFactorComparisonsTable",
show.cgOneFactorComparisonsTable <- function(object) showDefault(object))
## Error Bar Graph
errorbargraph <- function(estimates, centralvar,
critpoint,
endptscale="log",
analysisname="",
endptname="",
alpha=0.05,
digits=NULL,
approxstamp=FALSE,
titlestamp=TRUE,
offset=NULL,
ticklabels=NULL, ...) {
##
## PURPOSE: Construct an error bar graph based on pairwise
## multiple comparisons (method of Andrews, Sarner, and Snee, 1980)
##
## argument handling
estimates <- validEstimates(estimates)
endptscale <- validArgMatch(endptscale, c("log","original"))
validAlpha(alpha)
validBoolean(approxstamp)
validBoolean(titlestamp)
options(warn=-1)
curpar <- par(new=FALSE, mgp=c(3,0.25,0), tck=-0.010)
options(warn=0)
on.exit(par(curpar))
numberofgrps <- length(estimates)
grpnames <- names(estimates)
errorbarlength <- critpoint*sqrt(2*centralvar)/2
bardata <- list(lower=estimates - errorbarlength,
upper=estimates + errorbarlength)
logscale <- ifelse(endptscale=="log", TRUE, FALSE)
rangedata <- if(logscale) exp(unlist(bardata)) else unlist(bardata)
if(logscale) {
parmar <- par(mar=c(5, 4, 4, 3) + 0.1)
curpar$mar <- parmar$mar
errbar(1:numberofgrps, estimates/log(10),
yplus=bardata$upper/log(10),
yminus=bardata$lower/log(10),
ylab=endptname, xlab="",
pch=16,
xlim=c(0.5, numberofgrps + 0.5), axes=FALSE)
mtext("log-spaced", side=2, line=2.25, cex=0.7)
tickmarks <- setupAxisTicks(rangedata,
logscale=logscale, digits=digits,
offset=offset)
if(!is.null(ticklabels)) {
tickmarks <- makeTickMarks(ticklabels, tickmarks,
offset=offset)
}
axis(2, at=log10(tickmarks), labels=names(tickmarks),
cex.axis=0.8, adj=1, las=1)
log10endpt <- unlist(bardata)/log(10)
axis(4, at=pretty(log10endpt), pretty(log10endpt), cex.axis=0.7,
adj=0, las=1, tck=-0.0075)
mtext(side=4, text=if(is.expression(endptname)) {
catCharExpr("Log10 scale of", endptname)
} else { paste("Log10 scale of", endptname, sep=" ") },
line=2, cex=0.7, adj=0)
}
else {
errbar(1:numberofgrps, estimates,
yplus=bardata$upper,
yminus=bardata$lower,
ylab=endptname, xlab="",
pch=16,
xlim=c(0.5, numberofgrps + 0.5), axes=FALSE)
tickmarks <- setupAxisTicks(rangedata,
logscale=logscale, digits=digits)
if(!is.null(ticklabels)) {
tickmarks <- makeTickMarks(ticklabels, tickmarks,
offset=offset)
}
axis(2, at=tickmarks, labels=names(tickmarks),
cex.axis=0.8, adj=1, las=1)
}
## Axes Customization
grpnameticksettings <- setupGrpNameTicks(grpnames, 1:numberofgrps)
plotGrpNameTicks((grpnames), settings=grpnameticksettings)
minmaxTicks(if(logscale) exp(unlist(bardata)) else
unlist(bardata),
theaxis="y", logscale=logscale, digits=digits,
offset=offset)
## Annotations
if(titlestamp) {
title(main=paste("Error Bar Graph\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(alphapercent=100*alpha)
}
if(approxstamp) errorBarGraphApproximateStamp()
box()
invisible()
}
setMethod("errorBarGraph", "cgOneFactorFit",
errorBarGraph.cgOneFactorFit <-
function(fit, mcadjust=FALSE,
alpha=0.05, cgtheme=TRUE, device="single", ...) {
##
## PURPOSE: Construct an error bar graph based on pairwise
## multiple comparisons (method of Andrews, Sarner, and Snee, 1980)
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names=c("model", "ticklabels"))
if(class(fit@uvfit)[1]=="gls" || class(fit@aftfit)[1]=="survreg") {
stop(cgMessage("There is no errorBarGraph method",
"defined for a fitted model that allowed",
"unequal variances or censored observations.",
"You may wish to look",
"at the errorbargraph function",
"to construct a graph if your standard errors",
"are not too different."))
}
validBoolean(cgtheme)
validAlpha(alpha)
validBoolean(mcadjust)
ticklabelsarg <- getDotsArgName(dots, "ticklabels")
if(!is.na(ticklabelsarg)) {
ticklabels <- eval(parse(text=paste("dots$", ticklabelsarg, sep="")))
validList(ticklabels, names=c("mod","marks"),
argname="ticklabels")
}
else {
ticklabels <- NULL
}
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly", "rronly"))
}
else {
model <- "both"
}
device <- validArgMatch(device, c("single", "multiple", "ask"))
settings <- fit@settings
offset <- settings$addconstant
rrfit <- fit@rrfit
olsfit <- fit@olsfit
digits <- settings$digits
analysisname <- settings$analysisname
endptlabel <- makeEndptLabel(settings$endptname, settings$endptunits)
endptscale <- settings$endptscale
stamps <- settings$stamps
alphapercent <- round(100*alpha, 0)
confquantile <- 1 - alpha/2
grpnames <- settings$grpnames
numberofgrps <- length(grpnames)
if(numberofgrps==2) { mcadjust <- FALSE }
L <- contrastMatrix(list(grpnames), "pairwise")
ols <- rr <- FALSE ## Initializations
if(class(rrfit)[1]=="rlm" && model!="olsonly") {
rr <- TRUE
}
if(class(olsfit)[1]=="lm" && model!="rronly") {
ols <- TRUE
}
df.residual <- olsfit$df.residual
thetitle <- "Error Bar Graph"
dfru <- olsfit$dfru
if(ols) {
n <- with(dfru, tapply(endpt, grpf, length))
## Calculate the harmonic mean to account for unequal sample sizes
nharmonic <- n ## default all n equal
if(length(unique(n)) > 1) {
nharmonic <- with(olsfit$dfru,
length(n) / sum(tapply(endpt, grpf,
function(x) {
1/length(x)} )
))
}
olsestimates <- olsfit$coef
olscentralstderr <- summary(olsfit)$sigma/sqrt(nharmonic)
olscentralvar <- olscentralstderr^2
olscritpoint <-
if(!mcadjust) { qt(confquantile, df.residual) }
else {
set.seed(17)
multcompInform()
glhtobj <- list(object=NULL, linfct=L, coef=olsestimates,
vcov=vcov(olsfit),
type="user-defined", alternative="two.sided",
df=df.residual,
rhs=rep(0, nrow(L)))
class(glhtobj) <- "glht"
## assign this value:
attr(confint(glhtobj, level=1-alpha)$confint, "calpha")
}
if(mcadjust) { multcompDone("Classical Least Squares") }
}
if(rr) {
rrestimates <- rrfit$coef
summ.rrfit <- summary(rrfit, method="XtWX", correlation=FALSE)
varcov.rrfit <- (summ.rrfit$stddev^2) * summ.rrfit$cov.unscaled
## a robust / resistant estimate of center of the std errors
rrcentralvar <- rlm(diag(varcov.rrfit) ~ 1)[[1]]
rrcentralstderr <- sqrt(rrcentralvar)
rrcritpoint <-
if(!mcadjust) { qt(confquantile, df.residual) }
else {
set.seed(17)
multcompInform()
glhtobj <- list(object=NULL, linfct=L, coef=olsestimates,
vcov=varcov.rrfit,
type="user-defined", alternative="two.sided",
df=df.residual,
rhs=rep(0, nrow(L)))
class(glhtobj) <- "glht"
## assign this value:
attr(confint(glhtobj, level=1-alpha)$confint, "calpha")
}
if(mcadjust) { multcompDone("Resistant & Robust") }
}
if(rr && ols && is.element(model, "both") && device=="single") {
ols.dfr <- data.frame(type=rep("Classical", numberofgrps),
grpf=factorInSeq(grpnames), estimate=olsestimates,
lower=olsestimates -
olscritpoint*olscentralstderr/sqrt(2),
upper=olsestimates +
olscritpoint*olscentralstderr/sqrt(2))
rr.dfr <- data.frame(type=rep("Resistant & Robust", numberofgrps),
grpf=factorInSeq(grpnames), estimate=rrestimates,
lower=rrestimates -
rrcritpoint*rrcentralstderr/sqrt(2),
upper=rrestimates +
rrcritpoint*rrcentralstderr/sqrt(2))
all.dfr <- rbind(ols.dfr, rr.dfr)
thetitle <- paste(thetitle, "s", sep="")
theapproxmsg <-
if(n[1]!=nharmonic[1]) {
"Error bars are approximate."
}
else {
"Error bars in Resistant & Robust panel are approximate."
}
cgDevice(cgtheme=cgtheme)
trellispanelstg <- trellis.par.get("clip")$panel
trellis.par.set("clip", list(panel="off"))
on.exit(trellis.par.set("clip", list(panel=trellispanelstg)),
add=TRUE)
trellisparstg2 <- trellis.par.get("layout.widths")$ylab
trellis.par.set("layout.widths", list(ylab=3))
on.exit(trellis.par.set("layout.widths", list(panel=trellisparstg2)),
add=TRUE)
trellisparstg3 <- trellis.par.get("axis.components")
trellis.par.set("axis.components",
list(left=list(pad1=0.2, pad2=2),
bottom=list(pad1=0.2, pad2=1),
top=list(pad1=1, pad2=1),
right=list(pad1=0.5, pad2=1)
))
on.exit(trellis.par.set("axis.components", trellisparstg3),
add=TRUE)
ally <- unlist(all.dfr[, c("estimate","lower","upper")])
if(endptscale=="log") {
thegraph <- xyplot(Cbind(log10(exp(estimate)),
log10(exp(lower)),
log10(exp(upper))) ~ as.numeric(grpf)
| type,
data=all.dfr,
digits=digits,
panel = function(x, y, ...) {
grpnameticksettings <- setupGrpNameTicks(grpnames,
grplocation=
1:length(grpnames),
cexinit=0.8,
cexthreshold=0.5,
grid=TRUE)
plotGrpNameTicks(grpnames, grpnameticksettings, grid=TRUE)
if(is.element(panel.number(),
c(1, nlevels(all.dfr$type)))) {
panelside <- ifelse(panel.number()==1,
"left", "right")
ycex <- 0.6
tickmarks <- setupAxisTicks(exp(ally),
logscale=TRUE,
grid=TRUE,
digits=digits,
offset=offset,
ycex=ycex)
if(!is.null(ticklabels)) {
tickmarks <- makeTickMarks(ticklabels, tickmarks,
offset=offset)
}
panel.axis(side=panelside,
at=log10(tickmarks),
labels=names(tickmarks),
tck=0.20, text.cex=ycex,
rot=0,
outside=TRUE)
}
panel.xYplot(x, y, subscripts=FALSE,
method="bars",
label.curves=FALSE, ...)
yrange <- range(unwind(attr(y, "other")))
panel.text(x=rep(0, 2), y=yrange,
labels=paste("",
c(fround(min(10^yrange), digits),
fround(max(10^yrange), digits))),
col="blue", adj=0, cex=0.6)
},
layout=c(2,1), aspect=1, pch=16,
col="black",
as.table=TRUE,
xlim=c(0, numberofgrps + 0.5),
ylim=rangeExtend(log10(exp(ally))),
ylab=list(cex=0.7,
label=endptlabel),
xlab="",
scales=list(x=list(at=1:numberofgrps,
labels=NULL, tck=c(0.15, 0), axs="i"),
y=list(labels=NULL, tck=0)),
main=list(label=paste(thetitle, "\n",
settings$analysisname, sep=""), cex=1.1),
par.strip.text=list(cex=0.7)
)
}
else {
thegraph <- xyplot(Cbind(estimate,
lower,
upper) ~ as.numeric(grpf) | type,
data=all.dfr,
digits=digits,
panel = function(x, y, ...) {
grpnameticksettings <- setupGrpNameTicks(grpnames,
grplocation=
1:length(grpnames),
cexinit=0.8,
cexthreshold=0.5,
grid=TRUE)
plotGrpNameTicks(grpnames, grpnameticksettings, grid=TRUE)
if(is.element(panel.number(),
c(1, nlevels(all.dfr$type)))) {
panelside <- ifelse(panel.number()==1,
"left", "right")
ycex <- 0.6
tickmarks <- setupAxisTicks(ally,
logscale=FALSE,
grid=TRUE,
digits=digits,
offset=offset,
ycex=ycex)
if(!is.null(ticklabels)) {
tickmarks <- makeTickMarks(ticklabels, tickmarks,
offset=offset)
}
panel.axis(side=panelside,
at=tickmarks,
labels=names(tickmarks),
tck=0.20, text.cex=ycex,
rot=0,
outside=TRUE)
}
panel.xYplot(x, y, subscripts=FALSE,
method="bars",
label.curves=FALSE, ...)
yrange <- range(unwind(attr(y, "other")))
panel.text(x=rep(0, 2), y=yrange,
labels=paste("",
c(fround(min(yrange), digits),
fround(max(yrange), digits))),
col="blue", adj=0, cex=0.6)
},
layout=c(2,1), aspect=1, pch=16,
col="black",
as.table=TRUE,
xlim=c(0, numberofgrps + 0.5),
ylim=rangeExtend(ally),
ylab=list(cex=0.7,
label=endptlabel),
xlab="",
scales=list(x=list(at=1:numberofgrps,
labels=NULL, tck=c(0.15, 0), axs="i"),
y=list(labels=NULL, tck=0)),
main=list(label=paste(thetitle, "\n",
settings$analysisname, sep=""), cex=1.1),
par.strip.text=list(cex=0.7)
)
}
print(thegraph) ## , position=c(0,0,0.95,1))
if(stamps) graphStampCG()
errorBarGraphStamp(mcadjust, alphapercent, grid=TRUE)
errorBarGraphApproximateStamp(grid=TRUE, msg=theapproxmsg)
}
else if((model=="olsonly" || (ols && !rr && model=="both")) &&
device=="single") {
errorbargraph(olsestimates, olscentralvar,
olscritpoint,
endptscale,
analysisname,
endptname=endptlabel,
alpha=alpha,
digits=digits,
approxstamp={nharmonic[1]!=n[1]},
titlestamp=FALSE, offset=offset,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Error Bar Graph, Classical analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(mcadjust, alphapercent)
}
else if(model=="rronly" && rr && device=="single") {
errorbargraph(rrestimates, rrcentralvar,
rrcritpoint,
endptscale,
analysisname,
endptname=endptlabel,
alpha=alpha,
digits=digits,
approxstamp=TRUE,
titlestamp=FALSE,
offset=offset,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Error Bar Graph, Resistant & Robust analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(mcadjust, alphapercent)
}
else if(rr && ols &&
is.element(device, c("ask","multiple")) &&
is.element(model, "both")) {
device <- validArgMatch(device, c("multiple", "ask"))
if(device=="ask") {
op <- par(ask = TRUE)
on.exit(par(op), add=TRUE)
}
errorbargraph(olsestimates, olscentralvar,
olscritpoint,
endptscale,
analysisname,
endptname=endptlabel,
alpha=alpha,
digits=digits,
approxstamp={nharmonic[1]!=n[1]},
titlestamp=FALSE, offset=offset,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Error Bar Graph, Classical analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(mcadjust, alphapercent)
if(device=="multiple") {
if(!is.null(dots$model)) dots$model <- NULL ## since we only
## want dots arguments for trellis.device in next call
do.call("cgDevice", c(list(new=TRUE), dots))
cat(cgMessage("A new graphics device has been generated",
"to hold the Resistant & Robust",
"errorBarGraph version.",
"The Classical Least Squares version is on the previous",
"device.\n",
warning=TRUE))
}
errorbargraph(rrestimates, rrcentralvar,
rrcritpoint,
endptscale,
analysisname,
endptname=endptlabel,
alpha=alpha,
digits=digits,
approxstamp=TRUE,
titlestamp=FALSE,
offset=offset,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Error Bar Graph, Resistant & Robust analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
errorBarGraphStamp(mcadjust, alphapercent)
}
else {
stop(cgMessage("The chosen device and model arguments",
"are not compatible either with each other",
"or with the fitted model(s) in the fit object.",
seeHelpFile("errorBarGraph")))
}
invisible()
}
)
grpsummary <- function(estimates,
varcovmatrix,
n, errordf=Inf,
endptscale="log", mcadjust=FALSE,
alpha=0.05,
offset=NULL,
analysisname="", endptname="",
digits=NULL, display="print", ...) {
##
## PURPOSE: Function for computing group summary results.
## The output matrix is designed to be
## used as input to formatting or further processing
validBoolean(mcadjust)
n <- validN(n, estimates)
estimates <- validEstimates(estimates)
endptscale <- validArgMatch(endptscale, c("log","original"))
validAlpha(alpha)
display <- validArgMatch(display, c("print","none","show"))
if(!is.numeric(errordf) || any(errordf < 1)) {
stop(cgMessage("The errordf argument needs to be numeric",
"and (each value) set to be 1 or greater."))
}
ngrps <- length(estimates)
grpnames <- names(estimates)
## If only two groups were analyzed, then multiplicity adjustment has
## no effect. So we ensure mcadjust is set to FALSE
if(ngrps==2) { mcadjust <- FALSE }
if(mcadjust) {
if(length(errordf) > 1) {
stop(cgMessage("The errordf can only be a single value",
"when mcadjust=FALSE.",
seeHelpFile("grpsummary")))
}
set.seed(17)
multcompInform()
glhtobj <- list(object=NULL, linfct=diag(ngrps),
coef=estimates,
vcov=varcovmatrix,
type="user-defined", alternative="two.sided",
df=errordf,
rhs=rep(0, nrow(diag(ngrps))))
class(glhtobj) <- "glht"
mcpci.fit <- confint(glhtobj, level=1-alpha)
mcptest.fit <- summary(glhtobj)
se <- as.vector(mcptest.fit)$test$sigma
lowerci <- mcpci.fit$confint[,"lwr"]
upperci <- mcpci.fit$confint[,"upr"]
}
else if(!mcadjust) {
se <- sqrt(diag(varcovmatrix))
tcrit <- qt(1 - alpha/2, errordf)
lowerci <- estimates - tcrit * se
upperci <- estimates + tcrit * se
}
thetable <- data.frame(n=n,
estimate=estimates, se=se,
lowerci=lowerci, upperci=upperci)
row.names(thetable) <- grpnames
## Handle log case if needed
## Transformations to Percent Difference Expression Scale
if(endptscale=="log") {
logscale.est <- thetable[, 2]
logscale.sd <- thetable[, 3]
thetable[, 2] <- exp(logscale.est) # Point Estimate
thetable[, 3] <- exp(logscale.est) * logscale.sd # Std Err
if(any(logscale.sd > 0.50)) {
warning(cgMessage("There is at least one",
"group standard error in the log scale",
"that exceeds 0.50, so the",
"estimated geometric mean standard",
"errors may be nonsensical and should be",
"cautiously regarded.")
)
}
thetable[, 4:5] <- exp(thetable[, 4:5]) # Confidence Limits
}
if(!is.null(offset)) {
validNumeric(offset, positive=TRUE)
thetable[, c(2,4,5)] <- (thetable[, c(2,4,5)] - offset)
}
if(display=="print") {
fmt.table <- thetable
cat(paste("Group Summary Table",
if(analysisname!="") paste("for", analysisname), "\n"))
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
metric <- "Means"
if(endptscale=="log") {
metric <- paste("Geometric", metric)
}
else {
metric <- paste("Arithmetic", metric)
}
cat(metric, "\n")
if(!is.null(digits)) {
fmt.table$estimate <- fround(fmt.table$estimate, digits)
fmt.table$se <- fround(fmt.table$se, digits)
fmt.table$lowerci <- fround(fmt.table$lowerci, digits)
fmt.table$upperci <- fround(fmt.table$upperci, digits)
}
fmt.table$n <- fround(fmt.table$n, 0)
cat(paste(round(100*(1-alpha), 0), "% Confidence ",
"(alpha of ", alpha, ")",
if(mcadjust) {", Multiplicity Adjusted"},
"\n",
sep=""))
curwidth <- getOption("width")
on.exit(options(width=curwidth), add=TRUE)
if(curwidth < 500) { options(width=500) }
print(fmt.table, quote=FALSE)
return(invisible(thetable))
}
else if (display=="show") {
showDefault(thetable)
}
## else show nothing
invisible(thetable)
}
setClass("cgOneFactorGrpSummaryTable",
representation(ols.grps="dataframeOrNULL",
rr.grps="dataframeOrNULL",
aft.grps="dataframeOrNULL",
uv.grps="dataframeOrNULL",
settings="list"),
prototype(ols.grps=NULL, rr.grps=NULL, aft.grps=NULL, uv.grps=NULL,
settings=list()))
setMethod("grpSummaryTable", "cgOneFactorFit",
function(fit,
mcadjust=FALSE,
alpha=0.05, display="print", ...) {
##
## PURPOSE: Summarize each group distribution based
## on the model fit.
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
settings <- fit@settings
validAlpha(alpha)
validBoolean(mcadjust)
display <- validArgMatch(display, c("print","none","show"))
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly", "rronly"))
}
else {
model <- "both"
}
##
## initializations
aft <- ols <- rr <- uv <- FALSE
aft.grps <- ols.grps <- rr.grps <- uv.grps <- NULL
##
endptscale <- settings$endptscale
rrfit <- fit@rrfit
olsfit <- fit@olsfit
aftfit <- fit@aftfit
uvfit <- fit@uvfit
grpnames <- settings$grpnames
offset <- settings$addconstant
if(class(aftfit)[1]=="survreg") {
aft <- TRUE
validArgModel(...)
}
else if(class(uvfit)[1]=="gls") {
uv <- TRUE
validArgModel(...)
}
if(class(rrfit)[1]=="rlm" && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(class(olsfit)[1]=="lm" && model!="rronly" && !aft && !uv) {
ols <- TRUE
if(!rr) model <- "olsonly"
}
n <- with(olsfit$dfru, sapply(split(endpt, grpf), length))
if(rr || ols) {
df.residual <- olsfit$df.residual
}
if(ols) {
olsestimates <- olsfit$coef
varcov.olsfit <- vcov(olsfit)
ols.grps <- grpsummary(olsestimates,
varcov.olsfit,
n=n, errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
offset=offset, display="none", ...)
if(mcadjust) { multcompDone("Classical Least Squares") }
}
if(rr) {
rrestimates <- rrfit$coef
summ.rrfit <- summary(rrfit, method="XtWX")
stddev <- summ.rrfit$stddev
varcov.rrfit <- (stddev^2) * summ.rrfit$cov.unscaled
rr.grps <- grpsummary(rrestimates,
varcov.rrfit,
n=n, errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
offset=offset, display="none", ...)
if(mcadjust) { multcompDone("Resistant & Robust") }
}
if(aft) {
aftestimates <- aftfit$coef
df.residual <- aftfit$df.residual
varcov.aftfit <- chop.matrix(vcov(aftfit))
aft.grps <- grpsummary(aftestimates,
varcov.aftfit,
n=n, errordf=df.residual,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
offset=offset, display="none",...)
}
if(uv) {
if(mcadjust) {
stop(cgMessage("The mcadjust argument cannot be set",
"to TRUE when the fitted model allowed",
"for unequal variances."))
}
uv.grps <- grpsummary(uvfit$coef,
vcov(uvfit),
n=n, errordf=n-1,
endptscale=endptscale,
mcadjust=mcadjust,
alpha=alpha,
offset=offset, display="none",...)
}
settings <- list(endptscale=settings$endptscale,
analysisname=settings$analysisname,
endptname=settings$endptname,
alpha=alpha,
mcadjust=mcadjust,
digits=settings$digits,
sandaft=if(aft) TRUE else FALSE)
returnObj <- new("cgOneFactorGrpSummaryTable",
ols.grps=ols.grps, rr.grps=rr.grps,
aft.grps=aft.grps, uv.grps=uv.grps,
settings=settings)
if(display=="print") {
print(returnObj)
}
else if(display=="show") {
show(returnObj)
}
## else display=="none"
invisible(returnObj)
})
setMethod("print", "cgOneFactorGrpSummaryTable",
print.cgOneFactorGrpSummaryTable <-
function(x, digits=NULL, title=NULL, endptname=NULL, ...) {
##
## PURPOSE: Semi-formatted print version of Comparisons Table
##
## NOTE: Had to use x as an argument because of the system defined
## generic. I would have preferred to use object; hence the first
## statement below.
##
object <- x
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names="model")
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
model <- "both"
}
olsgrps <- object@ols.grps
rrgrps <- object@rr.grps
aftgrps <- object@aft.grps
uvgrps <- object@uv.grps
ols <- rr <- uv <- aft <- FALSE ## Initializations
if(!is.null(aftgrps)) {
aft <- TRUE
validArgModel(...)
}
else if(!is.null(uvgrps)) {
uv <- TRUE
validArgModel(...)
}
if(!is.null(rrgrps) && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(!is.null(olsgrps) && (model!="rronly") && !aft && !uv) {
ols <- TRUE
}
settings <- object@settings
alpha <- settings$alpha
mcadjust <- settings$mcadjust
if(is.null(digits)) {
digits <- settings$digits
}
else {
digits <- validArgDigits(digits)
}
curscipen <- getOption("scipen")
on.exit(options(scipen=curscipen), add=TRUE)
options(scipen=9)
if(is.null(title)) {
title <- paste("Group Summary Table of", settings$analysisname)
}
else {
validCharacter(title)
}
if(is.null(endptname)) {
endptname <- settings$endptname
if(!is.character(endptname)) {
endptname <- ""
}
}
else {
validCharacter(endptname)
}
cat(title,"\n")
if(endptname!="") { cat(paste("Endpoint:", endptname, "\n")) }
metric <- "Means"
if(settings$endptscale=="log") {
metric <- paste("Geometric", metric)
}
else {
metric <- paste("Arithmetic", metric)
}
cat(metric, "\n")
fmtdig <- function(x, digits) {
x$estimate <- fround(x$estimate, digits)
x$se <- fround(x$se, digits)
x$lowerci <- fround(x$lowerci, digits)
x$upperci <- fround(x$upperci, digits)
x$n <- fround(x$n, 0)
return(x)
}
informConfidence <- function() {
cat(paste(round(100*(1-alpha), 0), "% Confidence ",
"(alpha of ", alpha, ")",
if(mcadjust) {", Multiplicity Adjusted"},
"\n",
sep=""))
}
curwidth <- getOption("width")
on.exit(options(width=curwidth), add=TRUE)
if(curwidth < 500) { options(width=500) }
if(ols) {
cat("\nClassical Least Squares Model Fit\n")
informConfidence()
print(fmtdig(olsgrps, digits), quote=FALSE)
}
if(rr) {
cat("\nResistant & Robust Model Fit\n")
informConfidence()
print(fmtdig(rrgrps, digits), quote=FALSE)
}
if(aft) {
afttitle <- "\nAccelerated Failure Time Model Fit\n"
if(settings$sandaft) {
cat(paste(afttitle,
"with Sandwich Variance-Covariance Estimate\n",
sep=""))
}
else {
cat(afttitle)
}
informConfidence()
print(fmtdig(aftgrps, digits), quote=FALSE)
}
if(uv) {
cat("\nUnequal Variances Model Fit\n")
informConfidence()
print(fmtdig(uvgrps, digits), quote=FALSE)
}
invisible()
})
setMethod("show", "cgOneFactorGrpSummaryTable",
show.cgOneFactorGrpSummaryTable <- function(object) showDefault(object))
comparisonsgraph <- function(compstable,
difftype,
analysisname="",
endptname="",
alpha=0.05,
digits=NULL,
titlestamp=TRUE,
explanation=TRUE,
wraplength=20,
cex.comps=0.7,
ticklabels=NULL,
...) {
##
## PURPOSE: Plot Confidence Intervals of Comparisons;
## use horizontal error bars
##
## compstable needs to be a dataframe of proper format
difftype <- validArgMatch(difftype, c("percent","amount","simple"))
validAlpha(alpha)
validBoolean(titlestamp)
validBoolean(explanation)
digits <- validArgDigits(digits)
options(warn=-1)
curpar <- par(new=FALSE, mgp=c(3,0.25,0), tck=-0.010, mar=c(5,7,4,2)+0.1)
options(warn=0)
on.exit(par(curpar), add=TRUE)
numberofcomprs <- nrow(compstable)
comprnames <- row.names(compstable)
est <- compstable[, "estimate"]
lower <- compstable[, "lowerci"]
upper <- compstable[, "upperci"]
compys <- numberofcomprs:1
if(difftype=="percent") {
est <- log(pctToRatio(est))
lower <- log(pctToRatio(lower))
upper <- log(pctToRatio(upper))
}
allx <- c(lower, upper, est)
## We include zero because we always want the
## "no-difference" reference line to be visible
xmin <- min(c(0, allx))
xmax <- max(c(0, allx))
if(difftype=="percent") {
## snippet from Hmisc's errbar function
plot(est/log(10), compys, ylim = c(0.5, numberofcomprs + 0.5),
xlim = c(xmin, xmax)/log(10),
ylab="",
xlab = "",
pch = 16, axes=FALSE, ...)
mtext(side=1, text="log-spaced", line=4, cex=0.7)
ycoord <- par()$usr[3:4]
segments(lower/log(10), compys, upper/log(10), compys)
smidge <- 0.015 * (ycoord[2] - ycoord[1])/2
segments(lower/log(10), compys - smidge, lower/log(10), compys + smidge)
segments(upper/log(10), compys - smidge, upper/log(10), compys + smidge)
## zero reference line
abline(v=0, lty=3)
## Axes Customization
log10endpt <- log10(exp(allx))
xratioticks <- setupAxisTicks(exp(allx), ratio=TRUE, percent=TRUE,
axis="x", digits=digits)
if(!is.null(ticklabels)) {
xratioticks <- makeTickMarks(ticklabels, xratioticks,
percent=TRUE)
}
axis(1,at=log10(xratioticks),
labels=names(xratioticks),
tck=-.010, cex.axis=0.8)
xlabchar <- "Percent Difference"
}
else { ## difftype = "amount" or "simple" (Simple Difference)
plot(est, compys, ylim = c(0.5, numberofcomprs + 0.5),
xlim = c(xmin, xmax),
ylab="",
xlab="",## paste("Difference in\n ", endptname, sep=""),
pch=16, axes=FALSE, ...)
ycoord <- par()$usr[3:4]
segments(lower, compys, upper, compys)
smidge <- 0.015 * (ycoord[2] - ycoord[1])/2
segments(lower, compys - smidge, lower, compys + smidge)
segments(upper, compys - smidge, upper, compys + smidge)
## zero reference line
abline(v=0, lty=3)
## Axes Customization
xdiffticks <- setupAxisTicks((allx), ratio=FALSE, difference=TRUE,
logscale=FALSE,
axis="x", digits=digits)
if(!is.null(ticklabels)) {
xdiffticks <- makeTickMarks(ticklabels, xdiffticks,
percent=FALSE)
}
axis(1, at=xdiffticks,
labels=names(xdiffticks),
tck=-.010, cex=0.8)
xlabchar <- "Difference"
}
mtext(side=1, text=xlabchar, line=2)
if(is.expression(endptname)) {
mtext(side=1, text=catCharExpr("in ", endptname), line=3)
}
else {
mtext(side=1, text=paste("in", endptname), line=3)
}
## Wrap long comparison labels
for(i in seq(along=comprnames)) {
##if(nchar(comprnames[i]) > wraplength) {
## comprnames[i] <- gsub("vs.","vs.\n", comprnames[i])
##}
if(nchar(comprnames[i]) > wraplength) {
comprnames[i] <- paste(strwrap(comprnames[i], width=wraplength,
exdent=1), collapse="\n")
}
}
axis(2, at=compys, labels=comprnames,
adj=1, cex.axis=cex.comps, las=1)
box()
## Titles
par(curpar)
if(explanation) {
comparisonsGraphStamp(alphapercent=100*alpha)
}
if(titlestamp) {
title(main=paste("Comparisons Graph\n",
analysisname, sep=""), line=2, cex.main=1.1)
}
## min-max annotation
bardata <- c(est, lower, upper)
minmaxTicks(if(difftype=="percent") exp(bardata) else bardata,
theaxis="x", percent=if(difftype=="percent") TRUE else FALSE,
ratio=if(difftype=="percent") TRUE else FALSE,
logscale=if(difftype=="percent") TRUE else FALSE,
digits=digits)
## Add 0
text(y=par("usr")[3], x=0,
labels="0\n", col="blue", adj=0, cex=0.7)
invisible()
}
setMethod("comparisonsGraph", "cgOneFactorComparisonsTable",
comparisonsGraph.cgOneFactorComparisonsTable <-
function(compstable, cgtheme=TRUE, device="single",
wraplength=20, cex.comps=0.7, ...) {
##
## PURPOSE: create a graph of CI's on differences
## to judge comparisons
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names=c("model", "ticklabels"))
validBoolean(cgtheme)
settings <- compstable@settings
difftype <- if(settings$endptscale=="log") "percent" else "simple"
alpha <- settings$alpha
alphapercent <- round(100*alpha, 0)
mcadjust <- settings$mcadjust
analysisname <- settings$analysisname
endptlabel <- settings$endptlabel
stamps <- settings$stamps
digits <- settings$digits
rrcomprs <- compstable@rr.comprs
olscomprs <- compstable@ols.comprs
aftcomprs <- compstable@aft.comprs
uvcomprs <- compstable@uv.comprs
ticklabelsarg <- getDotsArgName(dots, "ticklabels")
if(!is.na(ticklabelsarg)) {
ticklabels <- eval(parse(text=paste("dots$", ticklabelsarg, sep="")))
validList(ticklabels, names=c("mod","marks"),
argname="ticklabels")
}
else {
ticklabels <- NULL
}
modelarg <- getDotsArgName(dots, "model")
if(!is.na(modelarg)) {
model <- eval(parse(text=paste("dots$", modelarg, sep="")))
model <- validArgMatch(model, choices=c("both", "olsonly","rronly"))
}
else {
dots$model <- "both"
model <- "both"
}
aft <- ols <- rr <- uv <- FALSE ## Initializations
if(!is.null(aftcomprs)) {
aft <- TRUE
validArgModel(...)
}
else if(!is.null(uvcomprs)) {
uv <- TRUE
validArgModel(...)
}
if(!is.null(rrcomprs) && model!="olsonly" && !aft && !uv) {
rr <- TRUE
}
if(!is.null(olscomprs) && (model!="rronly") && !aft && !uv) {
ols <- TRUE
}
device <- validArgMatch(device, c("single","multiple", "ask"))
thetitle <- "Comparisons Graph"
if(rr && ols && is.element(model, "both") && device=="single") {
all.dfr <- rbind(olscomprs, rrcomprs)
all.dfr$typef <- factorInSeq(c(rep("Classical",
nrow(olscomprs)),
rep("Resistant & Robust",
nrow(rrcomprs))))
thetitle <- paste(thetitle, "s", sep="")
numberofcomprs <- nrow(olscomprs)
comprnames <- row.names(olscomprs)
cgDevice(cgtheme=cgtheme)
trellispanelstg <- trellis.par.get("clip")$panel
trellis.par.set("clip", list(panel="off"))
on.exit(trellis.par.set("clip", list(panel=trellispanelstg)),
add=TRUE)
trellisparstg2 <- trellis.par.get("layout.widths")
ymargin <- max(nchar(comprnames))/2
trellis.par.set("layout.widths", list(left.padding=ymargin,
ylab=ymargin))
on.exit(trellis.par.set("layout.widths", trellisparstg2),
add=TRUE)
trellisparstg3 <- trellis.par.get("axis.components")
trellis.par.set("axis.components",
list(left=list(pad1=0.2, pad2=2),
bottom=list(pad1=0.2, pad2=1),
top=list(pad1=1, pad2=1),
right=list(pad1=0.5, pad2=1)
))
on.exit(trellis.par.set("axis.components", trellisparstg3),
add=TRUE)
allx <- unlist(all.dfr[, c("estimate","lowerci","upperci")])
est <- all.dfr[, "estimate"]
lower <- all.dfr[, "lowerci"]
upper <- all.dfr[, "upperci"]
compys <- numberofcomprs:1
if(difftype=="percent") {
est <- log(pctToRatio(est))
lower <- log(pctToRatio(lower))
upper <- log(pctToRatio(upper))
allx <- log(pctToRatio(allx))
xmin <- min(c(0, allx))
xmax <- max(c(0, allx))
## ideas in the next call adapted from errbar function
thegraph <- xyplot(compys ~ Cbind(est,
lower,
upper)
| typef,
data=all.dfr,
digits=digits,
panel = function(x, y, ...) {
xlower <- attr(x, "other")[, 1]
xupper <- attr(x, "other")[, 2]
xcex <- 0.6
xratioticks <-
setupAxisTicks(exp(c(xlower, x, xupper)),
ratio=TRUE, percent=TRUE,
axis="x",
grid=TRUE, xcex=xcex,
digits=digits)
if(!is.null(ticklabels)) {
xratioticks <-
makeTickMarks(ticklabels,
xratioticks,
percent=TRUE)
}
panel.axis(side="bottom",
at=log(xratioticks),
labels=names(xratioticks),
text.cex=xcex,
tck=0.20, outside=TRUE, rot=0)
panel.points(x, compys, pch=16, col="black")
lsegments(xlower, compys, xupper, compys)
##ycoord <- par()$usr[3:4]
##smidge <- 0.015 * (ycoord[2] - ycoord[1])
smidge <- 0.015 * max(compys)
lsegments(xlower, compys - smidge,
xlower, compys + smidge)
lsegments(xupper, compys - smidge,
xupper, compys + smidge)
## zero reference line
panel.abline(v=0, lty=3)
## Add 0 if not present
if(all(xratioticks!=1)) {
panel.axis(side="bottom",
at=0, labels="0", text.cex=0.6,
tck=0.20, outside=TRUE)
}
if(panel.number()==1) {
## Wrap long comparison labels
for(i in seq(along=comprnames)) {
if(nchar(comprnames[i]) > wraplength) {
comprnames[i] <-
paste(strwrap(comprnames[i],
width=wraplength,
exdent=1),
collapse="\n")
}
}
panel.axis(side="left",
at=compys,
labels=comprnames,
tck=0.20, text.cex=cex.comps,
rot=0,
outside=TRUE)
}
xrange <- range(c(xlower, xupper))
panel.text(y=rep(0.5, 2), x=xrange,
labels=paste(
c(fmtRatioToPercent(exp(min(xrange))),
fmtRatioToPercent(exp(max(xrange)))), "\n"),
col="blue", adj=0, cex=0.6)
},
layout=c(2,1), aspect=1, pch=16,
col="black",
as.table=TRUE,
ylim = c(0.5, numberofcomprs + 0.5),
xlim=rangeExtend(c(xmin, xmax), pct=10),
xlab="Percent Difference",
ylab="",
scales=list(y=list(at=1:numberofcomprs,
labels=NULL, tck=c(0.15, 0), axs="i"),
x=list(labels=NULL, tck=0)),
main=list(label=paste(thetitle, "\n",
analysisname, sep=""), cex=1.1),
par.strip.text=list(cex=0.7)
)
}
else {
xmin <- min(c(0, allx))
xmax <- max(c(0, allx))
## ideas in the next call adapted from Hmisc::errbar function
thegraph <- xyplot(compys ~ Cbind(est,
lower,
upper)
| typef,
data=all.dfr,
digits=digits,
panel = function(x, y, ...) {
xlower <- attr(x, "other")[, 1]
xupper <- attr(x, "other")[, 2]
xcex <- 0.6
xdiffticks <-
setupAxisTicks(c(xlower, x, xupper),
difference=TRUE,
axis="x",
logscale=FALSE,
grid=TRUE,
xcex=0.6,
digits=digits)
if(!is.null(ticklabels)) {
xdiffticks <-
makeTickMarks(ticklabels,
xdiffticks,
percent=FALSE)
}
panel.axis(side="bottom",
at=xdiffticks,
labels=names(xdiffticks),
text.cex=xcex,
tck=0.20, outside=TRUE, rot=0)
panel.points(x, compys, pch=16, col="black")
lsegments(xlower, compys, xupper, compys)
## ycoord <- par()$usr[3:4]
## smidge <- 0.015 * (ycoord[2] - ycoord[1])
smidge <- 0.015 * max(compys)
lsegments(xlower, compys - smidge,
xlower, compys + smidge)
lsegments(xupper, compys - smidge,
xupper, compys + smidge)
## zero reference line
panel.abline(v=0, lty=3)
## Add 0 if not present
if(all(xdiffticks!=0)) {
panel.axis(side="bottom",
at=0, labels="0", text.cex=0.6,
tck=0.20, outside=TRUE)
}
if(panel.number()==1) {
## Wrap long comparison labels
for(i in seq(along=comprnames)) {
if(nchar(comprnames[i]) > wraplength) {
comprnames[i] <-
paste(strwrap(comprnames[i],
width=wraplength,
exdent=1),
collapse="\n")
}
}
panel.axis(side="left",
at=compys,
labels=comprnames,
tck=0.20, text.cex=cex.comps,
rot=0,
outside=TRUE)
}
xrange <- range(c(xlower, xupper))
panel.text(y=rep(0.5, 2), x=xrange,
labels=paste(
c(fmtDifference(min(xrange),
digits=digits),
fmtDifference(max(xrange),
digits=digits)), "\n"),
col="blue", adj=0, cex=0.6)
},
layout=c(2,1), aspect=1, pch=16,
col="black",
as.table=TRUE,
ylim=c(0.5, numberofcomprs + 0.5),
xlim=rangeExtend(c(xmin, xmax), pct=10),
xlab="Difference",
ylab="",
scales=list(y=list(at=1:numberofcomprs,
labels=NULL, tck=c(0.15, 0), axs="i"),
x=list(labels=NULL, tck=0)),
main=list(label=paste(thetitle, "\n",
analysisname, sep=""), cex=1.1),
par.strip.text=list(cex=0.7)
)
}
print(thegraph)
seekViewport(trellis.vpname("xlab"))
grid.text(catCharExpr("in ", endptlabel), y = unit(-1, "lines"))
upViewport(0)
comparisonsGraphStamp(mcadjust, alphapercent, grid=TRUE,
desc=settings$type)
if(stamps) graphStampCG()
}
else if((model=="olsonly" || (ols && !rr && model=="both")) &&
device=="single") {
comparisonsgraph(olscomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Classical analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else if((model=="rronly" || (!ols && rr && model=="both")) &&
device=="single") {
## else if(model=="rronly" & rr & device=="single") {
comparisonsgraph(rrcomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Resistant & Robust analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else if(rr && ols &&
is.element(device, c("ask","multiple")) &&
is.element(model, "both")) {
device <- validArgMatch(device, c("multiple", "ask"))
if(device=="ask") {
op <- par(ask = TRUE)
on.exit(par(op), add=TRUE)
}
comparisonsgraph(olscomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Classical analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
if(device=="multiple") {
## dots$model <- NULL ## avoid unused argument error
eval(parse(text=paste("dots$", modelarg, " <- NULL", sep="")))
do.call("cgDevice", c(list(new=TRUE), dots))
cat(cgMessage("A new graphics device has been generated",
"to hold the Resistant & Robust",
"Comparisons Graph version.",
"The Classical Least Squares version is on the previous",
"device.\n",
warning=TRUE))
}
comparisonsgraph(rrcomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Resistant & Robust analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else if(aft && device=="single") {
comparisonsgraph(aftcomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
explanation=FALSE, ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, AFT analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else if(uv && device=="single") {
comparisonsgraph(uvcomprs, difftype, analysisname,
endptlabel, alpha, digits, titlestamp=FALSE,
explanation=FALSE,
wraplength=wraplength,
cex.comps=cex.comps,
ticklabels=ticklabels)
if(stamps) graphStampCG(grid=FALSE)
## Text Annotations
title(main=paste("Comparisons Graph, Unequal Variances analysis\n",
analysisname, sep=""), line=2, cex.main=1.1)
comparisonsGraphStamp(mcadjust, alphapercent,
desc=settings$type)
}
else {
stop(cgMessage("The chosen device and model arguments",
"are not compatible either with each other",
"or with the fitted model(s) in the fit object.",
seeHelpFile("comparisonsGraph")))
}
invisible()
}
)
validErrorDf <- function(errordf, varcovmatrix,
n, mcadjust) {
##
if(errordf=="approx") {
## First check that estimates are uncorrelated, i.e
## the varcovmatrix is diagonal.
if(sum(varcovmatrix[lower.tri(varcovmatrix)],
varcovmatrix[upper.tri(varcovmatrix)])!=0) {
stop(cgMessage("The varcovmatrix argument appears to have",
"nonzero off-diagonal elements. Currently,",
"the \"approx\" setting for the errordf",
"argument will only work for uncorrelated estimates."))
}
if(!isAllEqual(c(length(diag(varcovmatrix)),
length(n)))) {
stop(cgMessage("The lengths of the n argument and",
"the number of diagonal elements in the varcovmatrix",
"argument appear to differ. They need to be the same."))
}
if(mcadjust) {
stop(cgMessage("The mcadjust argument must be FALSE when the",
"errordf argument is \"approx\"."))
}
}
else if(length(errordf)!=1 || !is.numeric(errordf) || errordf < 1) {
stop(cgMessage("The errordf argument needs to be length 1",
"and set to be 1 or greater or the value \"approx\"."))
}
else return(TRUE)
}
validComparisonType <- function(type, errordf) {
type <- try(match.arg(type, c("pairwisereflect","pairwise",
"allgroupstocontrol", "custom")))
if(class(type)=="try-error") {
stop(cgMessage("The type argument must be one of \"pairwisereflect\",",
"\"pairwise\", \"allgroupstocontrol\",",
"or \"custom\"."))
}
if(errordf=="approx" && !any(is.element(type,
c("pairwisereflect","pairwise",
"allgroupstocontrol")))) {
stop(cgMessage("The type argument must be one of \"pairwisereflect\",",
"\"pairwise\", or \"allgroupstocontrol\"",
"when the errordf argument is \"approx\"."))
}
else if(!any(is.element(type, c("pairwisereflect","pairwise",
"allgroupstocontrol", "custom")))) {
stop(cgMessage("The type argument must be one of \"pairwisereflect\",",
"\"pairwise\", \"allgroupstocontrol\",",
"or \"custom\"."))
}
else return(TRUE)
}
validEstimates <- function(x) {
if(!is.numeric(x) || length(x) < 2 ) {
stop(cgMessage("The estimates argument needs to be numeric",
"and be of length two or more."
))
}
if(is.null(thenames <- names(x))) {
stop(cgMessage("The estimates arguments needs to have",
"a names attribute to use for",
"thes group labels."))
}
if(any(is.na(thenames)) || any(thenames=="")) {
stop(cgMessage("One or more of the groups names in the",
"names attribute of estimates",
"appears to be missing."))
}
else return(x)
}
validGrpNames <- function(x) {
if(!is.character(x) || length(x) < 2 ) {
stop(cgMessage("The grpnames argument needs to be character",
"and be of length two or more.",
))
}
if(any(is.na(x)) || any(x=="")) {
stop(cgMessage("One or more of the groups names in the",
"names attribute of estimates",
"appears to be missing.",
))
}
else return(x)
}
validN <- function(n, estimates) {
if(!is.numeric(n) || any(n < 1)) {
stop(cgMessage("The n argument vector needs to be numeric",
"and all elements must be one or greater."))
}
else if(length(n)!=length(estimates)) {
stop(cgMessage("The n argument vector needs to be of the",
"same length as the estimates vector."))
}
else return(as.integer(round(n, 0)))
}
validCnames <- function(type, x, reqlength) {
if(length(x)==1 && x=="derive") return(TRUE)
if(type!="custom") {
stop(cgMessage("The cnames argument must be left at its",
"default value of \"derive\" when type is",
"not \"custom\"."))
}
if(length(x) != reqlength) {
stop(cgMessage("The cnames argument must be either the",
"value \"derive\", ",
"or must be of length", reqlength,
", the number of rows of the",
"contrastmatrix argument."))
}
return(TRUE)
}
validAddPct <- function(addpct, endptscale) {
if(addpct && endptscale=="log") {
stop(cgMessage("The addpct argument cannot be TRUE,",
"and must be set to FALSE,",
"when the \"endptscale\" argument is",
"\"log\"."))
}
return(TRUE)
}
Try the cg package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.