R/p05SampleSizeOneFactorData.R In cg: Compare Groups, Analytically and Graphically

Documented in samplesizesamplesizegraphvalidDenDf

```## \$Id: p05SampleSizeOneFactorData.R 6053 2015-02-22 20:23:45Z bpikouni \$
## One-Factor Unpaired Groups Case

## Sample Size Calculations for One-Factor Unpaired Groups Data

samplesize <- function(sigmaest, endptscale,
planningname="", endptname="",
ngrps=2, direction="increasing",
mmdvec, dendf, ncp=NULL,
alpha=0.05, power=0.80,
nmax=1000, display="print", ...) {
##
## PURPOSE: Compute sample sizes based on variability and the other usually
## needed specifications. Idea of minimum difference with global
## F test is implemented.
##
## Input argument checking
validNumeric(ngrps, positive=TRUE, integer=TRUE)
if(ngrps < 2) {
stop(cgMessage("The ngrps argument needs to be an integer of 2 or greater."))
}
validNumeric(sigmaest, positive=TRUE)
validCharacter(planningname)
endptscale <- validArgMatch(endptscale, c("log", "original"))
direction <- validArgMatch(direction, c("increasing", "decreasing"))
validNumeric(mmdvec, positive=TRUE)
validAlpha(alpha)
validPower(power)
validDenDf(dendf)
validNcp(ncp)
validNumeric(nmax, positive=TRUE, integer=TRUE)
if(nmax < 2) {
stop(cgMessage("The nmax argument needs to be an integer of 2 or greater."))
}

if(endptscale=="log") {
if(direction=="decreasing" && any(mmdvec >= 100)) {
stop(cgMessage("Percent Change Values of 100% or greater reduction",
"seem to have been entered.  Since this is",
"numerically impossible, please check and delete any such",
"values."))
}
}

display <- validArgMatch(display, c("print","none","show"))
dots <- list(...)
validDotsArgs(dots, names=c(""))

nsolution <- vector("numeric", length=length(mmdvec))
numdf <- ngrps - 1

## main loop to calculate sample size
for(i in seq(along=mmdvec)) {
n <- 2 # initialization
mmdveci <- mmdvec[i]

if(endptscale=="log") {
if(direction=="decreasing") {
sgn <- -1
}
else { sgn <- 1 }
mmdveci <- log(sgn*mmdveci/100 + 1)
}

repeat {
dendf.result <- do.call("dendf", list(n, ngrps))

if(is.null(ncp)) {
## canonical case of one-factor unpaired groups
ncp.result <- ( n * (mmdveci)^2 ) / ( 2 * sigmaest^2 )
}
else {
## take the argument for specialized ncp function
## which requires this specific order of argument
## specification: sigmaest, mmdveci, and n
ncp.result <- do.call("ncp", list(sigmaest, mmdveci, n))
}
pwr <- 1 - pf(qf(1 - alpha, numdf, dendf.result),
numdf, dendf.result, ncp.result)
if (pwr > power || n >= nmax) {
nfinal <- n
break
}
else {
n <- n + 1
}
}
nsolution[i] <- nfinal
}

if(any(nsolution==nmax)) {
cat(cgMessage(" With the variability set at ",
signif(sigmaest, 4),
", the nmax threshold specified at ",
nmax,
"was reached for at least ",
"one of the specified differences",
"\n\n",
warning=TRUE))
}

computedss <- cbind(mmdvec, n=nsolution, N=ngrps*nsolution)

dimnames(computedss)[[2]] <- c("mmd","n","N")
n <- computedss[,"n"]
N <- computedss[,"N"]

if(display=="print") {
fmtcomputedss <- as.data.frame(computedss)
fmtcomputedss\$mmd <- fround(fmtcomputedss\$mmd, getNumDigits(fmtcomputedss\$mmd))
rownames(fmtcomputedss) <- fmtcomputedss\$mmd
fmtcomputedss\$mmd <- NULL
names(fmtcomputedss) <- c("n per group", "N Total")
fmtcomputedss\$"n per group" <- makeCensored(fmtcomputedss\$"n per group",
nmax)
fmtcomputedss\$"N Total" <- makeCensored(fmtcomputedss\$"N Total",
ngrps*nmax)

cat(paste("Sample Size Table",
if(planningname!="") paste("for", planningname), "\n"))
if(is.expression(endptname) || (endptname!="")) {
cat(paste("Endpoint:", endptname, "\n"))
}

## Taken from base:::chartr help file
.simpleCap <- function(x) {
s <- strsplit(x, " ")[[1]]
paste(toupper(substring(s, 1,1)), substring(s, 2),
sep="", collapse=" ")
}

diffmetric <- "Differences"
if(endptscale=="log") {
diffmetric <- paste("Percent", .simpleCap(direction), diffmetric)
}
else {
diffmetric <- paste(.simpleCap(direction),  diffmetric)
}
cat(diffmetric, "\n")

cat(paste(round(100*power,0),"% Power and ",
round(100*alpha,0),
"% Significance Level\n", sep=""))

cat("Variability Estimate ",
if(endptscale=="log") {
"(Log scale) "
},
"of ",
signif(sigmaest, 4),
"\n",
ngrps, " Groups",
"\n",
sep="")

curwidth <- getOption("width")
if(curwidth < 500) { options(width=500) }

print(fmtcomputedss, quote=FALSE)

}
else if(display=="show") {
showDefault(computedss)
}
## else show nothing
invisible(computedss)
}

samplesizegraph <- function(sstable,
Nscale="log",
mmdscale="log",
difftype,
direction,
analysisname="",
endptname="",
alpha=0.05,
power=0.80,
sigmaest,
nmax,
Nticklabels=NULL,
mmdticklabels=NULL,
ylab=NULL, ylabright=NULL,
titlestamp=TRUE,
explanation=TRUE, ...) {
##
## PURPOSE: Graph Sample Size as a function of Difference
##
## NOTE: argument sstable needs to be a dataframe of proper format
##
difftype <- validArgMatch(difftype, c("percent","amount","simple"))
direction <- validArgMatch(direction, c("increasing","decreasing"))
Nscale <- validArgMatch(Nscale, c("log","original"))
mmdscale <- validArgMatch(mmdscale, c("log","original"))
validAlpha(alpha)
validBoolean(titlestamp)
validBoolean(explanation)
dots <- list(...)
validDotsArgs(dots, names=c(""))

if(length(mmdvec <- sstable[, 1])==1) {
stop(cgMessage("No graph is produced since only one",
"difference value was specified in the sample size",
"table"))
}

options(warn=-1)
curpar <- par(new=FALSE, mgp=c(3,0.25,0), tck=-0.010)
options(warn=0)
on.exit(par(curpar))

n <- sstable[, 2]
N <- sstable[, 3]
numberofgrps <- round((N/n)[1], 0)

N.scaled <- scaleVar(N, endptscale=Nscale)
N.logscale <- if(Nscale=="log") TRUE else FALSE

mmd.logscale <- if(mmdscale=="log") TRUE else FALSE
logandpct <- if(mmd.logscale &&  difftype=="percent") TRUE else FALSE
mmd.scaled <- scaleVar(mmdvec, endptscale=mmdscale,
percent=logandpct)

parmar <- par(mar=c(5, 4, 4, 3) + 0.1)
curpar\$mar <- parmar\$mar

plot(mmd.scaled, N.scaled, type="n", pch=1,
xlab="",
ylab=if(is.null(ylab)) paste("Total Sample Size from",
numberofgrps, "groups") else ylab,
axes=FALSE,
xlim=rangeExtend(range(mmd.scaled),
pct=list(minside=6, maxside=0)))
grid(lty=1)
lines(mmd.scaled, N.scaled, type="b", pch=1)

xlabchar <- paste(if(difftype=="percent") {
## "Percent Difference"
if(direction=="decreasing") {
"Percent REDUCTION"
}
else {
"Percent INCREASE"
}
}
else {
## "Simple Difference"
if(direction=="decreasing") {
"Simple REDUCTION"
}
else {
"Simple INCREASE"
}
},
": ", if(numberofgrps > 2) {"Minimum "} ,"Detectable",
" Difference", sep="")

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)
}

if(titlestamp) {
title(main=paste("Sample Size Graph\n"), line=2, cex.main=1.1)
}

if(explanation) {
mtext(side=3, line=0.25,
text=paste(round(100*power,0)," % Power ; ",
round(100*alpha,0),
" % Significance Level ; ",
"Variability Estimate ",
if(difftype=="percent") {
"(Log scale) "
},
"of ",
signif(sigmaest, 4),
sep=""),
cex=0.7)
}

## Filling in axes tick marks and such.
names(N.scaled) <- as.character(N)
N.tickmarks <- setupAxisTicks(N, logscale=N.logscale,
digits=0)

if(!is.null(Nticklabels)) {
validList(Nticklabels, names=c("mod","marks"),
argname="Nticklabels")
N.tickmarks <- makeTickMarks(Nticklabels, N.tickmarks)
}

n.ticklabels <- paste("", round(N.tickmarks/numberofgrps, 0))
ycex <- yTicksCex(N.tickmarks, cexinit=0.8)

axis(2, at=scaleVar(N.tickmarks, endptscale=Nscale),
axis(4, at=scaleVar(N.tickmarks, endptscale=Nscale),
minmaxTicks(N, logscale=N.logscale, digits=0)

mmdvec.ticks <- if(logandpct) {
pctToRatio(mmdvec)
}
else {
mmdvec
}

mmdDigits <- getNumDigits(mmdvec)
names(mmd.scaled) <- as.character(mmdvec)

mmd.tickmarks <- setupAxisTicks(
mmdvec.ticks,
axis="x",
logscale=mmd.logscale,
digits=mmdDigits,
ratio=logandpct,
percent=logandpct,
remticks=TRUE)

if(length(mmd.tickmarks) < 3) {
mmd.tickmarks <- setupAxisTicks(
mmdvec.ticks,
axis="x",
logscale=mmd.logscale,
digits=getNumDigits(mmdvec.ticks),
ratio=logandpct,
percent=logandpct,
remticks=TRUE)
}

mmd.tickmarks <- mmd.tickmarks[names(mmd.tickmarks)!="0"]

if(!is.null(mmdticklabels)) {
validList(mmdticklabels, names=c("mod","marks"),
argname="mmdticklabels")
mmd.tickmarks <- makeTickMarks(mmdticklabels, mmd.tickmarks,
percent=logandpct)
}

xcex <- xTicksCex(mmd.tickmarks, cexinit=0.8)
axis(1, at=scaleVar(mmd.tickmarks, endptscale=mmdscale),
las=1)

box()
mtext(side=4, text=if(is.null(ylabright)) paste("Sample size per each of the",
numberofgrps, "groups") else ylabright, line=1.8, cex=0.9, adj=0)

if(any(n >= nmax)) {
.R. <- TRUE
coordinates <- Hmisc::largest.empty(mmd.scaled, N.scaled,
width=0.30*(diff(range(
c(min(mmd.scaled),max(mmd.scaled))))),
height=0.25*(diff(range(
c(min(N.scaled),max(N.scaled))))),
method="area")
xjust=0, yjust=1, bty="n",
legend=paste("For at least one difference\n",
"the sample size\n",
"calculations were truncated\n",
"at",  paste(nmax, "per group."),
sep=" "),
text.col="red",
cex=0.70)
}

invisible()

}

setClass("cgOneFactorSampleSizeTable",
representation(ols.sstable="dataframeMatrixOrNULL",
rr.sstable="dataframeMatrixOrNULL",
settings="list"),
prototype(ols.sstable=NULL,
rr.sstable=NULL,
settings=list()))

setMethod("samplesizeTable", "cgOneFactorFit",
samplesizeTable.cgOneFactorFit <-
function(fit, direction, mmdvec, power=0.80, alpha=0.05,
nmax=1000, display="print", ...) {
##
## PURPOSE:  Compute sample sizes based on variance
## estimation from One Factor Fit and other specifications
##
## Input arguments check
if(class(fit@uvfit)[1]=="gls" || class(fit@aftfit)[1]=="survreg") {
stop(cgMessage("There is no samplesizeTable method",
"defined for a fitted model that allowed",
"unequal variances or censored observations."))
}
dots <- list(...)
validDotsArgs(dots, names=c("model","ngrps"))
direction <- validArgMatch(direction, c("increasing", "decreasing"))
validNumeric(mmdvec, positive=TRUE)
validAlpha(alpha)
validPower(power)
validNumeric(nmax, positive=TRUE, integer=TRUE)
if(nmax < 2) {
stop(cgMessage("The nmax argument needs to be an integer of 2 or greater."))
}
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"
}
ngrpsarg <- getDotsArgName(dots, "ngrps")
if(!is.na(ngrpsarg)) {
ngrps <- eval(parse(text=paste("dots\$", ngrpsarg, sep="")))
}
else {
ngrps <- 2
}

validNumeric(ngrps, positive=TRUE, integer=TRUE)
if(ngrps < 2) {
stop(cgMessage("The ngrps argument needs to be an integer of 2",
"or greater."))
}

## initializations
settings <- fit@settings
ols <- rr <- FALSE
ols.sstable <- rr.sstable <- NULL
ols.sigmaest <- rr.sigmaest <- NULL
##
endptscale <- settings\$endptscale
planningname <- settings\$analysisname
endptname <- settings\$endptname

rrfit <- fit@rrfit
olsfit <- fit@olsfit

if(class(rrfit)[1]=="rlm" && model!="olsonly") {
rr <- TRUE
}
if(class(olsfit)[1]=="lm" && model!="rronly") {
ols <- TRUE
if(!rr) model <- "olsonly"
}

dendf <- function(n, ngrps) ngrps*(n-1)

if(ols) {
ols.sigmaest <- summary(olsfit)\$sigma
ols.sstable <- samplesize(sigmaest=ols.sigmaest,
endptscale=endptscale,
planningname=planningname,
endptname=endptname,
ngrps=ngrps, direction=direction,
mmdvec=mmdvec,
dendf=dendf,
alpha=alpha, power=power,
nmax=nmax,
display="none")
}
if(rr) {
rr.sigmaest <- summary(rrfit)\$stddev
rr.sstable <- samplesize(sigmaest=rr.sigmaest,
endptscale=endptscale,
planningname=planningname,
endptname=endptname,
ngrps=ngrps, direction=direction,
mmdvec=mmdvec,
dendf=dendf,
alpha=alpha, power=power,
nmax=nmax,
display="none")
}

settings <- c(settings,
list(sigmaest=list(ols=ols.sigmaest, rr=rr.sigmaest),
planningname=planningname,
ngrps=ngrps,
direction=direction,
alpha=alpha, power=power, nmax=nmax)
)

returnObj <- new("cgOneFactorSampleSizeTable",
rr.sstable=rr.sstable,
ols.sstable=ols.sstable,
settings=settings)

if(display=="print") {
print(returnObj)
}
else if(display=="show") {
show(returnObj)
}
## else display=="none"
invisible(returnObj)
}
)

setMethod("print", "cgOneFactorSampleSizeTable",
print.cgOneFactorSampleSizeTable <-
function(x, title=NULL, endptname=NULL, ...) {
##
## PURPOSE: Semi-formatted print version of Sample Size Table
##
## NOTE: Had to use x as an argument because of the system defined
## generic. 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"
}

ols <- rr <- FALSE  ## Initializations
if(!is.null(rr.sstable <- object@rr.sstable) && model!="olsonly") {
rr <- TRUE
rr.sstable <- as.data.frame(rr.sstable)
}
if(!is.null(ols.sstable <- object@ols.sstable) && (model!="rronly")) {
ols <- TRUE
ols.sstable <- as.data.frame(ols.sstable)
}

settings <- object@settings
alpha <- settings\$alpha
power <- settings\$power
endptscale <- settings\$endptscale
ngrps <- settings\$ngrps
nmax <- settings\$nmax

if(is.null(title)) {
title <- paste("Sample Size Table from", settings\$planningname)
}
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")) }

## Taken from base:::chartr help file
.simpleCap <- function(x) {
s <- strsplit(x, " ")[[1]]
paste(toupper(substring(s, 1,1)), substring(s, 2),
sep="", collapse=" ")
}

diffmetric <- "Differences"
if(settings\$endptscale=="log") {
diffmetric <- paste("Percent", .simpleCap(settings\$direction), diffmetric)
}
else {
diffmetric <- paste(.simpleCap(settings\$direction),  diffmetric)
}
cat(diffmetric, "\n")

cat(paste(round(100*power,0),"% Power and ",
round(100*alpha,0),
"% Significance Level\n", sep=""))

fmtdig <- function(x) {
x\$mmd <- fround(x\$mmd, getNumDigits(x\$mmd))
rownames(x) <- x\$mmd
x\$mmd <- NULL
names(x) <- c("n per group", "N Total")
x\$"n per group" <- makeCensored(x\$"n per group",
nmax)
x\$"N Total" <- makeCensored(x\$"N Total",
ngrps*nmax)
x
}

informSettings <- function(sigmaest) {
cat("Variability Estimate ",
if(endptscale=="log") {
"(Log scale) "
},
"of ",
signif(sigmaest, 4),
"\n",
ngrps, " Groups",
"\n",
sep="")
}

curwidth <- getOption("width")
if(curwidth < 500) { options(width=500) }

if(ols) {
cat("\nClassical Least Squares Based\n")
informSettings(settings\$sigmaest\$ols)
print(fmtdig(ols.sstable), quote=FALSE)
}

if(rr) {
cat("\nResistant & Robust Based\n")
informSettings(settings\$sigmaest\$rr)
print(fmtdig(rr.sstable), quote=FALSE)
}

invisible()
})

setMethod("show", "cgOneFactorSampleSizeTable",
show.cgOneFactorSampleSizeTable <- function(object) showDefault(object))

setMethod("samplesizeGraph", "cgOneFactorSampleSizeTable",
samplesizeGraph.cgOneFactorSampleSizeTable <-
function(sstable, Nscale = "log", mmdscale = "log",
## cgtheme=TRUE, device="single",
...) {
##
## PURPOSE: create a graph of samplesize calculations
##
## Input arguments check
dots <- list(...)
validDotsArgs(dots, names=c("model",
"mmdticklabels","Nticklabels",
"cgtheme","device"))

settings <- sstable@settings
difftype <- if(settings\$endptscale=="log") "percent" else "simple"
alpha <- settings\$alpha
power <- settings\$power
planningname <- settings\$planningname
endptlabel <- makeEndptLabel(settings\$endptname, settings\$endptunits)
ngrps <- settings\$ngrps
direction <- settings\$direction
nmax <- settings\$nmax
stamps <- settings\$stamps
sigmaest <- settings\$sigmaest

rr.sstable <- sstable@rr.sstable
ols.sstable <- sstable@ols.sstable

Nticklabelsarg <- getDotsArgName(dots, "Nticklabels")
if(!is.na(Nticklabelsarg)) {
Nticklabels <- eval(parse(text=paste("dots\$", Nticklabelsarg, sep="")))
validList(Nticklabels, names=c("mod","marks"),
argname="Nticklabels")
}
else {
Nticklabels <- NULL
}

mmdticklabelsarg <- getDotsArgName(dots, "mmdticklabels")
if(!is.na(mmdticklabelsarg)) {
mmdticklabels <- eval(parse(text=paste("dots\$", mmdticklabelsarg, sep="")))
validList(mmdticklabels, names=c("mod","marks"),
argname="mmdticklabels")
}
else {
mmdticklabels <- 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"
}

cgthemearg <- getDotsArgName(dots, "cgtheme")
if(!is.na(cgthemearg)) {
cgtheme <- eval(parse(text=paste("dots\$", cgthemearg, sep="")))
validBoolean(cgtheme)
}
else {
dots\$cgtheme <- TRUE
cgtheme <- TRUE
}

devicearg <- getDotsArgName(dots, "device")
if(!is.na(devicearg)) {
device <- eval(parse(text=paste("dots\$", devicearg, sep="")))
}
else {
dots\$device <- "single"
device <- "single"
}

ols <- rr <- FALSE  ## Initializations
if(!is.null(rr.sstable) && model!="olsonly") {
rr <- TRUE
}
if(!is.null(ols.sstable) && model!="rronly") {
ols <- TRUE
}

thetitle <- "Sample Size Graph"

if(rr && ols && is.element(model, "both") && device=="single") {
all.dfr <- as.data.frame(rbind(ols.sstable, rr.sstable))
all.dfr\$typef <- factorInSeq(c(rep("Classical",
nrow(ols.sstable)),
rep("Resistant & Robust",
nrow(rr.sstable))))
thetitle <- paste(thetitle, "s", sep="")

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)),
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)),
trellisparstg3 <- trellis.par.get("axis.components")
trellis.par.set("axis.components",
))
on.exit(trellis.par.set("axis.components", trellisparstg3),

all.N <- all.dfr\$N
all.mmd <- all.dfr\$mmd

N.logscale <- if(Nscale=="log") TRUE else FALSE
all.N.scaled <- scaleVar(all.N, endptscale=Nscale)

mmd.logscale <- if(mmdscale=="log") TRUE else FALSE
logandpct <- if(mmd.logscale &&  difftype=="percent") TRUE else FALSE
all.mmd.scaled <- scaleVar(all.mmd, endptscale=mmdscale,
percent=logandpct)

grid.newpage()
lvp <- viewport(x=0, width=unit(1, "npc") - unit(2, "lines"),
name="lvp", just="left")
tvp <- viewport(x=1, width=unit(2, "lines"),
name="tvp", just="right")

thegraph <- xyplot(all.N.scaled  ~ all.mmd.scaled | typef,
data=all.dfr,
panel=function(x, y, ...) {
panel.grid(h=-1, v=-1)
panel.xyplot(x, y, type="b", pch=1,
col="black", ...)

N.tickmarks <- setupAxisTicks(all.N,
logscale=N.logscale,
digits=0,
grid=TRUE,
ycex=0.7)
if(!is.null(Nticklabels)) {
validList(Nticklabels, names=c("mod","marks"),
argname="Nticklabels")
N.tickmarks <- makeTickMarks(Nticklabels, N.tickmarks)
}

n.ticklabels <- paste("", round(N.tickmarks/ngrps, 0))
ycex <- yTicksCex(N.tickmarks,
cexinit=0.7,
cexthreshold = 0.7,
grid=TRUE)

all.mmd.ticks <- if(logandpct) {
pctToRatio(all.mmd)
}
else {
all.mmd
}

mmdDigits <- getNumDigits(all.mmd)
mmd.tickmarks <-
setupAxisTicks(
all.mmd.ticks,
axis="x",
logscale=mmd.logscale,
digits=mmdDigits,
ratio=logandpct,
percent=logandpct,
remticks=TRUE,
grid=TRUE,
xcex=0.70)

if(length(mmd.tickmarks) < 3) {
mmd.tickmarks <- setupAxisTicks(
all.mmd.ticks,
axis="x",
logscale=mmd.logscale,
digits=getNumDigits(all.mmd.ticks),
ratio=logandpct,
percent=logandpct,
remticks=TRUE,
grid=TRUE,
xcex=0.70)
}

mmd.tickmarks <- mmd.tickmarks[names(mmd.tickmarks)!="0"]

if(!is.null(mmdticklabels)) {
validList(mmdticklabels, names=c("mod","marks"),
argname="mmdticklabels")
mmd.tickmarks <- makeTickMarks(mmdticklabels, mmd.tickmarks,
percent=logandpct)
}
xcex <- xTicksCex(mmd.tickmarks,
cexinit=0.7,
cexthreshold = 0.7,
grid=TRUE)

if(panel.number()==1) {
panel.axis(side="left",
at=scaleVar(N.tickmarks,
endptscale=Nscale),
labels=names(N.tickmarks),
tck=0.15, text.cex=ycex,
rot=0,
outside=TRUE)
}
else if(panel.number()==2) {
panel.axis(side="right",
at=scaleVar(N.tickmarks,
endptscale=Nscale),
labels=n.ticklabels,
tck=0.15, text.cex=ycex,
rot=0,
outside=TRUE)
}

panel.axis(side="bottom",
at=scaleVar(mmd.tickmarks,
mmdscale),
labels=names(mmd.tickmarks),
tck=0.15, rot=0,
text.cex=xcex,
outside=TRUE)
N.range <- range(y)
if(N.logscale) {
N.range <- round(10^N.range, 0)
}
minmaxTicks(N.range,
theaxis="y", logscale=N.logscale,
digits=0,
grid=TRUE)
if(any(round(N.range/ngrps, 0) >= nmax)) {

.R. <- TRUE
coordinates <-  Hmisc::largest.empty(x, y,
width=0.30*(diff(range(
c(min(x),
max(x))))),
height=0.25*(diff(range(
c(min(y),
max(y))))),
grid=TRUE,
method="area")
panel.text(x=coordinates\$rect\$x[4], y=coordinates\$rect\$y[4],
label=paste("For at least one difference\n",
"the sample size\n",
"calculations were truncated\n",
"at",  paste(nmax, "per group."),
sep=" "),
col="red",
cex=0.6,
}
},
layout=c(2,1), aspect=1,
as.table=TRUE,
xlim=rangeExtend(range(all.mmd.scaled),
pct=list(minside=20, maxside=10)),
xlab=paste(if(difftype=="percent") {
## "Percent Difference"
if(direction=="decreasing") {
"Percent REDUCTION"
}
else {
"Percent INCREASE"
}
}
else {
## "Simple Difference"
if(direction=="decreasing") {
"Simple REDUCTION"
}
else {
"Simple INCREASE"
}
},
": ", if(ngrps > 2) {"Minimum "} ,"Detectable",
" Difference",
sep=""),
ylab=paste("Total Sample Size from",
ngrps, "groups"),
scales=list(y=list(labels=NULL, tck=0),
x=list(labels=NULL, tck=0)),
main=list(label=paste(thetitle, "\n",
planningname,
sep=""), cex=1.1),
par.strip.text=list(cex=0.7)
)

##on.exit(lattice:::lattice.setStatus(print.more = FALSE),
pushViewport(lvp)
print(thegraph, newpage=FALSE)
upViewport()
pushViewport(tvp)
grid.text(paste("Sample size per each of the",
ngrps, "groups"),
x = unit(1, "lines"), rot=90, gp=gpar(cex=0.75),
just="center")

seekViewport(trellis.vpname("xlab"))
grid.text(catCharExpr("in ", endptlabel), y = unit(-1, "lines"))
upViewport(0)

seekViewport(trellis.vpname("main"))
grid.text(paste(round(100*power,0)," % Power ; ",
round(100*alpha,0),
" % Significance Level ; ",
"Variability Estimate ",
if(difftype=="percent") {
"(Log scale) "
},
"of\n",
signif(sigmaest\$ols, 4),
" for Classical and ",
signif(sigmaest\$rr, 4),
" for Resistant & Robust",
sep=""),
y=unit(-1, "lines"), rot=0, gp=gpar(cex=0.60),
just="center")
upViewport(0)
usedgrid <- TRUE
}

else if(((model=="olsonly" && ols) || (ols && !rr && model=="both")) &&
device=="single") {
samplesizegraph(ols.sstable,
mmdscale,
Nscale,
difftype,
direction,
planningname,
endptlabel,
alpha,
power,
sigmaest=sigmaest\$ols,
nmax=nmax,
Nticklabels,
mmdticklabels,
titlestamp=FALSE,
explanation=FALSE)
title(line=2,
main=paste("Sample Size Graph\n",
planningname, sep=""), cex=1.1)
mtext(side=3, line=0.5,
text=paste(round(100*power,0)," % Power ; ",
round(100*alpha,0),
" % Significance Level ; ",
"Classical Variability Estimate ",
if(difftype=="percent") {
"(Log scale) "
},
"of ",
signif(sigmaest\$ols, 4),
sep=""), cex=0.60)
usedgrid <- FALSE

}

else if(((model=="rronly" && rr) || (!ols && rr && model=="both")) &&
device=="single") {
samplesizegraph(rr.sstable,
mmdscale,
Nscale,
difftype,
direction,
planningname,
endptlabel,
alpha,
power,
sigmaest=sigmaest\$rr,
nmax=nmax,
Nticklabels,
mmdticklabels,
titlestamp=FALSE,
explanation=FALSE)
title(line=2,
main=paste("Sample Size Graph\n",
planningname, sep=""), cex=1.1)
mtext(side=3, line=0.5,
text=paste(round(100*power,0)," % Power ; ",
round(100*alpha,0),
" % Significance Level ; ",
"Resistant & Robust Variability Estimate ",
if(difftype=="percent") {
"(Log scale) "
},
"of ",
signif(sigmaest\$rr, 4),
sep=""), cex=0.60)
usedgrid <- FALSE

}
else if(rr && ols &&
is.element(model, "both")) {

}

samplesizegraph(ols.sstable,
mmdscale,
Nscale,
difftype,
direction,
planningname,
endptlabel,
alpha,
power,
sigmaest=sigmaest\$ols,
nmax=nmax,
Nticklabels,
mmdticklabels,
titlestamp=FALSE,
explanation=FALSE)
title(line=2,
main=paste("Sample Size Graph\n",
planningname, sep=""), cex=1.1)
mtext(side=3, line=0.5,
text=paste(round(100*power,0)," % Power ; ",
round(100*alpha,0),
" % Significance Level ; ",
"Classical Variability Estimate ",
if(difftype=="percent") {
"(Log scale) "
},
"of ",
signif(sigmaest\$ols, 4),
sep=""), cex=0.60)
if(stamps) graphStampCG(grid=FALSE)

if(device=="multiple") {
## since we only want dots arguments for trellis.device in
## next call we need some housekeeping
if(!is.null(dots\$model)) dots\$model <- NULL
if(!is.null(dots\$device) &&
dots\$device<- NULL
}
if(!is.null(dots\$cgtheme)) dots\$cgtheme <- NULL
if(!is.null(dots\$mmdticklabels)) dots\$mmdticklabels <- NULL
if(!is.null(dots\$Nticklabels)) dots\$Nticklabels <- NULL

do.call("cgDevice", c(list(new=TRUE), dots))
cat(cgMessage("A new graphics device has been generated",
"to hold",
"SampleSize graph version based on",
"the Resistant & Robust estimate.",
"The version based on the Classical Least Squares",
"estimate is on the previous",
"device.\n",
warning=TRUE))
}
samplesizegraph(rr.sstable,
mmdscale,
Nscale,
difftype,
direction,
planningname,
endptlabel,
alpha,
power,
sigmaest=sigmaest\$rr,
nmax=nmax,
Nticklabels,
mmdticklabels,
titlestamp=FALSE,
explanation=FALSE)
title(line=2,
main=paste("Sample Size Graph\n",
planningname, sep=""), cex=1.1)
mtext(side=3, line=0.5,
text=paste(round(100*power,0)," % Power ; ",
round(100*alpha,0),
" % Significance Level ; ",
"Resistant & Robust Variability Estimate ",
if(difftype=="percent") {
"(Log scale) "
},
"of ",
signif(sigmaest\$rr, 4),
sep=""), cex=0.60)
usedgrid <- FALSE
}
else {
stop(cgMessage("The chosen device and model arguments",
"are not compatible either with each other",
"or with the fitted model(s) in the SampleSize table object.",
seeHelpFile("SampleSizeTable")))
}

if(stamps) graphStampCG(grid=usedgrid)
invisible()

}
)

validDenDf <- function(dendf) {
if(!is.function(dendf)) {
stop(cgMessage("The object needs to be a function",
"with arguments n and ngrps."))
}
if(!identical(sort(names(formals(dendf))),
c("n","ngrps"))) {
stop(cgMessage("The dendf function needs to include",
"the exact argument names of n and ngrps."))
}
return(TRUE)
}

validNcp <- function(ncp) {
if(is.null(ncp)) return(TRUE)
if(!is.function(ncp)) {
stop(cgMessage("The object needs to be a function."))
}
theargnames <- names(formals(ncp))
if(sum(theargnames %in% c("sigmaest", "mmdvec", "n")) < 3) {
stop(cgMessage("The ncp function needs to include",
"the argument names of sigmaest, mmdvec, n."))
}
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 30, 2017, 12:12 a.m.