#Author: Sarah P. Flanagan
#Purpose: Plot genome-wide statistics and perform other population genomics visualizations
#Package: gwscaR
#' Plot genome-wide statistics from a data frame.
#' @param fst.dat The data.frame containing at least three columns: the statistic to be plotted, the chromosome ID, and the BP ID for each locus. Each row is a locus.
#' @param pt.cols The color of the points. The default is c("darkgrey","lightgrey").
#' @param fst.name The name of the column containing the statistic to be plotted. Default is "Fst".
#' @param chrom.name The name of the column containing the chromosome information for each locus. Default is "Chrom".
#' @param bp.name The name of the column containing the basepair information for each locus. Default is "BP".
#' @param axis.size The value of cex.axis passed to the plotting of the y-axis. Default is 0.5. Set axis.size = 0 to suppress plotting of the y-axis.
#' @param scaffs.to.plot A vector or list containing the names of the chromosomes/scaffolds you want to include in the order in which you would like them to be plotted.
#' @param y.lim The limits for the y-axis.
#' @param pt.cex The size of the points
#' @param xlabels Either a boolean (TRUE) or a vector of labels for the x-axis
#' @param xlab.indices If the xlabels are not the same as the names of the chromosomes in fst.dat, then provide xlab.indices, which is a vector of integers that are the indices in scaffs.to.plot for which chromosome sets you would like to label.
#' @param scaffold.widths A data.frame with first column naming the chromosomes/linkage groups, second column having the last number of loci on the chromosome for indexing purposes.
#' e.g.:
#' Chrom End
#' LG1 3945850
#' LG2 435215
#' This parameter is used to make multiple plots from the same overall dataset have the same widths.
#' @param ... plot parameters
#' @return The fst.dat data frame with new values in BP, scaled to be the sequence in which points are plotted based on their position on the chromosome and the order the chromosomes are plotted in.
#' @seealso Flanagan & Jones 2017
#' @export
fst.plot<-function(fst.dat,scaffold.widths=NULL,scaffs.to.plot=NULL,
fst.name="Fst", chrom.name="Chrom", bp.name="BP",
y.lim=NULL,xlabels=NULL,xlab.indices=NULL,
axis.size=0.5,pt.cols=c("darkgrey","lightgrey"),pt.cex=0.5,...){
#make sure fsts are numeric
fst.dat<-fst.dat[!is.na(fst.dat[,fst.name]),]
fst.dat[,fst.name]<-as.numeric(as.character(fst.dat[,fst.name]))
#if scaffold widths aren't provided, we add them ourselves
if(is.null(scaffold.widths)){
scaffold.widths<-data.frame(Chrom=levels(as.factor(fst.dat[,chrom.name])),
Max=tapply(as.numeric(as.character(fst.dat[,bp.name])),
fst.dat[,chrom.name],max))
}
colnames(scaffold.widths)<-c("Chrom","Max")
#if no specific scaffolds are provided, we plot them all
if(is.null(scaffs.to.plot)){
scaffs.to.plot<-scaffold.widths[,1]#unique(fst.dat[,chrom.name])
}
#set up new dataframe
new.dat<-data.frame(stringsAsFactors = F)
last.max<-0
for(i in 1:length(scaffs.to.plot)){
#pull out the data for this scaffold
if(nrow(scaffold.widths[scaffold.widths$Chrom %in% scaffs.to.plot[i],])>0){ #sanity check
chrom.dat<-fst.dat[fst.dat[,chrom.name] %in% scaffs.to.plot[i],]
if(nrow(chrom.dat)>0){
chrom.dat$plot.pos<-as.numeric(as.character(chrom.dat[,bp.name]))+last.max
new.dat<-rbind(new.dat,chrom.dat)
#last.max<-max(chrom.dat$plot.pos)+
# as.numeric(scaffold.widths[scaffold.widths[,1] %in% scaffs.to.plot[i],2])
}
last.max<-last.max+
as.numeric(scaffold.widths[scaffold.widths$Chrom %in% scaffs.to.plot[i],2])
}else{
print(paste(scaffs.to.plot[i], "has no designated width."))
}
}
#make sure everything is the correct class
new.dat$plot.pos<-as.numeric(as.character(new.dat$plot.pos))
new.dat[,chrom.name]<-as.factor(as.character(new.dat[,chrom.name]))
scaffs.to.plot<-as.factor(as.character(scaffs.to.plot))
#determine the axis limits
if(last.max>=max(new.dat$plot.pos)){
x.max<-last.max
}else{
x.max<-max(new.dat$plot.pos)#max(scaffold.widths[scaffold.widths[,1] %in% scaffs.to.plot,2],na.rm=T)
}
x.min<-min(new.dat$plot.pos,na.rm=T)
if(is.null(y.lim)){
y.max<-max(fst.dat[,fst.name])+0.1*max(fst.dat[,fst.name])
y.min<-min(fst.dat[,fst.name])-0.1*min(fst.dat[,fst.name])
if(min(fst.dat[,fst.name]) < 0) {
y.min<-min(fst.dat[,fst.name]) - 0.1*min(fst.dat[,fst.name])
} else {
y.min<-0
}
y.lim<-c(y.min,y.max)
}
#plot
colors<-data.frame(lg=as.character(new.dat[,chrom.name]),col=rep(pt.cols[1],nrow(new.dat)),
stringsAsFactors = F)
colors[as.numeric(factor(colors$lg,levels=scaffs.to.plot))%%2==0,"col"]<-pt.cols[2] #defining the levels maintains the order
graphics::plot(new.dat$plot.pos,new.dat[,fst.name],
xlim=c(x.min,x.max),ylim=y.lim,...,
cex=pt.cex, col=colors$col,
bty="n",axes=F, xlab="", ylab="")
#optionally add the y-axis
if(axis.size>0){
graphics::axis(2, at = seq(round(y.lim[1],2),round(y.lim[2],2),
round((y.lim[2]-y.lim[1])/2, digits=2)),
ylim =y.lim, pos=0,
labels=seq(round(y.lim[1],2),round(y.lim[2],2),
round((y.lim[2]-y.lim[1])/2, digits=2)),
las=1, xlab="", ylab="", cex.axis=axis.size)
}
#optionally add the x-axis
if(!is.null(xlabels)){
#is this a boolean or a list?
if(class(xlabels)=="logical"){
xlabels<-as.character(scaffs.to.plot)
}
#if there are not indices, they should just be the ones that match scaffolds
if(is.null(xlab.indices)){
xlab.indices<-which(scaffs.to.plot %in% xlabels)
}
for(i in 1:length(xlabels)){
graphics::text(x=stats::median(new.dat[new.dat[,chrom.name] %in% scaffs.to.plot[xlab.indices[i]],"plot.pos"]),
y=(y.lim[1]-((y.lim[2]-y.lim[1])/20)),
labels=xlabels[i],adj=1,xpd=TRUE,cex=as.numeric(axis.size))
}
}
return(new.dat)
}
#' Plot genome-wide statistics from a data frame with rectangles.
#' @note This fst.plot function is more efficient but this will work for most cases (might be buggy in some cases)
#' @param fst.dat The data.frame containing at least three columns: the statistic to be plotted, the chromosome ID, and the BP ID for each locus. Each row is a locus.
#' @param ci.dat A vector containing two values, upper and lower cutoff values (in that order).
#' Points above or below the cutoffs will be colored based on sig.col. Default is NULL, which turns this option off.
#' @param sig.col A vector containing two color values, the first for points above the upper cutoff value and the second for points below the lower cutoff value.
#' The defaults are red and yellow.
#' @param pt.col The color of the points. The default is grey7.
#' @param pt.cex The size of the points. Default is 0.5.
#' @param fst.name The name of the column containing the statistic to be plotted. Default is "Fst".
#' @param chrom.name The name of the column containing the chromosome information for each locus. Default is "Chrom".
#' @param bp.name The name of the column containing the basepair information for each locus. Default is "BP".
#' @param axis.size The value of cex.axis passed to the plotting of the y-axis. Default is 0.5. Set axis.size = 0 to suppress plotting of the y-axis.
#' @param scaffold.order A vector or list containing the names of the chromosomes/scaffolds in the order in which you would like them to be plotted.
#' @param groups A vector indicating which chromosomes to plot (generally to exclude those scaffolds not found in this particular set of statistics due to pruning/filters)
#' @param print.names A TRUE/FALSE value indicating whether chromosome/scaffold IDs should be printed beneath the x-axis.
#' @param y.lim The limits for the y-axis.
#' @param group.boundaries A data.frame with first column naming the chromosomes/linkage groups, second column having the last number of loci on the chromosome for indexing purposes.
#' e.g.:
#' Chrom End
#' LG1 3945850
#' LG2 435215
#' This parameter is used to make multiple plots from the same overall dataset have the same widths.
#' @param ... plot parameters
#' @return The fst.dat data frame with new values in BP, scaled to be the sequence in which points are plotted based on their position on the chromosome and the order the chromosomes are plotted in.
#' @seealso Flanagan & Jones 2017
#' @export
fst.plot.rect<-function(fst.dat,ci.dat=NULL, sig.col=c("red","yellow"),pt.col="grey7",
fst.name="Fst", chrom.name="Chrom", bp.name="BP",axis.size=0.5,
scaffold.order=NULL,groups=NULL,print.names=FALSE,y.lim=NULL,
group.boundaries=NULL,pt.cex=0.5,...){
if(!is.null(scaffold.order)){
scaff.ord<-scaffold.order$component_id
lgs<-scaffold.order$object
} else{
if(!is.null(group.boundaries)){
scaff.ord<-unique(group.boundaries[,1])
lgs<-scaff.ord
}else{
scaff.ord<-unique(fst.dat[,chrom.name])
lgs<-scaff.ord
}
}
if(!is.null(groups)){
lgs<-groups
scaff.ord<-groups
}
fst.dat[,fst.name]<-as.numeric(as.character(fst.dat[,fst.name]))
these.scaffs<-scaff.ord[scaff.ord %in% unique(fst.dat[,chrom.name])]
all.scaff<-split(fst.dat, factor(fst.dat[,chrom.name]))
#keep only the scaffolds you're plotting
if(!is.null(groups))
{
all.scaff<-all.scaff[names(all.scaff) %in% groups]
}
group.boundaries<-group.boundaries[scaff.ord,]
last.max<-0
rect.xs<-NULL
addition.values<-0
xlist<-NULL
#set up the rectangles
for(i in 1:length(scaff.ord)){
if(is.null(group.boundaries)){
new.max<-last.max+nrow(all.scaff[[as.character(scaff.ord[i])]])
}else{
new.max<-as.numeric(group.boundaries[(group.boundaries[1,]) %in% scaff.ord[i],2])+
as.numeric(as.character(last.max))
}
rect.xs<-rbind(rect.xs,c(last.max, new.max))
rownames(rect.xs)[i]<-scaff.ord[i]
addition.values<-c(addition.values, new.max)
last.max<-new.max
}
rownames(rect.xs)<-scaff.ord
for(i in 1:length(all.scaff)){
if(is.null(group.boundaries)){
#then this is done by relative position only
all.scaff[[i]]<-all.scaff[[i]][order(all.scaff[[i]][,bp.name]),]
all.scaff[[i]][,bp.name]<-
seq(rect.xs[rownames(rect.xs) %in% names(all.scaff)[i],1]+1,
rect.xs[rownames(rect.xs) %in% names(all.scaff)[i],1]+nrow(all.scaff[[i]]),1)
}else{
#then it takes genome space into account
all.scaff[[i]]<-
all.scaff[[i]][order(all.scaff[[i]][,bp.name]),]
all.scaff[[i]][,bp.name]<-
as.numeric(rect.xs[rownames(rect.xs) %in% names(all.scaff)[i],1])+
as.numeric(as.character(all.scaff[[i]][,bp.name]))
}
}
#change BP to plot
x.max<-max(rect.xs,na.rm=T)
x.min<-min(rect.xs,na.rm=T)
if(is.null(y.lim)){
y.max<-max(fst.dat[,fst.name])+0.1*max(fst.dat[,fst.name])
y.min<-min(fst.dat[,fst.name])-0.1*min(fst.dat[,fst.name])
if(min(fst.dat[,fst.name]) < 0) {
y.min<-min(fst.dat[,fst.name]) - 0.1*min(fst.dat[,fst.name])
} else {
y.min<-0
}
y.lim<-c(y.min,y.max)
}
displacement<-y.lim[1]-((y.lim[2]-y.lim[1])/30)
graphics::plot(c(x.min,x.max),y.lim,xlim=c(x.min,x.max),
ylim=y.lim, col=pt.col,...,
bty="n",type="n", axes=F, xlab="", ylab="")
for(i in 1:nrow(rect.xs)){
if(i%%2 == 0) {
rect.color<-"white"
} else {
rect.color<-"gray75"
}
graphics::rect(rect.xs[i,1],y.lim[1],rect.xs[i,2],y.lim[2],
col=rect.color, border=NA)
if(print.names==T){
graphics::text(x=mean(rect.xs[i,]),
y=displacement,labels=rownames(rect.xs)[i],
adj=1,xpd=T,srt=45)
}
}
for(i in 1:length(all.scaff)){
graphics::points(all.scaff[[i]][,bp.name],
all.scaff[[i]][,fst.name],
pch=19, cex=pt.cex,col=pt.col,
xlim=c(x.min,x.max),ylim=y.lim)
if(!is.null(ci.dat)){
temp.sig<-all.scaff[[i]][all.scaff[[i]][,fst.name] >= ci.dat[1],]
graphics::points(temp.sig[,bp.name], temp.sig[,fst.name],
col=sig.col[1], pch=19, cex=0.5)
temp.sig<-all.scaff[[i]][all.scaff[[i]][,fst.name] <= ci.dat[2],]
graphics::points(temp.sig[,bp.name], temp.sig[,fst.name],
col=sig.col[2], pch=19, cex=0.5)
}
}
if(axis.size>0){
graphics::axis(2, at = seq(round(y.lim[1],2),round(y.lim[2],2),
round((y.lim[2]-y.lim[1])/2, digits=2)),
ylim =y.lim, pos=0,
labels=seq(round(y.lim[1],2),round(y.lim[2],2),
round((y.lim[2]-y.lim[1])/2, digits=2)),
las=1,tck = -0.01, xlab="", ylab="", cex.axis=axis.size)
}
xes<-do.call("rbind",all.scaff)
return(xes)
}
#' This plots the output of structure
#' @param structure.out A data.frame containing the structure output.
#' @param k The value of K used in Structure
#' @param pop.order The order in which the populations are plotted.
#' @param filename The name of the jpeg for the output file.
#' @param make.file A boolean (TRUE/FALSE) indicating whether the jpeg should be created
#' @param plot.new A boolean (TRUE/FALSE) indicating whether this is being added to an existing plot
#' @param colors A list of colors for the structure colors (if not provided, defaults to rainbow palette)
#' @param xlabel A boolean (TRUE/FALSE) indicating whether x labels should be plotted (default is TRUE)
#' @param ylabel An optional label for the y-axis.
#' @param xlabcol An optional vector of colors for x-axis labels
#' @param lab.cex Size of the labels
#' @export
plotting.structure<-function(structure.out, k, pop.order,
filename=paste("str.k",k,".jpeg",sep=""),make.file=TRUE,lab.cex=1,
plot.new=TRUE,colors=NULL,xlabel=TRUE,ylabel=NULL,xlabcol=NULL){
str.split<-split(structure.out,structure.out[,1])
if(is.null(colors)){
bar.colors<-grDevices::rainbow(k,s=0.5)
} else {
bar.colors<-colors
}
if(is.null(xlabcol)){
xlabcol<-rep("black",length(str.split))
}
if(make.file==TRUE){
grDevices::jpeg(filename,width=7, height=1.25, units="in", res=300)
graphics::par(mfrow=c(1,length(str.split)),mar=c(1,0,0,0), oma=c(1,0,0,0),cex=0.5)
} else {
if(plot.new==TRUE){
graphics::par(mfrow=c(1,length(str.split)),mar=c(1,0,0,0), oma=c(1,0,0,0),cex=0.5)
}
}
for(i in 1:length(str.split)){
pop.index<-pop.order[i]
graphics::barplot(height=as.matrix(t(str.split[[pop.index]][,-1])),
beside=FALSE, space=0, border=NA, col=bar.colors,
xlab="", ylab="", xaxt='n', yaxt='n')#, new=plot.new)
if(xlabel==TRUE){
graphics::mtext(pop.index, 1, line=0.5, cex=lab.cex, outer=F,col=xlabcol[i])}
if(!is.null(ylabel)){
if(i == 1) { graphics::mtext(ylabel,2,cex=lab.cex) }
}
}
if(make.file==TRUE) {grDevices::dev.off()}
}
#' Plot covariance matrix from treemix
#' @param stem The filename basic stem for that run
#' @param poporder The list of populations in the order to plot them.
#' @param colors Default to blue, yellow, and red.
#' @param ... print parameters
#' @return cp The covariance matrix in the plotting order.
#' @export
treemix.cov.plot<-function(stem,poporder,colors=c("blue","yellow","red"),...){
if("package:lattice" %in% search() == FALSE || "package:grid" %in% search()==FALSE){
stop("ERROR: packages lattice and grid must be loaded")
}
cov<-utils::read.table(gzfile(paste(stem, ".cov.gz", sep = "")), as.is = T, head = T, quote = "", comment.char = "")
#reorder
covplot = data.frame(matrix(nrow = nrow(cov), ncol = ncol(cov)))
for(i in 1:length(poporder)){
for( j in 1:length(poporder)){
covplot[i, j] = cov[which(names(cov)==poporder[i]), which(names(cov)==poporder[j])]
}
}
rownames(covplot)<-poporder
colnames(covplot)<-poporder
cp<-as.matrix(covplot)
cp[lower.tri(cp)]<-NA
cp[upper.tri(cp)]<-covplot[upper.tri(covplot)]
#set up the colors
pal<-grDevices::colorRampPalette(colors)
ncol=80
cols<-pal(ncol)
cpl<-lattice::levelplot(cp,col.regions=cols,alpha.regions=0.7,
scales = list(x=list(rot=90),tck = 0),xlab="",ylab="")
print(cpl,...)
lattice::trellis.focus("legend", side="right", clipp.off=TRUE, highlight=FALSE)
grid::grid.text("covariance", 0.2, 0, hjust=0.5, vjust=1.2)
lattice::trellis.unfocus()
return(cp)
}
#' Add a legend to the margins of a multiplot figure
#' @param ... legend parameters
#' @note Modified from https://stackoverflow.com/questions/3932038/plot-a-legend-outside-of-the-plotting-area-in-base-graphics
#' @examples
#' plot(rnorm(50), rnorm(50), col=c("steelblue", "indianred"), pch=c(15,19))
#' outer.legend("top",legend=c("A","B"),pch=c(15,19),col=c("steelblue", "indianred"),ncol=2)
#' @export
outer.legend <- function(...) {
opar <- graphics::par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0),
mar=c(0, 0, 0, 0), new=TRUE)
on.exit(graphics::par(opar))
graphics::plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
graphics::legend(...)
}
#' Create violin plots with
#' @param x the data to plot
#' @param ... legend parameters
#' @param range The range for plotting
#' @param h The height of the density estimator (in sm.density), if ommitted h will be set to an optimum
#' @param ylim The y-axis limits
#' @param names Names to plot (default is NULL)
#' @param horizontal Whether the plot should be horizontal (defaults to FALSE)
#' @param col Color for violin plots
#' @param border Border color
#' @param lty Line type
#' @param lwd line width
#' @param rectCol rectangle color
#' @param colMed median point color
#' @param pchMed median point size
#' @param at The location for things
#' @param add Whether the plot gets added to others (defaults to FALSE)
#' @param wex Size of something
#' @param drawRect Whether to add rectangles to violin plots (defaults to TRUE)
#' @param plot.axes A logical value that determines whether or not to add axes (defaults to TRUE)
#' @param axis.box A logical value that determines whether the axes should be a box (defaults to FALSE)
#' @param plot.ann A logical value that determines whether the annotations should be plotted (defaults to TRUE)
#' @param axis.bty The axis box type (default is L)
#' @return Optionally returns the summary statistics of the data.
#' @note Modified from vioplot() in library(vioplot)
#' @examples
#' library(sm)
#' mu<-2
#' si<-0.6
#' bimodal<-c(rnorm(1000,-mu,si),rnorm(1000,mu,si))
#' uniform<-runif(2000,-4,4)
#' normal<-rnorm(2000,0,3)
#' gwsca.vioplot(bimodal,uniform,normal,col=c("red","blue","green"))
#' @export
gwsca.vioplot <- function(x,...,range=1.5,h=NULL,ylim=NULL,names=NULL, horizontal=FALSE,
col="magenta", border="black", lty=1, lwd=1, rectCol="black", colMed="white", pchMed=19, at, add=FALSE, wex=1,
drawRect=TRUE,plot.axes=TRUE,axis.box=FALSE,plot.ann=TRUE,axis.bty="l")
{
requireNamespace("sm")
# process multiple datas
datas <- list(x,...)
n <- length(datas)
if(missing(at)) at <- 1:n
# pass 1
#
# - calculate base range
# - estimate density
#
# setup parameters for density estimation
upper <- vector(mode="numeric",length=n)
lower <- vector(mode="numeric",length=n)
q1 <- vector(mode="numeric",length=n)
q3 <- vector(mode="numeric",length=n)
med <- vector(mode="numeric",length=n)
base <- vector(mode="list",length=n)
height <- vector(mode="list",length=n)
baserange <- c(Inf,-Inf)
# global args for sm.density function-call
args <- list(display="none")
if (!(is.null(h)))
args <- c(args, h=h)
for(i in 1:n) {
data<-datas[[i]]
# calculate plot parameters
# 1- and 3-quantile, median, IQR, upper- and lower-adjacent
data.min <- min(data)
data.max <- max(data)
q1[i]<-stats::quantile(data,0.25)
q3[i]<-stats::quantile(data,0.75)
med[i]<-stats::median(data)
iqd <- q3[i]-q1[i]
upper[i] <- min( q3[i] + range*iqd, data.max )
lower[i] <- max( q1[i] - range*iqd, data.min )
# strategy:
# xmin = min(lower, data.min))
# ymax = max(upper, data.max))
#
est.xlim <- c( min(lower[i], data.min), max(upper[i], data.max) )
# estimate density curve
smout <- do.call("sm.density", c( list(data, xlim=est.xlim), args ) )
# calculate stretch factor
#
# the plots density heights is defined in range 0.0 ... 0.5
# we scale maximum estimated point to 0.4 per data
#
hscale <- 0.4/max(smout$estimate) * wex
# add density curve x,y pair to lists
base[[i]] <- smout$eval.points
height[[i]] <- smout$estimate * hscale
# calculate min,max base ranges
t <- range(base[[i]])
baserange[1] <- min(baserange[1],t[1])
baserange[2] <- max(baserange[2],t[2])
}
# pass 2
#
# - plot graphics
# setup parameters for plot
if(!add){
xlim <- if(n==1)
at + c(-.5, .5)
else
range(at) + min(diff(at))/2 * c(-1,1)
if (is.null(ylim)) {
ylim <- baserange
}
}
if (is.null(names)) {
label <- 1:n
} else {
label <- names
}
# setup colors and borders
if(length(col) < n){
col<-c(col,rep(col,n-length(col)))
}
if(length(border)<n){
border<-c(border,rep(border,n-length(border)))
}
if(length(colMed)<n){
colMed<-c(colMed,rep(colMed,n-length(colMed)))
}
boxwidth <- 0.05 * wex
# setup plot
if(!add)
graphics::plot.new()
if(!horizontal) {
if(!add){
graphics::plot.window(xlim = xlim, ylim = ylim)
if(plot.axes){
if(plot.ann){
graphics::axis(2)
graphics::axis(1,at = at, label=label )
}else{
graphics::axis(2,labels=F)
graphics::axis(1,at = at, labels=F )
}
}
}
if(axis.box){ graphics::box() }
for(i in 1:n) {
# plot left/right density curve
graphics::polygon( c(at[i]-height[[i]], rev(at[i]+height[[i]])),
c(base[[i]], rev(base[[i]])),
col = col[i], border=border[i], lty=lty, lwd=lwd)
if(drawRect){
# plot IQR
graphics::lines( at[c( i, i)], c(lower[i], upper[i]) ,lwd=lwd, lty=lty)
# plot 50% KI box
graphics::rect( at[i]-boxwidth/2, q1[i], at[i]+boxwidth/2, q3[i], col=rectCol)
# plot median point
graphics::points( at[i], med[i], pch=pchMed, col=colMed )
}
}
}
else {
if(!add){
graphics::plot.window(xlim = ylim, ylim = xlim,bty=axis.bty)
if(plot.axes){
if(plot.ann){
graphics::axis(2)
graphics::axis(1,at = at, label=label )
}else{
graphics::axis(2,labels=F)
graphics::axis(1,at = at, labels=F )
}
}
}
if(axis.box){ graphics::box() }
for(i in 1:n) {
# plot left/right density curve
graphics::polygon( c(base[[i]], rev(base[[i]])),
c(at[i]-height[[i]], rev(at[i]+height[[i]])),
col = col[i], border=border[i], lty=lty, lwd=lwd)
if(drawRect){
# plot IQR
graphics::lines( c(lower[i], upper[i]), at[c(i,i)] ,lwd=lwd, lty=lty)
# plot 50% KI box
graphics::rect( q1[i], at[i]-boxwidth/2, q3[i], at[i]+boxwidth/2, col=rectCol)
# plot median point
graphics::points( med[i], at[i], pch=pchMed, col=colMed )
}
}
}
invisible (list( upper=upper, lower=lower, median=med, q1=q1, q3=q3))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.