histSpikeg <- function(formula=NULL, predictions=NULL, data, plotly=NULL,
lowess=FALSE, xlim=NULL, ylim=NULL,
side=1, nint=100,
frac=function(f) 0.01 + 0.02*sqrt(f-1) /
sqrt(max(f,2)-1),
span=3/4, histcol='black', showlegend=TRUE) {
## Raw data in data, predicted curves are in predictions
## If predictions is not given, side (1 or 3) is used
v <- all.vars(formula)
yv <- v[ 1]
xv <- v[-1]
X <- xv[1]
if(lowess) {
if(yv %nin% names(data))
stop(paste(yv, 'must be in data if lowess=TRUE'))
yval <- data[[yv]]
iter <- if(length(unique(yval[! is.na(yval)])) > 2) 3 else 0
lows <- function(x, y, iter, span) {
i <- ! is.na(x + y)
lowess(x[i], y[i], iter=iter, f=span)
}
}
if(length(xv) > 1 && ! (lowess || length(predictions)))
stop('predictions must be given or lowess=TRUE if formula has > 1 variable on the right')
x1 <- data[[X]]
if(length(xlim)) {
data <- data[x1 >= xlim[1] & x1 <= xlim[2], ]
x1 <- data[[X]]
}
if(length(unique(x1)) > nint) {
eps <- diff(range(x1, na.rm=TRUE)) / nint
x1 <- round(x1 / eps) * eps
data[[X]] <- x1
}
p <- predictions
# xr <- NULL
# if(length(p)) xr <- range(p[[X]], na.rm=TRUE)
# else if(lowess) xr <- range(data[[X]], na.rm=TRUE)
xr <- xlim
if(length(xv) == 1) {
tab <- as.data.frame(do.call(table, data[X]))
tab[[X]] <- as.numeric(as.character(tab[[X]]))
if(length(xr)) tab <- tab[tab[[X]] >= xr[1] & tab[[X]] <= xr[2], ]
if(lowess) {
p <- as.data.frame(lows(data[[X]], yval, iter, span))
names(p) <- c(X, yv)
}
if(length(p)) tab$.yy. <- approxExtrap(p[[X]], p[[yv]], xout=tab[[X]])$y
} else { ## grouping variable(s) present
tab <- as.data.frame(do.call(table, data[xv]))
tab <- subset(tab, Freq > 0)
tab[[X]] <- as.numeric(as.character(tab[[X]]))
if(length(xr)) tab <- tab[tab[[X]] >= xr[1] & tab[[X]] <= xr[2], ]
tab$.yy. <- rep(NA, nrow(tab))
gv <- xv[-1] # grouping variables
U <- unique(tab[gv]) # unique combinations
if(! lowess) {
for(k in 1 : nrow(U)) {
i <- rep(TRUE, nrow(tab))
j <- rep(TRUE, nrow(p))
for(l in 1 : ncol(U)) {
currgroup <- as.character(U[k, gv[l], drop=TRUE])
i <- i & (tab[[gv[l]]] == currgroup)
j <- j & (p[[gv[l]]] == currgroup)
} ## now all grouping variables intersected
tab$.yy.[i] <- approxExtrap(p[[X]][j], p[[yv]][j], xout=tab[[X]][i])$y
}
} ## end if(! lowess)
else { # lowess; need to compute p
p <- NULL
for(k in 1 : nrow(U)) {
i <- rep(TRUE, nrow(tab))
j <- rep(TRUE, nrow(data))
for(l in 1 : ncol(U)) {
currgroup <- as.character(U[k, gv[l], drop=TRUE])
i <- i & (tab[[gv[l]]] == currgroup)
j <- j & (data[[gv[l]]] == currgroup)
} ## now all grouping variables intersected
sm <- lows(data[[X]][j], data[[yv]][j], iter, span)
Sm <- sm; names(Sm) <- c(X, yv)
Uk <- U[k, , drop=FALSE]; row.names(Uk) <- NULL
p <- rbind(p, data.frame(Uk, Sm))
tab$.yy.[i] <- approxExtrap(sm, xout=tab[[X]][i])$y
}
} ## end lowess
} ## end grouping variables present
if(! length(ylim)) {
if(length(data) && yv %in% names(data))
ylim <- range(data[[yv]], na.rm=TRUE)
else if(length(p)) ylim <- range(p[[yv]], na.rm=TRUE)
}
if(! length(ylim)) stop('no way to infer ylim from information provided')
tab$.rf. <- frac(tab$Freq) * diff(ylim)
n <- nrow(tab)
tab$.ylo. <- if(length(p)) tab$.yy. - tab$.rf.
else if(side == 1) rep(ylim[1], n) else rep(ylim[2], n)
tab$.yhi. <- if(length(p)) tab$.yy. + tab$.rf.
else if(side == 1) ylim[1] + tab$.rf. else ylim[2] - tab$.rf.
hcol <- if(histcol == 'default') '' else sprintf(', col="%s"', histcol)
P <- plotly
if(! is.null(P) && requireNamespace("plotly")) {
tab$.xx. <- tab[[X]]
## Got NA/NaN argument, numerical expression has 8 elements
## when using add_segments
# P <- plotly:add_segments(P, data=tab,
# x=~ .xx., xend=~ .xx.,
# y=~ .ylo., yend=~ .yhi.,
# hoverinfo='none',
# name='Histogram', legendgroup='Histogram')
n <- 3 * nrow(tab)
.xx. <- .yy. <- rep(NA, n)
.xx.[seq(1, n, by=3)] <- tab[[X]]
.xx.[seq(2, n, by=3)] <- tab[[X]]
.yy.[seq(1, n, by=3)] <- tab$.ylo.
.yy.[seq(2, n, by=3)] <- tab$.yhi.
P <- plotly::add_lines(P, data=data.frame(.xx., .yy.),
x=~.xx., y=~.yy.,
line=list(color=histcol, width=1.4),
hoverinfo='none', showlegend=showlegend,
name='Histogram', legendgroup='Histogram')
if(lowess) {
af <- function(x) as.formula(paste('~', x))
P <- plotly::add_lines(P, data=p, x=af(X), y=af(yv),
hoverinfo='none', showlegend=showlegend,
name='loess', legendgroup='loess')
}
return(P)
}
res <- eval(parse(text=sprintf('ggplot2::geom_segment(data=tab,
aes(x=%s, xend=%s,
y=.ylo., yend=.yhi.), size=.25 %s)', X, X, hcol)))
if(lowess) res <- list(hist=res, lowess=eval(parse(text=
sprintf('ggplot2::geom_line(data=p, aes(x=%s, y=%s))', X, yv))))
res
}
utils::globalVariables(c('aes', 'Freq', '.ylo.', '.yhi.', 'x', 'y'))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.