# R convenience functions
#
# Author: mg14
###############################################################################
#' Function to transform P-values into symbols of significance (***)
#' @param s A vector of P-values
#' @param breaks The breaks
#' @param labels The labels
#' @return A character vector.
#'
#' @author mg14
#' @export
sig2star = function(s, breaks=c(-Inf, 0.001,0.01,0.05,0.1,1), labels=c("***","**","*",".","")) {
r <- s
r[] <- as.character(cut(s, breaks=breaks, labels))
r[is.na(r)] <- ""
r
}
#' Plot a donut
donut <- function(x, r0=0.5, r1=.75,...) {
pie(x, border=NA, labels = paste(names(x), " (",round(100*x),"%)", sep=""), radius=r1, ...)
polygon(cos(seq(0,2*pi,l=100))*r0, sin(seq(0,2*pi,l=100))*r0, col="white", border=NA)
}
#' Compute entropy (base e) from a vector of probabilities
#'
#' Given p, the entropy is computed as -sum(log(p^p)).
#' @param p Vector of probabilities
#' @return numeric
#'
#' @author mg14
#' @export
entropy <- function(p) {
stopifnot(abs(sum(p)-1) < length(p)*.Machine$double.eps)
-sum(log(p^p))
}
require(RColorBrewer)
set1 <- brewer.pal(9,"Set1")
plotcvnet = function(cvnet,X, main="", simple.annot=NULL, col=set1, col0="grey", cex=sqrt(abs(colMeans(X, na.rm=TRUE))[n]+0.02)*5,xlim=c(0.5,ncol(X)),...){
first=function(x) which(x)[1]
lollyplot <- function(cvnet, beta, ..., X, main, cex) {
r <- rank(cvnet$nzero[apply(beta!=0,1,first)], ties.method="first")
R <- 1 - cvnet$cvm/cvnet$cvm[1]
Rup <- 1 - cvnet$cvlo/cvnet$cvm[1]
Rlo <- 1 - cvnet$cvup/cvnet$cvm[1]
w <- apply(beta!=0,1,first)
beta <- beta[,which.min(cvnet$cvm)]
#n = names(beta)
#ix <- which(c(diff(cvnet$nzero)>0, TRUE) & cvnet$nzero > 0)
#xx <- c(0,cvnet$nzero[ix]) + 0.5
ix <- na.omit(w)
xx <- rank( ix, ties.method="first") + 0.5
o <- order(xx)
ix <- ix[o]
ix <- c(ix[1],ix)
xx <- c(0.5,xx[o])
#ix = c(ix[1],ix)
plot(xx, R[ix], type="S", xaxt="n", xlab="", ylim=c(0,max(Rup)), ylab=expression(paste("Explained variance ",R^2)), main = main, xlim=xlim, col=set1[1],...)
polygon(rep(c(xx,rev(xx)), each=2)[-2*length(ix) + c(0,-1)], rep(c(Rup[ix],rev(Rlo[ix])), each=2)[-c(1,4*length(ix))], border=NA, col = paste(set1[1],"22", sep=""))
lines(xx, cvnet$lambda[ix]/max(cvnet$lambda)*max(R), lty=1, type="S", col="black")
u = par("usr")[3:4]
#mtext(colnames(oncocyto), at=r, side=1, las=2, font=ifelse(grepl("[A-Z]", colnames(oncocyto)),3,1), cex=.7, line=0.5)
n <- colnames(X)#names(ix)[-1]
colLab <- ifelse(n %in% names(simple.annot), col[simple.annot[n]], "black")
colLab[w > which.max(R) | is.na(w)] <- "grey"
rotatedLabel(r, rep(0,length(r)), colnames(X), font=ifelse(grepl("^[A-Z]", n) &! grepl("PC", n),3,1), cex=0.66, col=colLab)
scale <- diff(range(beta, na.rm=TRUE))*1.05
if(scale ==0) scale=1
shift <- mean(range(beta, na.rm=TRUE))#median(c, na.rm=TRUE)
#lines(r,(beta - shift)/scale*diff(u)*.7 + mean(u), col="grey", type="h", lty=1)
abline(h=(0 - shift)/scale*diff(u)*.7 + mean(u), col="grey", lty=2)
points(r,(beta - shift)/scale*diff(u)*.7 + mean(u), col=colLab, cex=cex, pch=ifelse(beta==0,1,19))
axis(4, labels=pretty(beta),at = (pretty(beta) - shift)/scale*diff(u)*.7 + mean(u), col="darkgrey", col.axis="darkgrey")
mtext(side=4, line=3, "Coefficient", col="darkgrey", at=mean(u), las=0)
}
if(class(cvnet$glmnet.fit)[1] != "multnet"){
#r <- rank(cvnet$nzero[apply(cvnet$glmnet.fit$beta!=0,1,first)], ties.method="first")
#w = which.min(cvnet$cvm)
#beta = cvnet$glmnet.fit$beta[,w]
#c = c[c!=0]
lollyplot(cvnet = cvnet, beta=cvnet$glmnet.fit$beta, ... = ..., X = X, main=main, cex=cex)
}else {j = 1; for(beta in cvnet$glmnet.fit$beta){
#w = which.min(cvnet$cvm)
#beta = beta[,w]
#c = c[c!=0]
#r <- rank(cvnet$nzero[apply(beta!=0,1,first)], ties.method="first")
lollyplot(cvnet = cvnet, beta=beta, ... = ..., X = X, main=paste(main, ": ", names(cvnet$glmnet.fit$beta)[j], sep=""), cex=cex)
j = j+1
}
}
}
#' Plot a rotated axis label
#' @param x0 The x position
#' @param y0 The y position
#' @param labels The labels
#' @param pos The positioning relative to the coordinate. Default = 1 (below).
#' @param cex The character expansion factor. Default = 1.
#' @param srt The degree of rotation. Default = 45.
#' @param ...
#' @return NULL
#'
#' @author mg14
#' @export
rotatedLabel <- function(x0 = seq_along(labels), y0 = rep(par("usr")[3], length(labels)), labels, pos = 1, cex=1, srt=45, ...) {
w <- strwidth(labels, units="user", cex=cex)
h <- strheight(labels, units="user",cex=cex)
u <- par('usr')
p <- par('plt')
f <- par("fin")
xpd <- par("xpd")
par(xpd=NA)
text(x=x0 + ifelse(pos==1, -1,1) * w/2*cos(srt/360*2*base::pi), y = y0 + ifelse(pos==1, -1,1) * w/2 *sin(srt/360*2*base::pi) * (u[4]-u[3])/(u[2]-u[1]) / (p[4]-p[3]) * (p[2]-p[1])* f[1]/f[2] , labels, las=2, cex=cex, pos=pos, adj=1, srt=srt,...)
par(xpd=xpd)
}
myscale <- function(data){
apply(data,
2,
function(x){
if(all(x %in% c(0,1,NA))) # Don't scale categorical variables, only re-center
#x
scale(x, scale=FALSE)
else
scale(x)
#rank(x) / sum(!is.na(x)) - 0.5
}
)
}
#' Jitter by violins
#' @param x A vector of numbers
#' @param magnitude The overall magnitude of scatter
#' @return NULL
#'
#' @author mg14
#' @export
violinJitter <- function(x, magnitude=1){
d <- density(x)
data.frame(x=x, y=runif(length(x),-magnitude/2, magnitude/2) * approxfun(d$x, d$y)(x))
}
#' Violins plot with jittered data points
#' @param y the values to be plotted
#' @param x the levels to be jittered
#' @param col the colors to be used for the outline
#' @param col.pty the colors used for the points
#' @param magnitude the magnitude of the violins
#' @param ... additional parameters passed to plot()
#' @return NULL
#'
#' @author mg14
#' @examples
#' f <- rep(1:5, each=50)
#' violinJitterPlot(rnorm(length(f), f), factor(f))
#' @export
violinJitterPlot <- function(y, x=factor(rep(1, length(y))), col=1:nlevels(x), col.pty = colTrans(col), magnitude=1,xlab="",ylab="", plot.violins=TRUE,...){
yl <- split(y,x)
o <- order(x)
yj <- do.call("rbind",lapply(yl, violinJitter, magnitude = magnitude))
dt <- lapply(yl, density)
s <- max(sapply(dt, function(x) max(x$y)))
plot(na.omit(as.numeric(x)[o]) + yj$y/s, yj$x, col=col.pty[na.omit(as.numeric(x)[o])], xlab=xlab, ylab=ylab, xaxt="n",...)
axis(side=1, at=1:nlevels(x), labels=levels(x))
qt <- lapply(yl, quantile)
if(plot.violins)
for(i in 1:length(dt)){
lines(i+dt[[i]]$y*magnitude/2/s, dt[[i]]$x, col=col[i])
lines(i-dt[[i]]$y*magnitude/2/s, dt[[i]]$x, col=col[i])
for(q in 2:4){
w <- which.min(abs(dt[[i]]$x-qt[[i]][q]))
segments(i+dt[[i]]$y[w]*magnitude/2/s, dt[[i]]$x[w],i-dt[[i]]$y[w]*magnitude/2/s, dt[[i]]$x[w], lwd=ifelse(q==3,2,1))
}
}
}
#' Stars plot
#'
#' Modified from graphics::stars()
#' @param x data.frame or matrix to plot
#' @return NULL
#'
#' @author mg14
stars <- function (x, full = TRUE, scale = TRUE, radius = TRUE, labels = dimnames(x)[[1L]],
locations = NULL, nrow = NULL, ncol = NULL, len = 1, key.loc = NULL,
key.labels = dimnames(x)[[2L]], key.xpd = TRUE, xlim = NULL,
ylim = NULL, flip.labels = NULL, draw.segments = FALSE, col.segments = 1L:n.seg,
col.stars = NA, col.lines = NA, axes = FALSE, frame.plot = axes,
main = NULL, sub = NULL, xlab = "", ylab = "", cex = 0.8,
lwd = 0.25, lty = par("lty"), xpd = FALSE, mar = pmin(par("mar"),
1.1 + c(2 * axes + (xlab != ""), 2 * axes + (ylab != ""), 1, 0)), add = FALSE, plot = TRUE, density=NA, init.angle=0, ...)
{
if (is.data.frame(x))
x <- data.matrix(x)
else if (!is.matrix(x))
stop("'x' must be a matrix or a data frame")
if (!is.numeric(x))
stop("data in 'x' must be numeric")
n.loc <- nrow(x)
n.seg <- ncol(x)
if (is.null(locations)) {
if (is.null(nrow))
nrow <- ceiling(if (!is.numeric(ncol)) sqrt(n.loc) else n.loc/ncol)
if (is.null(ncol))
ncol <- ceiling(n.loc/nrow)
if (nrow * ncol < n.loc)
stop("'nrow * ncol' is less than the number of observations")
ff <- if (!is.null(labels))
2.3
else 2.1
locations <- expand.grid(ff * 1L:ncol, ff * nrow:1)[1L:n.loc,
]
if (!is.null(labels) && (missing(flip.labels) || !is.logical(flip.labels)))
flip.labels <- ncol * mean(nchar(labels, type = "c")) >
30
}
else {
if (is.numeric(locations) && length(locations) == 2) {
locations <- cbind(rep.int(locations[1L], n.loc),
rep.int(locations[2L], n.loc))
if (!missing(labels) && n.loc > 1)
warning("labels do not make sense for a single location")
else labels <- NULL
}
else {
if (is.data.frame(locations))
locations <- data.matrix(locations)
if (!is.matrix(locations) || ncol(locations) != 2)
stop("'locations' must be a 2-column matrix.")
if (n.loc != nrow(locations))
stop("number of rows of 'locations' and 'x' must be equal.")
}
if (missing(flip.labels) || !is.logical(flip.labels))
flip.labels <- FALSE
}
xloc <- locations[, 1]
yloc <- locations[, 2]
angles <- if (full)
seq.int(0, 2 * pi, length.out = n.seg + 1)[-(n.seg +
1)] + init.angle
else if (draw.segments)
seq.int(0, pi, length.out = n.seg + 1)[-(n.seg + 1)] + init.angle
else seq.int(0, pi, length.out = n.seg)
if (length(angles) != n.seg)
stop("length of 'angles' must equal 'ncol(x)'")
if (scale) {
x <- apply(x, 2L, function(x) (x - min(x, na.rm = TRUE))/diff(range(x,
na.rm = TRUE)))
}
x[is.na(x)] <- 0
mx <- max(x <- x * len)
if (is.null(xlim))
xlim <- range(xloc) + c(-mx, mx)
if (is.null(ylim))
ylim <- range(yloc) + c(-mx, mx)
deg <- pi/180
op <- par(mar = mar, xpd = xpd)
on.exit(par(op))
dev.hold()
on.exit(dev.flush(), add = TRUE)
if (plot && !add)
plot(0, type = "n", ..., xlim = xlim, ylim = ylim, main = main,
sub = sub, xlab = xlab, ylab = ylab, asp = 1, axes = axes)
if (!plot)
return(locations)
s.x <- xloc + x * rep.int(cos(angles), rep.int(n.loc, n.seg))
s.y <- yloc + x * rep.int(sin(angles), rep.int(n.loc, n.seg))
if (draw.segments) {
aangl <- c(angles, if (full) 2 * pi + init.angle else pi+init.angle)
for (i in 1L:n.loc) {
px <- py <- numeric()
for (j in 1L:n.seg) {
k <- seq.int(from = aangl[j], to = aangl[j +
1], by = 1 * deg)
px <- c(px, xloc[i], s.x[i, j], x[i, j] * cos(k) +
xloc[i], NA)
py <- c(py, yloc[i], s.y[i, j], x[i, j] * sin(k) +
yloc[i], NA)
}
polygon(px, py, col = col.segments, lwd = lwd, lty = lty, border=col.lines, density=density)
}
}
else {
for (i in 1L:n.loc) {
polygon(s.x[i, ], s.y[i, ], lwd = lwd, lty = lty,
border = col.lines[i], col = col.stars[i], density=density[i])
if (radius)
segments(rep.int(xloc[i], n.seg), rep.int(yloc[i],
n.seg), s.x[i, ], s.y[i, ], lwd = lwd, lty = lty)
}
}
if (!is.null(labels)) {
y.off <- mx * (if (full)
1
else 0.1)
if (flip.labels)
y.off <- y.off + cex * par("cxy")[2L] * ((1L:n.loc)%%2 -
if (full)
0.4
else 0)
text(xloc, yloc - y.off, labels, cex = cex, adj = c(0.5,
1))
}
if (!is.null(key.loc)) {
par(xpd = key.xpd)
key.x <- len * cos(angles) + key.loc[1L]
key.y <- len * sin(angles) + key.loc[2L]
if (draw.segments) {
px <- py <- numeric()
for (j in 1L:n.seg) {
k <- seq.int(from = aangl[j], to = aangl[j +
1], by = 1 * deg)
px <- c(px, key.loc[1L], key.x[j], len * cos(k) +
key.loc[1L], NA)
py <- c(py, key.loc[2L], key.y[j], len * sin(k) +
key.loc[2L], NA)
}
polygon(px, py, col = col.segments, lwd = lwd, lty = lty, border=col.lines, density=density)
}
else {
polygon(key.x, key.y, lwd = lwd, lty = lty)
if (radius)
segments(rep.int(key.loc[1L], n.seg), rep.int(key.loc[2L],
n.seg), key.x, key.y, lwd = lwd, lty = lty)
}
lab.angl <- angles + if (draw.segments)
(angles[2L] - angles[1L])/2
else 0
label.x <- 1.1 * len * cos(lab.angl) + key.loc[1L]
label.y <- 1.1 * len * sin(lab.angl) + key.loc[2L]
for (k in 1L:n.seg) {
text.adj <- c(if (lab.angl[k] < 90 * deg || lab.angl[k] >
270 * deg) 0 else if (lab.angl[k] > 90 * deg &&
lab.angl[k] < 270 * deg) 1 else 0.5, if (lab.angl[k] <=
90 * deg) (1 - lab.angl[k]/(90 * deg))/2 else if (lab.angl[k] <=
270 * deg) (lab.angl[k] - 90 * deg)/(180 * deg) else 1 -
(lab.angl[k] - 270 * deg)/(180 * deg))
text(label.x[k], label.y[k], labels = key.labels[k],
cex = cex, adj = text.adj)
}
}
if (frame.plot)
box(...)
invisible(locations)
}
ggPlot <- function(x, g, width=.8, xlab="", ylab="",...){
g <- factor(g)
s <- split(x, g)
o <- order(sapply(s[levels(g)], function(x) if(!all(is.na(x))) median(x, na.rm=TRUE) else NA), na.last=TRUE)
g <- factor(g, levels=levels(g)[o])
s <- split(x,g)
plot(unlist(sapply(s[levels(g)], function(x) (rank(x)-1)/(length(x)-1) - .5)) * width + sort(as.numeric(g)), unlist(s), xaxt="n", xlab=xlab, ylab=ylab, ...)
m <- sapply(s, median, na.rm=TRUE)
segments(1:nlevels(g)-width/2, m, 1:nlevels(g)+width/2, m)
abline(v=1:nlevels(g), lty=3)
axis(side=1, at=1:nlevels(g), labels = levels(g))
}
MoscowPlot <- function(x, sd, col=1:length(x), at = seq_along(x), width = 0.9, ...){
y <- seq(0,max(x + 3*sd), l=100)
plot(NA,NA, xlim=c(0,length(x)+1), ylim=range(y), xaxt="n", ...)
for(i in seq_along(x)){
xx <- width/2 * pnorm(y, x[i], sd[i], lower.tail = FALSE)
polygon(i +c(-xx, rev(xx)), c(y, rev(y)) , col=col[i], border=NA)
}
}
lasso <- function(..., theta=1, scale=TRUE) {
x <- cbind(...)
nvar <- ncol(x)
xname <- as.character(parse(text=substitute(cbind(...))))[-1]
vars <- apply(x, 2, function(z) var(z[!is.na(z)]))
class(x) <- 'coxph.penalty'
if (missing(theta))
stop("Missing theta")
if (scale)
pfun <- function(coef,theta, ndead, scale) {
list(penalty= sum(abs(coef) *scale)*theta,
first = theta*sign(coef)*scale,
second = 0,
flag=FALSE)
}
else
pfun <- function(coef,theta, ndead, scale) {
list(penalty= sum(abs(coef))*theta,
first = theta*sign(coef),
second = 0,
flag=FALSE)
}
temp <- list(pfun=pfun,
diag=TRUE,
cfun=function(parms, iter, history) {
list(theta=parms$theta, done=TRUE) },
cparm=list(theta= theta),
pparm= vars,
varname=paste('lasso(', xname, ')', sep=''))
attributes(x) <- c(attributes(x), temp)
x
}
#' A pseudo transparency to a color
#' @param col A color
#' @param f The amount of pseudotransparency
#' @return color
#'
#' @author mg14
#' @export
colTrans <- function(col, f=2){
hsv <- apply(col2rgb(col), 2, rgb2hsv)
tmp <- car::logit(hsv[2,]*.95 + 0.025) - f
hsv[2,] <- pmin(1,pmax(0,(1/(exp(-tmp) +1) - 0.025)/.95))
tmp <- car::logit(hsv[3,]*.95 + 0.025) + f
hsv[3,] <- pmin(1,pmax(0,(1/(exp(-tmp) +1) - 0.025)/.95))
apply(hsv,2, function(x) hsv(x[1],x[2],x[3]))
}
jitterBox <- function(x,y, xlim = c(1,nlevels(x)) + c(-.5,.5),...){
#boxplot(y ~ x, ..., notch=TRUE)
plot(jitter(as.numeric(x)),y, xlim = xlim, ..., xaxt="n", pch=16)
s <- split(y,x)
m <- sapply(s, mean, na.rm=TRUE)
n <- 1:length(m)
l <- sapply(s, length)
c <- qnorm(0.975,0,sapply(s, sd, na.rm=TRUE)/sqrt(l))
rect(n-.2,m-c,n+.2,m+c,border=NA, col="#88888844")
segments(n-.2,m,n+.2,m, lwd=2)
axis(side = 1, at=n, labels=levels(x))
if(length(unique(x))>1){
a <- anova(lm(y~x))
mtext(side=3, line=0, text= paste("P =", format(a$`Pr(>F)`[1],digits=2)))
}
}
survdiff <- function (formula, data, subset, na.action, rho = 0, continuity = FALSE)
{
call <- match.call()
m <- match.call(expand.dots = FALSE)
m$rho <- NULL
m$continuity <- NULL
if (!inherits(formula, "formula"))
stop("The 'formula' argument is not a formula")
Terms <- if (missing(data))
terms(formula, "strata")
else terms(formula, "strata", data = data)
m$formula <- Terms
m[[1]] <- as.name("model.frame")
if (is.R())
m <- eval(m, parent.frame())
else m <- eval(m, sys.parent())
y <- model.extract(m, "response")
if (!inherits(y, "Surv"))
stop("Response must be a survival object")
if (attr(y, "type") != "right")
stop("Right censored data only")
ny <- ncol(y)
n <- nrow(y)
offset <- attr(Terms, "offset")
if (!is.null(offset)) {
offset <- as.numeric(m[[offset]])
if (length(attr(Terms, "factors")) > 0)
stop("Cannot have both an offset and groups")
if (any(offset < 0 | offset > 1))
stop("The offset must be a survival probability")
expected <- sum(-log(offset))
observed <- sum(y[, ny])
if (rho != 0) {
num <- sum(1/rho - ((1/rho + y[, ny]) * offset^rho))
var <- sum(1 - offset^(2 * rho))/(2 * rho)
}
else {
var <- sum(-log(offset))
num <- var - observed + if(continuity) 1/2 else 0
}
chi <- num * num/var
rval <- list(n = n, obs = observed, exp = expected, var = var,
chisq = chi)
}
else {
strats <- attr(Terms, "specials")$strata
if (length(strats)) {
temp <- untangle.specials(Terms, "strata", 1)
dropx <- temp$terms
if (length(temp$vars) == 1)
strata.keep <- m[[temp$vars]]
else strata.keep <- strata(m[, temp$vars], shortlabel = TRUE)
}
else strata.keep <- rep(1, nrow(m))
if (length(strats))
ll <- attr(Terms[-dropx], "term.labels")
else ll <- attr(Terms, "term.labels")
if (length(ll) == 0)
stop("No groups to test")
else groups <- strata(m[ll])
fit <- survdiff.fit(y, groups, strata.keep, rho)
if (is.matrix(fit$observed)) {
otmp <- apply(fit$observed, 1, sum)
etmp <- apply(fit$expected, 1, sum)
}
else {
otmp <- fit$observed
etmp <- fit$expected
}
df <- (etmp > 0)
if (sum(df) < 2)
chi <- 0
else {
temp2 <- ((otmp - etmp - if(continuity) 1/2 else 0)[df])[-1]
vv <- (fit$var[df, df])[-1, -1, drop = FALSE]
chi <- sum(solve(vv, temp2) * temp2)
}
rval <- list(n = table(groups), obs = fit$observed, exp = fit$expected,
var = fit$var, chisq = chi)
if (length(strats))
rval$strata <- table(strata.keep)
}
na.action <- attr(m, "na.action")
if (length(na.action))
rval$na.action <- na.action
rval$call <- call
if (is.R())
class(rval) <- "survdiff"
else oldClass(rval) <- "survdiff"
rval
}
#' Merge levels of a factor
#' @param f A factor
#' @param mergeList A list of type list(newlevel1 = c(oldlevel1, oldlevel2), newlevel2=c(...),...)
#' @return A factor
#'
#' @author mg14
#' @export
mergeLevels <- function(f, mergeList){
oldLevels <- setdiff(levels(f), unlist(mergeList))
newLevels <- as.list(oldLevels)
names(newLevels) <- oldLevels
newLevels <- c(newLevels, mergeList)
newFactor <- f
levels(newFactor) <- newLevels
return(newFactor)
}
#' Replace NAs with zeros
#' @param X
#' @return same as X
#'
#' @author mg14
#' @export
na.zero <- function(X){
w <- is.na(X)
X[w] <- 0
attr(X, "na.action") <- which(w)
return(X)
}
#' Replace NAs with mean
#' @param X
#' @return same as X
#'
#' @author mg14
#' @export
na.mean <- function(X){
w <- is.na(X)
X[w] <- mean(X, na.rm=TRUE)
attr(X, "na.action") <- which(w)
return(X)
}
glmnet.formula <- function (formula, data, subset, ...)
{
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- quote(stats::model.frame)
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "any")
if (is.empty.model(mt)) {
stop()
}
x <- model.matrix(mt, mf, contrasts)
g <- glmnet(x[,-1],y,...)
g$y <- y
g$x <- x
return(g)
}
plotSurv <- function(f, col=1:(length(terms(f))-1)^2 , ...){
s <- survfit(f, ...)
c <- coxph(f, ...)
summary(c)
p <- 2*pnorm(abs(diff(c$coefficients)),sd=sqrt(sum(diag(c$var))), lower.tail=FALSE)
plot(s, col=col, mark=NA)
legend("topright", bty="", rownames(summary(s)$table), col=col, lty=1)
title(main=paste("P =",signif(p,2)), font=1)
}
splitSurv <- function(surv, timeTpl){
timeSurv <- surv[,1]
w <- which(timeTpl < timeSurv)
t1 <- c(rep(0,nrow(surv)), timeTpl[w])
t2 <- c(pmin(timeSurv, timeTpl, na.rm=TRUE), timeSurv[w])
event <- c(surv[,2], surv[w,2])
event[w] <- 0
Surv(t1,t2,event)
}
splitIndex <- function(surv, timeTpl){
timeSurv <- surv[,1]
w <- which(timeTpl < timeSurv)
c(1:nrow(surv), w)
}
makeTimeDependent <- function(dataFrame, timeTpl, timeSurv = dataFrame$time, time0Surv = rep(0, nrow(dataFrame)), event=dataFrame$event){
w <- which(timeTpl < timeSurv & timeTpl > time0Surv)
index <- c(1:nrow(dataFrame), w)
d <- dataFrame[index,]
d$index <- index
d$time1 <- c(time0Surv, timeTpl[w])
d$time2 <- c(pmin(timeSurv, timeTpl, na.rm=TRUE), timeSurv[w])
d$transplant <- c(rep(0,nrow(dataFrame)), rep(1, length(w)))
e <- c(event, event[w])
e[w] <- 0
d$event <- e
return(d)
}
pbound <- function(x, lower, upper){
pmin(pmax(x,lower),upper)
}
humanSvg <- '<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by grConvert v0.1-0 -->
<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" width="35pt" height="84pt" viewBox="0 0 35 84" version="1.1">
<g id="surface1">
<path style="fill-rule:nonzero;fill:rgb(100%,100%,100%);fill-opacity:1;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke:rgb(0%,0%,0%);stroke-opacity:1;stroke-miterlimit:10;" d="M 34.099163 47.801325 C 34.099163 50.101674 33.001535 51.199302 30.899693 51.199302 C 28.801744 51.199302 27.700223 50.101674 27.700223 47.801325 L 27.700223 25.000237 L 26.100488 25.000237 L 26.100488 79.901116 C 26.100488 82.201465 24.80046 83.299093 22.200405 83.299093 C 19.600349 83.299093 18.300321 82.201465 18.300321 79.901116 L 18.300321 47.400419 L 16.299679 47.400419 L 16.299679 79.901116 C 16.299679 82.201465 14.999651 83.299093 12.399595 83.299093 C 9.79954 83.299093 8.499512 82.201465 8.499512 79.901116 L 8.499512 25.000237 L 6.899777 25.000237 L 6.899777 47.801325 C 6.899777 50.101674 5.798256 51.199302 3.700307 51.199302 C 1.598465 51.199302 0.500837 50.000474 0.500837 47.801325 L 0.500837 21.699568 C 0.500837 17.499777 3.4006 15.300628 9.200126 15.300628 L 25.399874 15.300628 C 31.1994 15.300628 34.099163 17.398577 34.099163 21.699568 Z M 24.099847 7.298061 C 24.099847 9.201395 23.399233 10.898438 22.099205 12.400865 C 20.701869 13.8994 19.098242 14.701214 17.1988 14.701214 C 15.299358 14.701214 13.699623 13.8994 12.298396 12.400865 C 10.90106 10.898438 10.301646 9.201395 10.301646 7.298061 C 10.301646 5.398619 10.998368 3.798884 12.298396 2.498856 C 13.598424 1.198828 15.299358 0.498214 17.1988 0.498214 C 19.098242 0.498214 20.701869 1.198828 22.099205 2.498856 C 23.500432 3.798884 24.099847 5.398619 24.099847 7.298061 Z M 24.099847 7.298061 " transform="matrix(1.003584,0,0,1.003584,0.137993,0)"/>
</g>
</svg>'
human <- function(){
t <- tempfile(); writeLines(humanSvg, t); pic <- readPicture(t)
}
#' Function to convert factors to numeric
#' @param x A vector or data.frame
#' @return A data.frame with numeric vectors, or a factors (if x is non-numeric), or a single vector if x was a vector
#'
#' @author mg14
#' @export
numericalize <- function(x){
if(class(x)=="data.frame")
as.data.frame(sapply(x, .numericalize, simplify=FALSE))
else
.numericalize(x)
} # Convert factors to numeric
.numericalize <- function(x) if(class(x)!='factor') return(x) else if(all(is.na(as.numeric(levels(x))))) return(x) else return(as.numeric(as.character(x)))
asum <- function(x, dim, ...){
apply(x, setdiff(seq_along(dim(x)), dim), sum, ...)
}
nn.poisglm <- function(D, P, maxIter = 1e4, tol=1e-6) {
n <- nrow(D)
m <- ncol(D)
s <- ncol(P)
tP <- t(P)
rP <- rep(colSums(P), m)
D <- as.matrix(D)
E1 <- E2 <- matrix(runif(s * m, 1e-3, 1), ncol = m)
err <- 2*tol
KD <- function(D, X) sum(D * log(D/X) - D + X)
k1 <- Inf
k2 <- 0
iter <- 1
while ( (iter < maxIter & err > tol & k1 - k2 > tol) | iter <= 10) {
E1 <- E2
k1 <- k2
E2 <- E1 * (tP %*% (D/(P %*% (E1 + .Machine$double.eps))))/rP
k2 <- KD(D, P%*%E2)
iter <- iter + 1
err <- mean(abs(E2 - E1)/(E1+.Machine$double.eps), na.rm=TRUE)
if(iter %% 100 == 0) cat(round(-log10(err)))
}
cat("\n")
if(iter == maxIter) warning(paste("No convergence after",iter, "iterations."))
E2
}
mindist <- function(x, mindist=0.03){
y <- x
while( any((diff(y))<mindist) ) {
ydiff = (diff(y))
runs = rle(ydiff<mindist)
aux = which(ydiff<mindist)[1]
l = runs$lengths[which(runs$values)[1]]
cl = seq(aux,aux+l)
# New suggested values for the cluster chosen
centr = median(y[cl])
y[cl] = seq(0, mindist*(length(cl)-1)+0.001, length.out=length(cl))
y[cl] = y[cl] - median(y[cl]) + centr
}
return(y)
}
poisl <- function(X, theta, Y){
Yp <- X %*% theta
Y %*% log(Yp) - sum(Yp) - sum(lfactorial(Y))
}
poisdl <- function(X, theta, Y){
lambda <- as.numeric(X %*% theta)
t(Y / lambda - 1) %*% X
}
poisI <- function(X, theta, Y){
lambda <- as.numeric(X %*% theta)
t(Y/lambda^2 * X) %*% X
}
ssnmf <- function(D, S, s=0, maxIter=100, minE = 0, whichS = 1:ncol(D)){
n <- nrow(D)
o <- ncol(S)
m <- ncol(D)
P <- cbind(S, if(s>0) matrix(runif(n*s,0,1), ncol=s) else NULL)
P <- P/rep(colSums(P), each=nrow(P))
colnames(P) <- c(colnames(S), if(s >0) paste0("N.",1:s) else NULL)
E <- matrix(runif((s+o)*m, 0,1), ncol=m)
E <- E * (t(P)%*% (D / (P %*% E))) / rep(colSums(P), m)
iter <- 1
while(iter < maxIter){
P <- P * ((D / (P %*% E)) %*% t(E)) / rep(rowSums(E), each=n)
P[,1:o] <- S
P <- P/rep(colSums(P), each=nrow(P))
E <- E * (t(P)%*% (D / (P %*% E))) / rep(colSums(P), m)
E[E/rep(colSums(E), each=nrow(E)) < minE] <- 0
E[-(1:o),setdiff(1:ncol(E), whichS)] <- 0
E <- E * rep(colSums(D)/colSums(E), each=nrow(E))
iter <- iter + 1
}
list(E=E,P=P)
}
nmSolve <- function(D, P, maxIter = 500, tol=1e-3) {
n <- nrow(D)
m <- ncol(D)
s <- ncol(P)
tP <- t(P)
rP <- rep(colSums(P), m)
D <- as.matrix(D)
E1 <- E2 <- matrix(runif(s * m, 1e-3, 1), ncol = m)
err <- 2*tol
iter <- 1
while (iter < maxIter & err > tol) {
E1 <- E2
E2 <- E1 * (tP %*% (D/(P %*% (E1 + .Machine$double.eps))))/rP
iter <- iter + 1
err <- mean(abs(E2 - E1)/(E1+.Machine$double.eps), na.rm=TRUE)
if(iter %% 100 == 0) cat(round(-log10(err)))
}
cat("\n")
if(iter == maxIter) warning(paste("No convergence after",iter, "iterations."))
E2
}
roundProp <- function(x, p=100){
y <- floor(x)
d <- x-y
o <- order(d, decreasing=TRUE)
i <- 1
while(sum(y) < p){
y[o[i]] <- y[o[i]] + 1
i <- i+1
}
return(y)
}
scatterpie <- function(x,y, p, r, xlab="",ylab="", circles=FALSE, lwd.circle=rep(1, length(x)), lty.circle=rep(1, length(x)), add=FALSE ,...){
if(!add) plot(x,y, xlab=xlab, ylab=ylab, pch=NA)
for(i in seq_along(x)){
.pie(p[i,], x0=x[i], y0=y[i], radius=r[i], add=TRUE, ...)
if(circles){
u <- par("usr")
pr <- (u[2]-u[1])/(u[4]-u[3])
fr <- par("pin")[1]/par("pin")[2]
polygon(x[i] + cos(seq(0,2*pi, l=100)) * r[i], y[i]+sin(seq(0,2*pi, l=100)) * r[i] / pr * fr, col=NA, lty=lty.circle[i], lwd=lwd.circle[i])
}
}
}
.pie <- function (x, x0=0, y0=0, labels = names(x), edges = 200, radius = 0.8, clockwise = FALSE,
init.angle = if (clockwise) 90 else 0, density = NULL, angle = 45,
col = NULL, border = NULL, lty = NULL, main = NULL, add=FALSE,...)
{
if (!is.numeric(x) || any(is.na(x) | x < 0))
stop("'x' values must be positive.")
if (is.null(labels))
labels <- as.character(seq_along(x))
else labels <- as.graphicsAnnot(labels)
x <- c(0, cumsum(x)/sum(x))
dx <- diff(x)
nx <- length(dx)
if(!add) plot.new()
pin <- par("pin")
xlim <- ylim <- c(-1, 1) + c(x0,y0)
if (pin[1L] > pin[2L])
xlim <- (pin[1L]/pin[2L]) * xlim
else ylim <- (pin[2L]/pin[1L]) * ylim
if(!add) dev.hold()
on.exit(dev.flush())
if(!add) plot.window(xlim, ylim, "", asp = 1)
if (is.null(col))
col <- if (is.null(density))
c("white", "lightblue", "mistyrose", "lightcyan",
"lavender", "cornsilk")
else par("fg")
col <- rep(col, length.out = nx)
border <- rep(border, length.out = nx)
lty <- rep(lty, length.out = nx)
angle <- rep(angle, length.out = nx)
density <- rep(density, length.out = nx)
twopi <- if (clockwise)
-2 * pi
else 2 * pi
u <- par("usr")
pr <- (u[2]-u[1])/(u[4]-u[3])
fr <- par("pin")[1]/par("pin")[2]
t2xy <- function(t) {
t2p <- twopi * t + init.angle * pi/180
list(x = radius * cos(t2p), y = radius * sin(t2p)/pr*fr)
}
for (i in 1L:nx) {
n <- max(2, floor(edges * dx[i]))
P <- t2xy(seq.int(x[i], x[i + 1], length.out = n))
P$x <- P$x + x0
P$y <- P$y + y0
polygon(c(P$x, x0), c(P$y, y0), density = density[i], angle = angle[i],
border = border[i], col = col[i], lty = lty[i])
P <- t2xy(mean(x[i + 0:1]))
lab <- as.character(labels[i])
if (!is.na(lab) && nzchar(lab)) {
lines(c(1, 1.05) * P$x + x0, c(1, 1.05) * P$y + y0)
text(1.1 * P$x + x0, 1.1 * P$y + y0, labels[i], xpd = TRUE,
adj = ifelse(P$x < 0, 1, 0), ...)
}
}
title(main = main, ...)
invisible(NULL)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.