################# ################# kn <- function(x){ ( out <- knit(x) ) system(paste('pdflatex', out)) } ################# ################# show.colors <- function(x=22,y=30){ par(mfrow=c(1,1)) par(mai=c(.4,.4,.4,.4), oma=c(.2,0,0,.2)) plot(c(-1,x),c(-1,y), xlab='', ylab='', type='n', xaxt='n', yaxt='n', bty='n') for(i in 1:x){for(j in 1:y){ k <- y*(i-1)+j ; co <- colors()[k] rect(i-1, j-1, i, j, col=co, border=1)}} text(rep(-.5,y),(1:y)-.5, 1:y, cex=1.2-.016*y) text((1:x)-.5, rep(-.5,x), y*(0:(x-1)), cex=1.2-.022*x) title('col=colors()[n]') } ################# ################# add.leading0 <- function(v){ nc <- nchar(v) ; ncm <- max(nc) ; ncd <- ncm-nc to.add <- rep(0,L <- length(nc)) for(i in 1:L){ to.add[i] <- paste(rep(0,ncd[i]+1),collapse='') } substr(paste(to.add,v,sep=''),2,ncm) } ################# ################# within.group.id <- wgi <- function(v,g=NULL){ one.group <- function(v){ nl <- length(v) new <- which(diff(v)!=0) gs <- c(new[1],diff(c(new,nl))) rep(1:length(gs), gs) } out1 <- one.group(v) if(length(g)>0){ out1 <- NULL g <- factor(g) n.group <- length(unique(g)) for(i in 1:n.group){ vv <- v[g==unique(g)[i]] out1 <- c(out1,one.group(vv)) } } out2 <- NULL a <- c(which(diff(out1)!=0),length(v)) b <- c(a[1],diff(a)) for(i in 1:length(b)){out2 <- c(out2,1:b[i])} cbind(out1,out2) } ################# ################# logMain <- function(mn, mx, base){ lg <- seq(floor(log(mn, base)), ceiling(log(mx, base))) or <- base^lg data.frame(lg,or) } ################# ################# logBetween <- function(mn, mx, base){ # starting point is 10^exponent lo <- c(floor(log(mn, base)), ceiling(log(mx, base))) sp <- seq(min(lo), max(lo), by=1) bs <- as.list(base^sp) or <- unlist(lapply(bs, function(x) cumsum(rep(x,9))[-1])) lg <- log(or, base) data.frame(lg,or) } ################# ################# editRedCapData <- function(d){ ## Do d <- data ; rm(data) first. ## Don't forget the .r file that redcap creates; it has rm(list=ls()) in it! fac <- unlist(strsplit(names(d)[grep('factor', names(d))], '.factor')) for(i in fac){ nu <- which( names(d) == i) fa <- which( names(d) == paste(i,'factor',sep='.') ) la <- label(d[,nu]) d[,nu] <- d[,fa] label(d[,nu]) <- la names(d[,nu]) } d <- d[, -grep('factor', names(d))] names(d) <- capitalize(gsub('_+', '.', names(d))) d } ################# ################# groupSum <- function(v,g=NULL, Combined=TRUE, Test=FALSE, np=FALSE, MissingAsGroup=FALSE, r=2){ vg <- list(v) if(Test) Combined <- TRUE if(length(g) > 0) vg <- split(v, f=g) summaryNA <- function(x){ x <- x[ !is.na(x) ] sx <- summary(x)[1:6] tot <- sum(x) L <- length(x) s <- sd(x) se <- s/sqrt(L) c(sx,tot,L,s,se) } sm <- data.frame(do.call('rbind', lapply(vg, summaryNA))) if(Combined & length(g)>0){ comb <- summaryNA(v) sm <- rbind(sm,comb) row.names(sm)[nrow(sm)] <- 'Combined' } names(sm) <- c('Min','Q1','Med','Mean','Q3','Max','Total','N','SD','SE') sm <- subset(sm, select=c(N,Min,Q1,Med,Q3,Max,Mean,SD,SE)) out <- round(data.frame(sm), r) pv <- NULL if(Test){ if(length(vg)==2){ if(!np) pv <- c(NA,NA, t.test(vg[[1]], vg[[2]])$p.value) if( np) pv <- suppressWarnings( c(NA,NA, wilcox.test(vg[[1]], vg[[2]])$p.value) ) } out <- cbind(out, pv) } round(data.frame(out), r) } ################# ################# propPlot <- function(x, TrueFalse, howManyGroups=4, cutPoints=NULL){ # x is continuous. # TrueFalse if(is.null(cutPoints[1])) usr <- FALSE TrueFalse <- TrueFalse[!is.na(x)] x <- x[!is.na(x)] if(is.null(cutPoints[1])){ cox <- quantile(x, seq(0,1, by=1/howManyGroups)) co <- as.numeric(cox) co[1] <- co[1] - 1 xc <- cut(x,co) co[1] <- co[1] + 1 } if(!is.null(cutPoints[1])) co <- cutPoints xc <- cut(x,co) tab <- table( TrueFalse , xc) plot(range(x), c(0,1), xlab='', ylab='', type='n', xaxt='n', yaxs='i') # segments( min(x),0,min(x),1, col=grey(.9)) # segments( max(x),0,max(x),1, col=grey(.9)) # for(i in 1:(length(co)-2)) segments(co[i],0,co[i],1, col=grey(.9)) abline(h=(1:4)*0.2, col=grey(0.95)) cot <- formatC(co, format='f', digits=2) axis(1, at=co, label=cot) if(is.null(cutPoints[1])) axis(1, at=co, label=names(cox), line=1, tick=FALSE) midPoints <- colMeans(rbind( rev(rev(co)[-1]), co[-1])) width <- abs(midPoints-co[-1]) p <- as.numeric(prop.table(tab,2)[2,]) for(i in 1:length(midPoints)){ rect(midPoints[i]-width[i]/1.3,0,midPoints[i]+width[i]/1.3,p[i], col='royalblue', border='lightblue') text(midPoints[i], p[i], paste(tab[2,i], ' / ', colSums(tab)[i]), pos=3, cex=0.7 ) } box(bty='L') tab } ################# ################# multTime <- function(lap, multiplier, fmt=TRUE){ # lap is in 'Min.Sec' format. "8.47" means 8 minutes 47 seconds. Min <- floor(lap) ## interger portion. Sec <- (lap - Min) * 100 m <- Min + Sec/60 ## in decimal MINUTES. MM <- m * multiplier MIN <- floor(round(MM*100)/100) SEC <- round( 60*(MM - MIN)) out <- MIN+round(SEC/100,2) if(fmt) out <- paste(MIN, formatC(SEC, format='f', digit=0, width=2, flag=0), sep=':') out } ################# ################# addTime <- function(laps, fmt=TRUE){ # laps is in 'Min.Sec' format. "8.47" means 8 minutes 47 seconds. min <- floor(laps) sec <- (laps-min) * 100 MIN <- sum(min) SEC <- sum(sec) secM <- SEC %/% 60 secS <- round(SEC %% 60) multTime(MIN+secM+secS/100, 1, fmt) } ################# ################# FormatSum <- function(sumReverse, round.proportion=0, round.numeric=1){ ## Version 1.1 5/4/2015 # sumReverse is summary with method='reverse'. sum.cat <- function(outStats, round.proportion, round.numeric, var.name, pv, group.names){ # Categorical # outStats is sumReverse$stats[[i]] o <- outStats nm <- dimnames(outStats)[[2]] l <- length(group.names) if( ncol(outStats) != l){ o <- matrix(0, ncol=l, nrow=nrow(outStats)) row.names(o) <- row.names(outStats) for(k in 1:ncol(outStats)){ o[, which(group.names == nm[k]) ] <- outStats[,k] } } pt <- prop.table(o, margin=2) pp <- ifelse( any(pt>0.999, na.rm=TRUE), 3, 2) pp <- ifelse( round.proportion==0, pp, pp+1+round.proportion) P <- formatC( 100*c(pt), digits=round.proportion, width=pp, format='f') nn <- max(nchar(c(o))) N <- formatC( c(o), width=nn) o <- matrix(paste(N, ' (', P, '%)', sep=''), nrow=nrow(o)) o[, !group.names %in% nm] <- rep('', nrow(o)) o <- cbind(var.name, row.names(outStats), o, pv) if(nrow(o) == 2) o <- o[2,] if(!is.null(nrow(o))){ o[,1] <- c( var.name, rep('', nrow(o)-1)) o[,ncol(o)] <- c(pv, rep('', nrow(o)-1)) } o } sum.num <- function(outStats, quant, round.numeric, var.name, pv, group.names){ # Numerical # outStats is sumReverse$stats[[i]] # quant is sumReverse$quant o <- outStats l <- length(group.names) if( nrow(outStats) != l){ o <- matrix(0, ncol=ncol(outStats), nrow=l) row.names(o) <- group.names for(k in 1:nrow(outStats)){ o[ which(row.names(o) == row.names(outStats)[k]),] <- outStats[k,] } } med <- prettyNum( as.numeric( o[, quant==0.5 ]) )#, digits=round.numeric, format='f')#, width=wdt) loq <- prettyNum( as.numeric( o[, quant==0.25]) )#, digits=round.numeric, format='f')#, width=wdt) upq <- prettyNum( as.numeric( o[, quant==0.75]) )#, digits=round.numeric, format='f')#, width=wdt) num <- paste(med, ' (', loq, ', ', upq, ')', sep='') num[! row.names(o) %in% row.names(outStats) ] <- '' c(var.name, '', num, pv) } args2 <- list(round.proportion, sumReverse$quant) num.dig <- c(round.proportion, round.numeric) group.names <- names(sumReverse$group.freq) if(!is.null(sumReverse$testresults)){ pv <- sapply(sumReverse$testresults, function(x) x$P) pval <- paste('P=', formatC(pv, format='f', digits=3), sep='') pval[ pv<0.001 ] <- 'P<0.001' } else { pval <- rep(' ', length(sumReverse$stats) ) } ta <- NULL for(i in 1:length(sumReverse$stats)){ this.type <- sumReverse$type[i] this.fun <- list(sum.cat, sum.num)[[this.type]] ta0 <- this.fun(sumReverse$stats[[i]], args2[[this.type]], num.dig[this.type], sumReverse$labels[i], pval[i], group.names) ta <- rbind(ta, ta0) } ta <- rbind( c('', '', paste('(N=', as.numeric(sumReverse$group.freq), ')', sep=''), ''), ta) tad <- data.frame(ta, row.names=NULL) names(tad) <- c(' ', ' ', names(sumReverse$group.freq), 'P.value') tad } ################# ################# odds.ratio2 <- function(mat, alpha=.05, con.cor =T){ # mat is a 2x2 matirix # con.cor = continuity correction if(sum(dim(mat) == c(2,2))!=2){stop("\n", "Input must be a 2x2 matrix.", "\n")} cc <- c(0, .5)[con.cor+1] n11 <- mat[1,1] ; n12 <- mat[1,2] ; n21 <- mat[2,1] ; n22 <- mat[2,2] m11 <- n11 + cc ; m12 <- n12 + cc ; m21 <- n21 + cc ; m22 <- n22 + cc #theta <- (n11*n22)/(n12*n21) theta <- (m11*m22)/(m12*m21) # with continuity correction sighat <- sqrt(1/m11 + 1/m12 + 1/m21 + 1/m22) za <- qnorm(1-alpha/2) log.ci <- round(c(log(theta) - za * sighat, log(theta), log(theta) + za * sighat), 5) ori.ci <- round(exp(log.ci), 5) rev.or <- round((1/ori.ci)[c(3, 2, 1)], 5) Log.o.r. <- log(theta) ase <- sighat zz <- Log.o.r./ase pp <- 2*(1-pnorm(abs(zz))) summ <- data.frame(Log.Odds.Ratio=Log.o.r., ASE=ase, Z=zz, P=pp) o.r. <- data.frame(round(rbind(log.ci, ori.ci, rev.or), 4)) names(o.r.) <- c('LWR', 'Point', 'UPR') r.s <- apply(mat, 2, sum) c.s <- apply(mat, 1, sum) g.s <- sum(r.s) ans <- data.frame(cbind(rbind(mat, r.s), c(c.s, g.s)), row.names=c(" 1 ", " 2 ", "Total")) names(ans) <- c(" A ", " B ", "Total") list(summ, o.r., ans) } ################# ################# odds.ratio <- function(mat, alpha=.05, con.cor=T){ temp <- odds.ratio2(mat, alpha, con.cor) list(temp[[1]], temp[[2]])}