summaryS <- function(formula, fun=NULL,
data=NULL, subset=NULL, na.action=na.retain,
continuous=10,
...) {
formula <- Formula::Formula(formula)
Y <- if(length(subset))
model.frame(formula, data=data, subset=subset, na.action=na.action)
else
model.frame(formula, data=data, na.action=na.action)
X <- Formula::model.part(formula, data=Y, rhs=1)
Y <- Formula::model.part(formula, data=Y, lhs=1)
nY <- NCOL(Y)
nX <- NCOL(X)
namY <- names(Y)
namX <- names(X)
if(nX == 0) X <- data.frame(x=rep(1, NROW(Y)))
ux <- unique(X)
Z <- NULL
n <- nrow(X)
g <- function(y) {
y <- y[! is.na(y)]
if(is.character(y)) c(NA, NA)
else if(is.factor(y)) c(1, length(levels(y)))
else if(length(unique(y)) < continuous) range(y)
else if(! is.matrix(y)) quantile(y, c(0.01, 0.99))
}
ylim <- lapply(Y, g)
w <- reshape(cbind(X, Y), direction='long', v.names='y',
varying=namY, times=namY, timevar='yvar')
if(inherits(Y[[1]], 'Surv')) {
at <- attributes(Y[[1]])
at <- at[setdiff(names(at), c('dim', 'dimnames'))]
attributes(w$y) <- c(attributes(w$y), at)
}
w$yvar <- factor(w$yvar, namY)
funlabel <- NULL
if(length(fun)) {
by <- c('yvar', if(length(namX)) namX else 'x')
y <- w$y
w <- summarize(y, w[by], fun, type='matrix', keepcolnames=TRUE)
funlabel <- if(is.matrix(w$y)) colnames(w$y)[1]
}
gg <- function(x) if(is.character(x) || is.factor(x))
'categorical' else 'numeric'
## For some reason sapply is doubling names e.g. sbp.sbp
sapply2 <- function(data, ...) {
s <- sapply(data, ...)
names(s) <- names(data)
s
}
xlabels <- sapply2(X, label)
xlabels <- ifelse(xlabels == '', names(xlabels), xlabels)
ylabels <- sapply2(Y, label)
ylabels <- ifelse(ylabels == '', names(ylabels), ylabels)
for(n in names(w)) # takes care of R change stringsAsFactors=FALSE
if(is.character(w[[n]])) w[[n]] <- factor(w[[n]])
structure(w, class=c('summaryS', 'data.frame'),
formula=formula, fun=fun,
xnames=names(X), xlabels=xlabels, xunits=sapply2(X, units),
xtype=sapply(X, gg),
ynames=namY, ylabels=ylabels, yunits=sapply2(Y, units),
ylim=ylim, funlabel=funlabel)
}
plot.summaryS <-
function(x, formula=NULL, groups=NULL, panel=NULL,
paneldoesgroups=FALSE, datadensity=NULL, ylab='',
funlabel=NULL, textonly='n', textplot=NULL, digits=3, custom=NULL,
xlim=NULL, ylim=NULL, cex.strip=1, cex.values=0.5, pch.stats=NULL,
key=list(columns=length(groupslevels),
x=.75, y=-.04, cex=.9,
col=lattice::trellis.par.get('superpose.symbol')$col,
corner=c(0,1)),
outerlabels=TRUE, autoarrange=TRUE, scat1d.opts=NULL, ...)
{
sRequire('lattice')
sRequire('latticeExtra')
xtype <- attr(x, 'xtype')
nn <- sum(xtype == 'numeric')
if(nn > 1) stop('does not handle more than one numeric continuous x')
X <- x
at <- attributes(x)
Form <- at$formula
nX <- at$nX
nY <- at$nY
ylabels <- at$ylabels
yunits <- at$yunits
ylims <- at$ylim
xnames <- at$xnames
xlabels <- at$xlabels
xunits <- at$xunits
fun <- at$fun
funlabel <- if(length(at$funlabel)) at$funlabel else funlabel
Panel <- panel
ptype <- if(length(fun)) { # used to always be 'dot'
if(length(Panel)) 'xy.special' else 'dot'
} else 'xy'
if(ptype %in% c('xy', 'xy.special') && ! any(xtype == 'numeric'))
stop('must have a numeric x variable to make x-y plot')
groupslevels <- if(length(groups)) levels(x[[groups]])
condvar <- xnames[xtype == 'categorical']
## Reorder condvar in descending order of number of levels
numu <- function(x)
if(is.factor(x)) length(levels(x))
else
length(unique(x[! is.na(x)]))
if(autoarrange && length(condvar) > 1) {
nlev <- sapply(X[condvar], numu)
condvar <- condvar[order(nlev)]
}
form <- if(length(formula)) formula
else {
## Non-groups conditioning variables
ngcond <- setdiff(condvar, groups)
## Collapsed non-group conditioning variables
ccv <- paste('|', paste(c(ngcond, 'yvar'), collapse=' * '))
## Collapsed non-group cond var after the first
ccv1 <- if(length(ngcond) > 1)
paste('|', paste(c(ngcond[-1], 'yvar'), collapse=' * '))
f <- if(ptype %in% c('xy', 'xy.special'))
paste('y ~', xnames[xtype == 'numeric'], ccv, sep='')
else paste(ngcond[1], '~ y', ccv1, sep='')
as.formula(f)
}
pst <- list(cex=cex.strip)
yvarlev <- NULL
for(v in levels(X$yvar)) {
un <- yunits[v]
l <- if(ylabels[v] == v && un == '') v else
labelPlotmath(ylabels[v], un)
yvarlev <- c(yvarlev, l)
}
strip <- function(which.given, which.panel, var.name, factor.levels, ...) {
current.var <- var.name[which.given]
levs <- if(current.var == 'yvar') yvarlev else factor.levels
lattice::strip.default(which.given, which.panel, var.name,
factor.levels=levs, ...)
}
ylev <- levels(X$yvar)
## lims <- if(length(xlim)) xlim else ylims[ylev]
lims <- ylims[ylev]
## lims needs to be repeated according to layout
vars <- all.vars(form)
cond <- vars[- (1 : 2)]
if(! all(cond %in% 'yvar') && cond[1] != 'yvar') {
ngny <- setdiff(cond, 'yvar')
nr <- length(unique(do.call('paste', X[ngny])))
lims <- rep(lims, each=nr)
}
if(length(ylim)) lims <- ylim
d <- if(ptype == 'xy') {
pan <- if(! length(datadensity)) function(...) {}
else
function(x, y, subscripts, groups=NULL, ...) {
gp <- length(groups)
plot.line <-
lattice::trellis.par.get(if(gp) "superpose.line" else "plot.line")
col <- plot.line$col
gr <- if(gp) groups[subscripts] else factor(rep('', length(x)))
ypos <- unit(1, 'npc')
for(i in 1:length(levels(gr))) {
j <- which(gr == levels(gr)[i])
ypos <- ypos -
if(length(scat1d.opts) && 'nhistSpike' %in% names(scat1d.opts))
unit(2, 'mm') else unit(1, 'mm')
do.call('scat1d', c(list(x=x[j], y=ypos, col=col[i], grid=TRUE),
scat1d.opts))
}
}
scal <- list(y=list(relation='free', limits=lims, rot=0))
if(length(xlim)) scal$x <- list(limits=xlim)
xlab <- labelPlotmath(xlabels[xtype == 'numeric'],
xunits [xtype == 'numeric'])
if(! length(groups)) {
if(! length(Panel)) Panel <- lattice::panel.xyplot
lattice::xyplot(form, data=X,
panel=function(...) {Panel(...); pan(...)},
xlab=xlab, ylab=ylab, scales=scal, strip=strip,
par.strip.text=pst, ...)
} else {
panel.groups <- if(paneldoesgroups) NULL else
if(length(Panel)) Panel else lattice::panel.xyplot
if(! paneldoesgroups) Panel <- lattice::panel.superpose
g <- if(length(panel.groups))
"lattice::xyplot(form, groups=%s, data=X, scales=scal, panel=function(...) {if(length(Panel)) Panel(..., scat1d.opts=scat1d.opts); pan(...)}, panel.groups=panel.groups, auto.key=key, xlab=xlab, ylab=ylab, strip=strip, par.strip.text=pst, ...)"
else
"lattice::xyplot(form, groups=%s, data=X, scales=scal, panel=function(...) {if(length(Panel)) Panel(..., scat1d.opts=scat1d.opts); pan(...)}, auto.key=key, xlab=xlab, ylab=ylab, strip=strip, par.strip.text=pst, ...)"
eval(parse(text=sprintf(g, groups)))
}
} else { # Dot chart or xy.special
## If > 1 calculated statistics, y is a matrix.
## Save first column as y and save remaining columns that are not
## text only as yother
y <- X$y
yother <- NULL
yText <- NULL
if(is.matrix(y) && length(c(textonly, textplot))) {
ptext <- colnames(y) %in% c(textonly, textplot)
if(any(ptext)) yText <- y[, ptext, drop=FALSE]
ponly <- colnames(y) %in% textonly
if(any(ponly)) y <- y[, ! ponly]
}
if(is.matrix(y) && ncol(y) > 1) {
yother <- y[, -1, drop=FALSE]
X$y <- y[, 1]
}
# ylev <- levels(X$yvar)
# lims <- if(length(xlim)) xlim else ylims[ylev]
# ## lims needs to be repeated according to layout
# vars <- all.vars(form)
# cond <- vars[- (1 : 2)]
# if(! all(cond %in% 'yvar') && cond[1] != 'yvar') {
# ngny <- setdiff(cond, 'yvar')
# nr <- length(unique(do.call('paste', X[ngny])))
# lims <- rep(lims, each=nr)
# }
# if(length(ylim)) lims <- ylim
scal <- list(x=list(relation='free', limits=xlim)) # limits=lims))
if(ptype == 'xy.special') names(scal) <- 'y'
if(ptype == 'dot') {
pan <- function(x, y, subscripts, groups=NULL, yother=NULL,...) {
gp <- length(groups)
dot.symbol <-
lattice::trellis.par.get(if(gp)'superpose.symbol' else 'dot.symbol')
plot.line <-
lattice::trellis.par.get(if(gp) "superpose.line" else "plot.line")
pch = dot.symbol$pch
col = dot.symbol$col
cex = dot.symbol$cex
font = dot.symbol$font
segmnts <- function (x0, y0, x1, y1, ...)
grid::grid.segments(x0, y0, x1, y1, default.units = "native",
gp = grid::gpar(...))
if(length(yother)) {
snames <- colnames(yother)
nc <- ncol(yother)
yoth <- yother[subscripts,, drop=FALSE]
gr <- if(gp) groups[subscripts] else factor(rep('', length(x)))
if(length(yText)) {
yText <- yText[subscripts,, drop=FALSE]
if(length(custom)) {
k <- custom(yText)
pasted <- k$result
llong <- unit(1, 'strwidth', k$longest)
} else {
pasted <- rep('', length(y))
for(i in 1 : ncol(yText)) {
if(i > 1) pasted <- paste(pasted, ' ', sep='')
pasted <- paste(pasted, colnames(yText)[i], '=',
format(round(yText[, i], digits=digits)),
sep='')
}
llong <- unit(1, 'strwidth',
paste(' ', pasted[which.max(nchar(pasted))],
sep=''))
}
xpos <- unit(1, 'npc') - unit(1, 'mm') - length(levels(gr)) * llong
}
for(i in 1:length(levels(gr))) {
j <- which(gr == levels(gr)[i])
if(nc > 1)
segmnts(yoth[j, 1], y[j], yoth[j, nc], y[j],
lwd=2 * plot.line$lwd[1],
lty=plot.line$lty[i], col=plot.line$col[i])
if(nc == 4)
segmnts(yoth[j ,2], y[j], yoth[j ,3], y[j],
lwd=2*plot.line$lwd[1],
lty=plot.line$lty[i], col=plot.line$col[i])
for(k in 1 : nc) {
if(length(pch.stats)) {
p <- pch.stats[snames[k]]
if(!is.na(p))
lattice::lpoints(yoth[j, k], y[j], pch=p, cex=cex, col=col[i],
font=font)
}
}
## Show selected statistics just under dot lines
if(length(yText)) {
xpos <- xpos + llong
grid::grid.text(pasted[j], xpos,
unit(y[j], 'native') - unit(1.75, 'mm'), just='right',
gp=grid::gpar(cex=cex.values, col=col[i]))
}
}
if(gp) lattice::panel.superpose(x, y, groups=as.numeric(groups),
subscripts=subscripts,
pch=pch, col=col, cex=cex, font=font, ...)
else
lattice::panel.dotplot(x, y, subscripts=subscripts,
pch=pch, col=col, cex=cex, font=font, ...)
}
else {
if(gp)
lattice::panel.superpose(x, y, groups=as.numeric(groups),
subscripts=subscripts,
pch=pch, col=col, cex=cex,
font=font, ...)
else
lattice::panel.dotplot(x, y, subscripts=subscripts,
pch=pch, col=col, cex=cex, font=font, ...)
}
}
d <- if(!length(groups))
lattice::dotplot(form, data=X, panel=pan, strip=strip, par.strip.text=pst,
xlab=funlabel, scale=scal, yother=yother,...)
else eval(parse(text=
sprintf("lattice::dotplot(form, groups=%s, data=X, panel=pan, strip=strip, par.strip.text=pst, auto.key=key, xlab=funlabel, scale=scal, yother=yother, ...)", groups) ))
} # end ptype 'dot'
else { # ptype 'xy.special'
xl <- labelPlotmath(xlabels[1], xunits[1])
yl <- if(ylab == '') funlabel else ylab
d <- if(!length(groups))
lattice::xyplot(form, data=X, panel=Panel, strip=strip,
par.strip.text=pst, xlab=xl, ylab=yl, scale=scal, yother=yother,
...)
else {
panel.groups <- if(paneldoesgroups) NULL else Panel
if(! paneldoesgroups) panel <- lattice::panel.superpose
g <- if(length(panel.groups))
"lattice::xyplot(form, groups=%s, data=X, panel=Panel, panel.groups=panel.groups, strip=strip, par.strip.text=pst, auto.key=key, xlab=xl, ylab=yl, scale=scal, yother=yother, ...)"
else "lattice::xyplot(form, groups=%s, data=X, panel=Panel, strip=strip, par.strip.text=pst, auto.key=key, xlab=xl, ylab=yl, scale=scal, yother=yother, ...)"
eval(parse(text=sprintf(g, groups)))
}
}
}
if(outerlabels && length(dim(d)) == 2)
d <- latticeExtra::useOuterStrips(d, strip=strip, strip.left=strip)
d
}
plotp.summaryS <-
function(data, formula=NULL,
groups=NULL, sfun=NULL, fitter=NULL,
showpts=! length(fitter),
funlabel=NULL, digits=5,
xlim=NULL, ylim=NULL,
shareX=TRUE, shareY=FALSE,
autoarrange=TRUE,
...)
{
xtype <- attr(data, 'xtype')
nn <- sum(xtype == 'numeric')
if(nn > 1) stop('does not handle more than one numeric continuous x')
X <- data
at <- attributes(data)
Form <- at$formula
nX <- at$nX
nY <- at$nY
ylabels <- at$ylabels
yunits <- at$yunits
ylims <- at$ylim
xnames <- at$xnames
xlabels <- at$xlabels
xunits <- at$xunits
fun <- at$fun
funlabel <- if(! length(funlabel) && length(at$funlabel))
at$funlabel else funlabel
funlabel <- htmlTranslate(funlabel, greek=TRUE)
ly <- length(ylabels)
ylab <- ylabels
for(i in 1 : length(ylab))
ylab[i] <- labelPlotmath(ylabels[i], yunits[i], html=TRUE)
aform <- function(n) as.formula(paste('~', n))
fmt <- function(x) htmlSN(x, digits=digits)
ptype <- if(length(fun)) { # used to always be 'dot'
if(length(sfun)) 'xy.special' else 'dot'
} else 'xy'
if(length(sfun)) ptype <- 'xy.special'
if(ptype %in% c('xy', 'xy.special') && ! any(xtype == 'numeric'))
stop('must have a numeric x variable to make x-y plot')
groupslevels <- if(length(groups)) levels(data[[groups]])
condvar <- xnames[xtype == 'categorical']
## Reorder condvar in descending order of number of levels
numu <- function(x)
if(is.factor(x)) length(levels(x))
else
length(unique(x[! is.na(x)]))
if(autoarrange && length(condvar) > 1) {
nlev <- sapply(X[condvar], numu)
condvar <- condvar[order(nlev)]
}
form <- if(length(formula)) formula
else {
## Non-groups conditioning variables
ngcond <- setdiff(condvar, groups)
## Collapsed non-group conditioning variables
ccv <- paste('|', paste(c(ngcond, 'yvar'), collapse=' * '))
## Collapsed non-group cond var after the first
ccv1 <- if(length(ngcond) > 1)
paste('|', paste(c(ngcond[-1], 'yvar'), collapse=' * '))
f <- if(ptype %in% c('xy', 'xy.special'))
paste('y ~', xnames[xtype == 'numeric'], ccv, sep='')
else paste(ngcond[1], '~ y', ccv1, sep='')
as.formula(f)
}
yvarlev <- NULL
for(v in levels(X$yvar)) {
un <- yunits[v]
l <- if(ylabels[v] == v && un == '') v else
labelPlotmath(ylabels[v], un, html=TRUE)
yvarlev <- c(yvarlev, l)
}
ylev <- levels(X$yvar)
lims <- ylims[ylev]
vars <- all.vars(form)
cond <- vars[- (1 : 2)]
if(! all(cond %in% 'yvar') && cond[1] != 'yvar') {
ngny <- setdiff(cond, 'yvar')
nr <- length(unique(do.call('paste', X[ngny])))
}
if(length(ylim)) lims <- ylim
gp <- length(groups)
ylev <- levels(X$yvar)
nyvar <- length(ylev)
xn <- setdiff(xnames, groups)
if(length(xn) %nin% 1:2)
stop(paste('expecting at most two variables, found these:',
paste(xn, collapse=', ')))
if(length(xn) > 1) {
strata <- aform(xn[2])
xn <- xn[1]
}
else
strata <- NULL
statnames <- if(is.matrix(X$y)) colnames(X$y)
.txt. <- paste0(xn, ': ', fmt(X[[xn]]))
if(gp) .txt. <- paste0(.txt., '<br>', X[[groups]])
nstat <- length(statnames)
if(nstat == 0) .txt. <- paste0(.txt., '<br>',
X$yvar, ': ', fmt(X$y))
else
for(i in 1 : nstat) { ## ?? was if(i > 1)'<br>'
if(i == 2 && length(funlabel) && funlabel != '' && funlabel != ' ')
.txt. <- paste0(.txt., '<br>', funlabel)
.txt. <- paste0(.txt., '<br>', statnames[i], ': ',
fmt(X$y[, i]))
}
X$.txt. <- .txt.
xlab <- labelPlotmath(xlabels[xn], xunits[xn], html=TRUE)
gp <- length(groups)
gr <- if(gp) X[[groups]] else factor(rep('', nrow(X)))
if(ptype == 'xy') {
if(nstat > 0 && ! gp) X <- cbind(X, tracename=statnames[1])
p <- plotlyM(X, x = aform(xn), y = ~y, htext = ~.txt.,
color = if(gp) aform(groups),
multplot = ~yvar, strata = strata,
xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim,
fitter=fitter, showpts=showpts,
...)
} else { ## end ptype xy
## Dot chart or xy.special
## If > 1 calculated statistics, y is a matrix.
## Save first column as y and save remaining columns as yother
y <- X$y
yother <- NULL
if(is.matrix(y) && ncol(y) > 1) {
yother <- y[, -1, drop=FALSE]
X$y <- y[, 1]
}
if(ptype == 'dot') {
R <- X
R$yhi <- NA
if(length(statnames) && ! gp)
R <- cbind(R, tracename=statnames[1])
if(length(yother)) {
snames <- colnames(yother)
nc <- ncol(yother)
## if(nc > 1) yother[,1] yother[,nc]
if((all(c('Lower', 'Upper') %in% snames)) && nc < 4) {
S <- R
S$y <- yother[, 'Lower']
S$yhi <- yother[, 'Upper']
if(! gp) S$tracename <- 'C.L.'
R <- rbind(R, S)
}
if(nc == 4) {
S <- R
S$y <- yother[, 2]
S$yhi <- yother[, 3]
if(! gp) S$tracename <- paste(snames[2:3], collapse=' ')
R <- rbind(R, S)
}
# for(k in 1 : nc) {
# if(length(pch.stats)) {
# p <- pch.stats[snames[k]]
# if(!is.na(p))
# R <- rbind(R, data.frame(x=X$y, y=yother[, k], yhi=NA, gr=gr,
# tracename=snames[k]))
# ## TODO: need to communicate pch=p
# }
# }
}
p <- plotlyM(R, x=aform(xn), multplot=~ yvar,
color = if(gp) aform(groups),
htext = ~ .txt.,
rotate = TRUE,
ylab = ylab,
xlim = xlim, ylim=ylim,
shareX = shareX, shareY=shareY, ...)
} # end ptype 'dot'
else { # ptype 'xy.special'
xl <- labelPlotmath(xlabels[1], xunits[1], html=TRUE)
yl <- if(! length(ylab) || ylab[1] == '') funlabel else ylab
p <- sfun(X[[xn]], X$y, groups=if(gp) gr, yother=yother, yvar=X$yvar,
maintracename=statnames[1], xlab=xl, ylab=yl,
xlim=xlim, ylim=ylim, zeroline=FALSE, ...)
}
}
p
}
mbarclPanel <- function(x, y, subscripts, groups=NULL, yother, ...) {
sRequire('lattice')
gp <- length(groups)
plot.line <-
lattice::trellis.par.get(if(gp) "superpose.line" else "plot.line")
yother <- yother[subscripts, , drop=FALSE]
se <- if('se' %in% colnames(yother)) yother[, 'se']
yother <- yother[, colnames(yother) %nin% c('n', 'se'), drop=FALSE]
if(all(c('0.375', '0.625') %in% colnames(yother))) {
## If HD median estimate is not between 0.375 and 0.625 quantiles
## take it to be the closer of the two
y375 <- yother[, '0.375']
y625 <- yother[, '0.625']
y <- pmin(y, y625)
y <- pmax(y, y375)
}
yother <- cbind(y=y, yother)
gr <- if(gp) groups[subscripts] else factor(rep('', length(x)))
lev <- levels(gr)
for(i in 1 : length(lev)) {
j <- which(gr == levels(gr)[i])
multLines(x[j], yother[j,, drop=FALSE],
col=plot.line$col[i],
lty=plot.line$lty[i],
lwd=plot.line$lwd[i],
grid=TRUE, pos=c('left', 'right')[i])
}
if(length(lev) == 2 && length(se)) {
xu <- sort(unique(x))
j1 <- gr == lev[1]
j2 <- gr == lev[2]
Y <- matrix(NA, nrow=length(xu), ncol=4,
dimnames=list(as.character(xu), c('y1', 'y2', 'se1', 'se2')))
x1 <- as.character(x)[j1]
x2 <- as.character(x)[j2]
Y[x1, 'y1' ] <- y [j1]
Y[x1, 'se1'] <- se[j1]
Y[x2, 'y2' ] <- y [j2]
Y[x2, 'se2'] <- se[j2]
ymid <- (Y[, 'y1'] + Y[, 'y2']) / 2.
halfwidthci <- qnorm(0.975) * sqrt(Y[, 'se1']^2 + Y[, 'se2']^2)
col <- adjustcolor('black', alpha.f=0.7)
grid::grid.segments(xu, ymid - 0.5 * halfwidthci,
xu, ymid + 0.5 * halfwidthci,
default.units='native',
gp=grid::gpar(col=col, lwd=1.5))
}
}
medvPanel <-
function(x, y, subscripts, groups=NULL, violin=TRUE, quantiles=FALSE,
...) {
sRequire('lattice')
gp <- length(groups)
plot.line <-
lattice::trellis.par.get(if(gp) "superpose.line" else "plot.line")
sym <- lattice::trellis.par.get(if(gp) "superpose.symbol" else "plot.symbol")
quant <- function(y) {
probs <- c(0.05, 0.125, 0.25, 0.375)
probs <- sort(c(probs, 1. - probs))
y <- y[! is.na(y)]
if(length(y) < 3) {
if(quantiles) {
w <- c(median(y), rep(NA, 9), length(y))
names(w) <- c('Median', format(probs), 'se', 'n')
}
else w <- c(Median=median(y), se=NA, n=length(y))
return(w)
}
w <- if(quantiles) hdquantile(y, probs)
m <- hdquantile(y, 0.5, se=TRUE)
se <- as.numeric(attr(m, 'se'))
c(Median=as.numeric(m), w, se=se, n=length(y))
}
denpoly <- function(x, y, col, n=50, pos, ...) {
y <- y[! is.na(y)]
n <- length(y)
if(n < 2) return()
den <- density(y, n=n, ...)
d <- den$y
y <- den$x
## Scale density of 0-3 mm
d <- 3 * d / max(d)
d <- c(d, d[length(d)])
mm <- grid::convertUnit(unit(d, 'mm'), 'mm', typeFrom='dimension')
kol <- if(n < 5 ) adjustcolor(col, alpha.f=0.2)
else if(n < 10) adjustcolor(col, alpha.f=0.4)
else col
grid::grid.polygon(y=unit(c(y, y[1]), 'native'),
x=if(pos == 'left') unit(x, 'native') - mm
else unit(x, 'native') + mm,
gp=grid::gpar(col=FALSE, fill=kol))
}
gr <- if(gp) groups[subscripts] else factor(rep('', length(x)))
lev <- levels(gr)
W <- NULL
for(i in 1 : length(lev)) {
j <- which(gr == levels(gr)[i])
xj <- x[j]
yj <- y[j]
w <- summarize(yj, xj, quant, type='matrix', keepcolnames=TRUE)
Y <- w$yj
xu <- w$xj
lattice::lpoints(xu, Y[,'Median'], cex=sym$cex[i], pch=sym$pch[i], col=sym$col[i],
alpha=sym$alpha[i])
lattice::llines(xu, Y[,'Median'], col=plot.line$col[i], lty=plot.line$lty[i],
lwd=plot.line$lwd[i], alpha=plot.line$alpha)
col <- plot.line$col[i]
if(violin) for(xx in sort(unique(xj)))
denpoly(xx, yj[xj == xx],
col=adjustcolor(plot.line$col[i], alpha.f=0.4),
pos=c('left', 'right')[i])
if(quantiles)
multLines(xu, Y[, colnames(Y) %nin% c('se', 'n'), drop=FALSE],
col=plot.line$col[i],
lty=plot.line$lty[i],
lwd=plot.line$lwd[i],
grid=TRUE, pos=c('left', 'right')[i])
W <- rbind(W, cbind(gr=levels(gr)[i], w))
}
if(length(lev) == 2) {
x <- W$xj
xu <- sort(unique(W$xj))
j1 <- W$gr == lev[1]
j2 <- W$gr == lev[2]
Y <- matrix(NA, nrow=length(xu), ncol=4,
dimnames=list(as.character(xu), c('y1', 'y2', 'se1', 'se2')))
x1 <- as.character(x)[j1]
x2 <- as.character(x)[j2]
Y[x1, 'y1' ] <- W$yj[j1, 'Median']
Y[x1, 'se1'] <- W$yj[j1, 'se']
Y[x2, 'y2' ] <- W$yj[j2, 'Median']
Y[x2, 'se2'] <- W$yj[j2, 'se']
ymid <- (Y[, 'y1'] + Y[, 'y2']) / 2.
halfwidthci <- qnorm(0.975) * sqrt(Y[, 'se1']^2 + Y[, 'se2']^2)
col <- adjustcolor('black', alpha.f=0.7)
grid::grid.segments(xu, ymid - 0.5 * halfwidthci,
xu, ymid + 0.5 * halfwidthci,
default.units='native',
gp=grid::gpar(col=col, lwd=1.5))
}
}
mbarclpl <- function(x, y, groups=NULL, yother, yvar=NULL,
maintracename='y', xlim=NULL, ylim=NULL,
xname='x', alphaSegments=0.45, ...) {
gp <- length(groups)
gr <- if(gp) groups else rep('', length(x))
gr <- as.factor(gr)
color <- if(gp) ~ .g.
lev <- levels(gr)
if(! length(yvar)) yvar <- rep('', length(x))
yvar <- as.factor(yvar)
se <- if('se' %in% colnames(yother)) yother[, 'se']
## prn(se, 'mbarclpl', fi='/tmp/z')
cy <- colnames(yother)
n <- if('n' %in% cy) yother[, 'n']
yother <- yother[, cy %nin% c('n', 'se'), drop=FALSE]
if(all(c('0.375', '0.625') %in% cy)) {
## If HD median estimate is not between 0.375 and 0.625 quantiles
## take it to be the closer of the two
y375 <- yother[, '0.375']
y625 <- yother[, '0.625']
y <- pmin(y, y625)
y <- pmax(y, y375)
}
fmt <- function(x) htmlSN(x, digits=5)
xdel <- 0.01 * diff(range(x, na.rm=TRUE))
xtxt <- paste0(xname, ': ', fmt(x), '<br>')
if(length(n)) xtxt <- paste0(xtxt, 'n:', n, '<br>')
R <- data.frame(x=x, y=y, yhi=NA, .g.=gr, .yvar.=yvar,
tracename = maintracename, connect=TRUE,
txt = paste0(xtxt, maintracename, ': ',
fmt(y)))
p <- ncol(yother)
half <- p / 2
x0 <- x
for(i in 1 : half) {
i1 <- i
i2 <- p - i + 1
tn <- paste0(colnames(yother)[i1], ' - ', colnames(yother)[i2])
x0 <- ifelse(gr == lev[1],
x - i * xdel, x + i * xdel)
txt <- paste0(xtxt, tn, ': [', fmt(yother[, i1]),
', ', fmt(yother[, i2]), ']')
r <- data.frame(x=x0, y=yother[, i1], yhi=yother[, i2],
.g.=gr, .yvar.=yvar,
tracename=tn, connect=NA, txt=txt)
R <- rbind(R, r)
}
if(length(lev) == 2 && length(se)) {
rr <- NULL
for(yv in levels(yvar)) {
k <- yvar == yv
xk <- x[k]; yk <- y[k]; sek <- se[k]; grk <- gr[k]
xu <- sort(unique(xk))
j1 <- grk == lev[1]
j2 <- grk == lev[2]
Y <- matrix(NA, nrow=length(xu), ncol=4,
dimnames=list(as.character(xu),
c('y1', 'y2', 'se1', 'se2')))
x1 <- as.character(xk)[j1]
x2 <- as.character(xk)[j2]
Y[x1, 'y1' ] <- yk [j1]
Y[x1, 'se1'] <- sek[j1]
Y[x2, 'y2' ] <- yk [j2]
Y[x2, 'se2'] <- sek[j2]
ymid <- (Y[, 'y1'] + Y[, 'y2']) / 2.
halfwidthci <- qnorm(0.975) * sqrt(Y[, 'se1']^2 + Y[, 'se2']^2)
ydel <- Y[, 'y2'] - Y[, 'y1']
txt <- paste0(xname, ': ', fmt(xu),
'<br>\u0394: ', fmt(ydel),
'<br>0.95 C.I. for \u0394: [',
fmt(ydel - halfwidthci), ', ',
fmt(ydel + halfwidthci), ']')
r <- data.frame(x = xu,
y = ymid - 0.5 * halfwidthci,
yhi = ymid + 0.5 * halfwidthci,
.g. = paste0(lev[1], ' vs. ', lev[2]),
.yvar. = yv, txt=txt,
tracename='\u00BD 0.95 C.I. for \u0394', connect=FALSE)
## prn(cbind(yother, 2 * rep(halfwidthci, each=2)), fi='/tmp/z')
rr <- rbind(rr, r)
}
R <- rbind(R, rr)
}
plotlyM(R, multplot=~.yvar., color=color, htext=~txt,
xlim=xlim, ylim=ylim, alphaSegments=alphaSegments, ...)
}
medvpl <-
function(x, y, groups=NULL, yvar=NULL, maintracename='y',
xlim=NULL, ylim=NULL, xlab=xname, ylab=NULL, xname='x',
zeroline=FALSE, yother=NULL, alphaSegments=0.45,
dhistboxp.opts=NULL, ...) {
gp <- length(groups)
gr <- if(gp) groups else factor(rep('', length(x)))
lev <- levels(gr)
yvarpres <- length(yvar)
if(! length(yvar)) yvar <- rep('', length(x))
if(! length(ylab)) ylab <- structure(unique(yvar), names=yvar)
if(! length(names(ylab))) stop('ylab must have names')
fmt <- function(x) htmlSN(x, digits=5)
R <- NULL
for(yv in unique(yvar)) {
k <- yvar == yv
xk <- x[k]
yk <- y[k]
gk <- gr[k]
r <- do.call('dhistboxp',
c(list(yk, group=gk, strata=xk, xlab=ylab[yv]),
dhistboxp.opts))
r$strata <- NULL
ry <- r$y
r$y <- r$x
r$x <- ry
ryhi <- r$yhi
r$yhi <- r$xhi
r$xhi <- ryhi
r$yvar <- yv
R <- rbind(R, r)
xku <- unique(xk)
if(length(lev) == 2) {
for(xa in xku) {
Y <- matrix(NA, nrow=length(xku), ncol=4,
dimnames=list(as.character(xku), c('y1','y2','se1','se2')))
for(ga in lev) {
j <- xk == xa & gk == ga
ykj <- yk[j]
ykj <- ykj[! is.na(ykj)]
if(length(ykj) < 3) {
m <- median(ykj)
se <- NA
}
else {
m <- hdquantile(ykj, probs=0.5, se=TRUE)
se <- attr(m, 'se')
}
med <- as.vector(m)
if(ga == lev[1]) {
med1 <- med
se1 <- se
} else {
med2 <- med
se2 <- se
}
}
ydel <- med2 - med1
ymid <- (med1 + med2) / 2.
halfwidthci <- qnorm(0.975) * sqrt(se1 ^ 2 + se2 ^ 2)
txt <- paste0(xname, ': ', fmt(xa),
'<br>\u0394: ', fmt(ydel),
'<br>0.95 C.I. for \u0394: [',
fmt(ydel - halfwidthci), ', ',
fmt(ydel + halfwidthci), ']')
R <- rbind(R,
data.frame(x =xa, y =ymid - 0.5 * halfwidthci,
xhi=NA, yhi=ymid + 0.5 * halfwidthci,
group=paste0(lev[1], ' vs. ', lev[2]),
yvar = yv,
txt = txt,
type = '', connect=NA))
}
}
}
R$group <- paste(R$group, R$type)
R$type <- NULL
plotlyM(R, htext=~txt, multplot=~yvar, color=~group,
alphaSegments=alphaSegments,
xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim, zeroline=zeroline,
...)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.