Nothing
R/default.plot.mvformula.R
In mvabund: Statistical Methods for Analysing Multivariate Abundance Data
Defines functions
default.plot.mvformula
################################################################################
## PLOT.mvformula: Plot functions for mvabund objects #
## (multivariate abundance data) #
################################################################################
default.plot.mvformula <- function( x,
y,
type="p",
main,
xlab,
ylab="Abundances",
col=if(type=="bx") "white" else "black",
fg="grey",
pch=1,
las=1,
write.plot="show",
filename="plot.mvabund",
n.vars=if(any(is.na(var.subset))) 12 else length(var.subset),
overall.main="",
data=parent.frame(),
var.subset=NA,
xvar.select=TRUE,
n.xvars=NA,
# n.xvars=if(any(is.na(xvar.subset))) min(3, sum(!is.interaction)) else length(xvar.subset),
xvar.subset=NA,
scale.lab="ss",
t.lab="t",
mfrow=NULL,
mfcol=NULL,
shift=TRUE,
border="black",
all.labels=FALSE,
keep.window=FALSE,
ask, ... ) {
#Kill Function if a BoxPlot has been passed
if (type=="bx") { stop("\nERROR: It is not ideal to use a boxplot for this data type\n") }
# If there is only one response variable this is changed to TRUE
# (overwritten independent of passed value)
dev <- dev.list()
dev.name <- getOption("device")
if(is.null(dev.name))
stop("Make sure that the 'device' option has a valid value, e.g. 'options(device = 'windows')'. Allowed values here are 'windows', 'win.graph', 'x11', 'X11'.")
# if(!(any(dev.name == c("windows", "win.graph", "x11", "X11")) ) )
# stop("Make sure that the 'device' option has a valid value, e.g. 'options(device = 'windows')'.
# Allowed values here are 'windows', 'win.graph', 'x11', 'X11'.")
m <- match.call(expand.dots = FALSE)
dots <- m$...
logWarn <- FALSE
if(!is.null(dots$log)){
# dots$log <- NULL
if(regexpr("y", dots$log ) !=-1) {
dots$log <- sub("y","", dots$log)
logWarn <- TRUE
}
}
usex <- TRUE
if (length(dots)>0){
# Delete arguments in ... that are defined lateron and cannot be used twice in
# the plot function.
deactive <- c("axes", "cex.axis")
deactivate <- (1:length(dots))[names(dots) %in% deactive ]
for (i in length(deactivate):1) {
dots[ deactivate[i] ] <- NULL #fixed the [[]] compile problem, could reduace functionallity.
}
dots <- lapply(dots, eval, parent.frame())
if(!is.null(dots$formula)) {
x <- dots$formula
dots$formula <- NULL
usex <- FALSE
}
subset.char <- dots$subset
} else subset.char <- NULL
miss.varsubset <- missing(var.subset)
# Change logical var.subset to numerical var.subset, if necessary. Note that
# NA values are logical as well, but should be excluded here.
if(!miss.varsubset){
if(is.logical(var.subset) & any(!is.na(var.subset)))
var.subset <- which(var.subset[!is.na(var.subset)])
}
miss.varsubset <- !is.numeric(var.subset) # If this function is called within
# another, the missing function could be tricked out.
var.subset <- as.vector(var.subset)
miss.xvarsubset <- missing(xvar.subset)
if(!miss.xvarsubset){
if(is.logical(xvar.subset) & any(!is.na(xvar.subset)))
var.subset <- which(xvar.subset[!is.na(xvar.subset)])
}
miss.xvarsubset <- any(is.na(xvar.subset)) || any(is.null(xvar.subset))
if ( (missing(x) & usex) ) stop("'formula' is missing or incorrect")
if(!inherits(x, "formula") ) {
# if no formula is passed a formula of supposedly passed
# matrices/vectors/data.frames will be tried to build.
# If there is only one mvabund object, it is assumed to be the response.
# If there are two mvabund objects, x (the first argument passed) is the
# response, !note that this is different in plot.mvabund
# The response cannot be a data.frame, when both objects are data.frames an
# error will be produced.
# If none is an mvabund object, x is the response, if it is not a data.frame.
if (missing(y)) { stop("A response and independent variables need to be provided for a formula.") }
if( (is.mvabund(y) & ! is.mvabund(x) ) | is.data.frame(x)) {
if(is.data.frame(y)) stop("a data.frame is not accepted as response for building the formula")
targs <- match.call(call = sys.call(which = 1), expand.dots = FALSE)
mvabund.formula <- x <- as.formula(paste(deparse(targs$y,width.cutoff = 500),"~",deparse(targs$x,width.cutoff = 500), sep=""))
# That might result in troubles if either targs$y or targs$x has the name "x"
# try to remove x (formula) earlier after giving it the name mvabund.formula
} else {
targs <- match.call(call = sys.call(which = 1), expand.dots = FALSE)
mvabund.formula <- x <- as.formula( paste(deparse(targs$x,width.cutoff = 500),"~", deparse(targs$y,width.cutoff = 500), sep="") )
}
# if(length(x) != 3) stop("'formula' is incorrect")
} else mvabund.formula <- x
resp <- attr(terms(mvabund.formula),"response")
term.labels <- attr(terms(mvabund.formula),"term.labels")
if(resp) {
if(length(mvabund.formula) != 3) stop("'formula' is incorrect")
mvabund.object.1<-as.matrix(model.response(model.frame( mvabund.formula, data=data)))
nam <- paste(deparse( mvabund.formula[[2]]), ".", sep="")
if(!is.null(colnames( mvabund.object.1 ))){
if(any(is.na(suppressWarnings(as.numeric(sub( nam ,"", colnames( mvabund.object.1 ))))))){
colnames( mvabund.object.1 ) <- sub( nam ,"", colnames( mvabund.object.1 ))
}
}
} else stop("'formula' has no response")
mvabund.object.1 <- mvabund(mvabund.object.1, neg=TRUE)
N <- NROW(mvabund.object.1) # number of sites
p <- NCOL(mvabund.object.1) # number of organism types
# Automatically correct n.vars if > p
n.vars <- min(n.vars,p)
if(length(term.labels)!=0){
foo.ext <- extend.x.formula(mvabund.formula,return.interaction=TRUE, extend.term=TRUE)
mvabund.formula <- foo.ext$formula
is.interaction <- foo.ext$is.interaction
# the length of this shows the number of x variables, including interaction
# term and factors, but factors only appear once and not like in model.matrix
# with (one column for every level) - 1
} else {
is.interaction <- c()
}
# Exclude Interactions.
if(any(is.interaction)){
mvabund.formula <- terms(mvabund.formula)[-which(is.interaction)]
print(2)
# terms.foo <- attr(model.frame(mvabund.formula), "terms")
# facs <- attr(terms.foo,"factors")[,is.interaction, drop = FALSE]
# datClass <- attr(terms.foo,"dataClasses")
# for(faci in 1:NCOL(facs)){
# if(any(datClass[facs[,faci]>0] == "factor"))
# stop("Formulas with interactions that include factors cannot yet be plotted")
# }
}
mf <- model.frame(mvabund.formula, data=data)
expl.data <- as.data.frame(mf[, -1])
# This is a matrix, ie factors are split to dummy variables
# not used for the actual plot, but for getting assign, etc
pExpl <- length(is.interaction) # NCOL(expl.data) + sum(is.interaction)
if(pExpl == 0) stop("The model is empty. There's nothing to plot.")
miss.xlab <- missing(xlab)
if (miss.xlab) xlab <- NULL
miss.main <- missing(main)
if (miss.main) {
nam <- paste(deparse( mvabund.formula[[2]]), ".", sep="")
# main <- colnames(as.matrix( model.frame( mvabund.formula)))[1:p]
main <- colnames(mvabund.object.1)
if(!is.null(main)){
if(any(is.na(suppressWarnings(as.numeric( sub(nam,"", main) ))))){
main <- colnames(mvabund.object.1) <- sub(nam,"", main)
} else colnames(mvabund.object.1) <- main
}
# main <- colnames(mvabund.object.1) <- colnames(as.matrix(
# model.frame( mvabund.formula)))[1:p]
n.mn <- nchar(main, type="chars")
if(any(n.mn > 14)) main[n.mn > 14] <- paste(substr(main[n.mn > 14], 1, 7),"...\n...",
substr(main[n.mn > 14], n.mn[n.mn > 14]-7, n.mn[n.mn > 14]) )
} else if (length(main)==p){
colnames(mvabund.object.1) <- main
} else if (!all(is.na(var.subset)) & length(main)==length(var.subset)){
colnames(mvabund.object.1)[var.subset]<- main
main <- colnames(mvabund.object.1)
} else if (length(main)==1) {
main <- rep(main,times =p)
} else stop("'main' must be a vector consisting of one or ", as.character(p), " characters" )
if (length(main)==1) { main <- rep(main,times= p) }
######### BEGIN edit var.subset, n.vars and mvabund.objects #########
# subset allows double variables
var.subset <- as.vector(var.subset)
var.subset.dim <- length(var.subset)
if (miss.varsubset) {
if(!any(length(ylab)==c(1,p))) stop("the length of 'ylab' is not correct")
sum.mvabund.object.1 <-
t(mvabund.object.1) %*% matrix(1,ncol=1,nrow=N)
# Find abundance ranks OF MVABUND.OBJECT.1.
var.subset<-order(sum.mvabund.object.1, decreasing = TRUE)
# Ensure no more than n.vars in var.subset.
if (n.vars < length(var.subset)) var.subset <- var.subset[1:n.vars]
var.subset.dim <- length(var.subset)
} else if ( p<max(var.subset) ) {
# Some dimension checks for subset, n.vars.
stop ("You have given an invalid 'var.subset'.")
} else if (n.vars!=var.subset.dim) {
n.vars <- var.subset.dim
}
mvabund.object.1 <- mvabund.object.1[,var.subset, drop=FALSE]
if(length(ylab)==p) {
ylab <- ylab[var.subset]
} else if (length(ylab)==1) ylab <- rep(ylab, times=n.vars)
main<- main[var.subset]
######### END edit var.subset, n.vars and mvabund.objects #########
mvabund.formulaUni <- formulaUnimva( mvabund.formula, var.subset=var.subset, split.x=TRUE, intercept=0 )
######### BEGIN edit xvar.subset, n.xvars etc #########
#DW changes, 30/10/14. 17/8: moved this outside of variable selection loop
#only define default n.xvars here, once is.interaction has been defined
if(is.na(n.xvars))
{
if(any(is.na(xvar.subset)))
n.xvars = min(3, sum(!is.interaction))
else
n.xvars = length(xvar.subset)
}
if(xvar.select & pExpl>1) #begin varaible selection (if required)
{
if(is.list(pch))
{
pch2 <- list()
for(i in 1:length(xvar.subset))
{
pch2[[xvar.subset[i] ]] <- pch[[i]]
}
pch <- pch2
}
if(is.list(col))
{
col2 <- list()
for(i in 1:length(xvar.subset))
{
col2[[xvar.subset[i]]] <- pch[[i]]
}
col<- col2
}
if(miss.xvarsubset)
{
# DW change, 17/8/15: changed to a warning if no n.xvars passed.
if(any(names(dots)==n.xvars)==FALSE)
warning(paste("No sQuote(n.xvars) passed so subset selection will be for",n.xvars,"variables"))
# Find the n.xvars independent variables with the highest average R^2.
xvar.subset <- best.r.sq(formula= mvabund.formula, # mvabund.formulaUni,
var.subset=var.subset, n.xvars= min(n.xvars, pExpl))
xvar.subset=xvar.subset$xs #to keep the indices only
if (is.null(xvar.subset))
{
xvar.subset <- 1:pExpl
warning("xvar.subset could not be found")
}
}
else if (length(xlab)==length(xvar.subset))
{
warning("The values of 'xlab' are taken to be corresponding to 'xvar.subset'")
lab <- xlab
xlab <- rep("",times=pExpl)
xlab[xvar.subset] <- lab
# It is required that xlab has a name for every column of x,
# if xvar.subset is passed, xlab of the same length is accepted.
}
else xvar.subset <- 1:pExpl
}
else xvar.subset <- 1:pExpl
pExpl <- length(xvar.subset)
######### END edit xvar.subset, n.xvars etc #########
########### factor plot #################
is.all.factor <- TRUE
for(i in xvar.subset) {
if (!is.factor(expl.data[,i])) is.all.factor <- FALSE
}
if(is.all.factor) {
# If the indep. var are all factors, draw a factor plot.
m[[1]]<- NULL
m$... <- NULL
m$data <- NULL
m$all.labels <- NULL
m$x <- NULL
m$y <- NULL
m$xvar.select <- m$n.xvars <- m$xvar.subset <- m$var.subset <-
m$nvars <- NULL
m <- lapply(m, eval, parent.frame())
lterm <- length(attr(terms(mvabund.formula),"term.labels"))
if(is.null(dots$legendtitle) & ( lterm==1 | lterm == pExpl)) {
m$legend.title <- attr(terms(mvabund.formula),"term.labels")[xvar.subset]
} else if(!is.null(dots$title)){
if(length(dots$title)>1) dots$title <- dots$title[xvar.subset]
}
if(shift) message("Overlapping points were shifted along the y-axis to make them visible.")
cat("\n \tPIPING to 1st MVFACTOR PLOT \n")
do.call( "default.plotMvaFactor", c(list(x=mvabund.object.1, y=expl.data[,xvar.subset, drop=FALSE],transformation="no"), m, dots))
return(invisible())
}
########### END factor plot #################
# print a warning if the y-axis was wanted to be log-transformed - this should
# be done inside the formula
if(logWarn) warning("logarithmic axis is only possible on the x axis in 'plot.mvformula'")
########## whether to write the plot ##################################
if (write.plot!="show") {
if (write.plot=="eps" | write.plot=="postscript") {
postscript(paste(filename,".eps", sep="") )
} else if (write.plot=="pdf") {
pdf(paste(filename,".pdf", sep="") )
} else if (write.plot=="jpeg" ) {
jpeg(paste(filename,".jpeg", sep=""))
} else if (write.plot=="bmp" ){
bmp(paste(filename,".bmp", sep=""))
}
# Specify the window where to draw the plot.
dev.curr <- dev.cur()
on.exit( dev.off(which = dev.curr) )
}
###########
mvabund.formula <- mvabund.formulaUni
###########
if (miss.xlab) xlab <- attr(terms(mvabund.formula[[1]]),"term.labels")
else
{
llab <- length(attr(terms(mvabund.formula[[1]]),"term.labels"))
# xlab needs to include names for every column of x.
# not only the ones chosen with x.var.subset
if (length(xlab)!= llab & length(xlab)>1)
stop("'names.ind' must be a vector consisting of ",
as.character(llab), " characters")
}
############################################################################
# multiple plots are produced using the formula #
# according to R's defaults in the the univariate case #
############################################################################
if ((dim(expl.data)[1])!=N ) {stop ("You have not given independent data with appropriate dimension.")}
######### BEGIN establish row and column sizes #########
if (is.null(mfcol))
{
#DW changes, 17-8-15: define mfrow here now that n.xvars is defined, and make it length one if one variable to plot
if(length(n.xvars[xvar.subset])==1)
mfrow=n.vars
else
mfrow=c(min(5,n.vars), min(3, n.xvars[xvar.select]))
# end change
perwindow <- mfrow
mfr <- TRUE
} else {
perwindow <- mfcol
mfr <- FALSE
}
if(length(perwindow)==1) {
if(any(c(pExpl, n.vars)==1)){
columns <-ceiling(sqrt(perwindow))
rows <-columns-1
if (rows*columns<perwindow) rows<-columns
} else if (( (pExpl<7) & n.vars<5 & (pExpl*n.vars)<=perwindow ) | (pExpl>1 & perwindow%%pExpl == 0) ) {
columns <- pExpl
rows <- ceiling(perwindow/columns)
} else {
tim <- ceiling((pExpl*n.vars)/perwindow)
if (pExpl > tim){
rows <- ceiling(pExpl/tim)
columns <- ceiling(perwindow/rows)
} else {
columns <- ceiling(n.vars/tim)
rows <- ceiling(perwindow/columns)
}
}
} else {
rows <- perwindow[1]
columns <- perwindow[2]
perwindow <- rows*columns
}
if(any(c(pExpl, n.vars)==1)) {
if (n.vars==1) all.labels <- TRUE
if(missing(ask)) {
if (perwindow < pExpl*n.vars) { ask <- TRUE }
else ask <- FALSE
}
} else {
if(missing(ask)) {
if (n.vars>columns | pExpl>rows) { ask <- TRUE }
else ask <- FALSE
}
}
opar <- par( "mfrow", "mfcol", "mar", "ask")
if((keep.window | all(opar$mfrow==c(rows,columns)))& write.plot=="show")
opar$mfrow <- opar$mfcol<- NULL
if(write.plot=="show") on.exit(par(opar))
if(write.plot=="show" & is.null(dev)) {
if (columns > (rows+1)){
width <- 16
height <- max(rows*width/columns,5)
} else {
height <- 11
width <- max(height*columns/rows,4)
}
dev.off()
do.call(dev.name, args=list(height=height,width=width))
}
######### END establish row and column sizes #########
if (all(par("mar")== c(5.1, 4.1, 4.1, 2.1)) & (rows*columns>1) ) {
if (overall.main=="") side3 <- 1.7 else side3 <- 1
# par(mar=c(3,4,side3,1)+0.1)
par(mar=c(1.5,1.75,side3,1)+0.1)
}
#Set outer margin of the image (in text lines), and margin line on which labels are drawn.
#DW change, 4/5/15
#par(oma=c(1,1,2,1))
par(oma=c(2,2,2,1),mgp=c(1.75,0.5,0))
######### BEGIN plot #########
t.lab <- substr(t.lab,1,1)
if(!t.lab %in% c("o", "t")) stop("You have passed an invalid 't.lab'")
if (t.lab=="o") {
# Obtain the transformation-axis ticklabels and tickmark values.
# at the moment transAxis only works with the original name (so 'x' must be used)
# and not if there are indices used, then the normal axis is drawn
trAxis <- try(transAxis(x))
if(class(trAxis)=="try-error"){
axes <- TRUE
yaxislabs <- yaxvalues <- NULL
} else {
if (trAxis$trans ) {
axes <- FALSE
opp <- par("mgp")
par(mgp=c(3,0.2,0))
if(write.plot=="show") on.exit(par(opp), add =TRUE)
} else axes <- TRUE
yaxislabs <- trAxis$yaxlabs[var.subset]
yaxvalues <- trAxis$yaxvalues[var.subset]
}
} else if (t.lab=="t") {
axes <- TRUE
yaxislabs <- yaxvalues <- NULL
} else stop("undefined value in 't.lab'")
if (write.plot=="show") { aske <- ask }
else { aske <- FALSE }
if (!mfr) {
rowsmfr <- rows
colsmfr <- rows <- columns
columns <- rowsmfr
}
if(any(c(pExpl, n.vars)==1)) {
layoutmat <- 1:min(pExpl*n.vars, perwindow)
} else {
layoutmat <- (1:(min(n.vars,rows)*min(pExpl,columns)))
layouthelp <- layouthelp2 <-matrix(0, ncol=columns,nrow=rows)
}
if(any(c(pExpl, n.vars)==1)) {
layouthelp <- matrix( c(layoutmat,
rep.int(0,times=max(c(0,rows*columns-max(layoutmat)))) ), ncol=columns,nrow=rows, byrow=TRUE )
uselayout2 <- FALSE
} else {
layouthelp[1:min(n.vars,rows),1:min(pExpl,columns)] <-
matrix(layoutmat, nrow=min(n.vars,rows),ncol=min(pExpl,columns), byrow=TRUE )
if(pExpl>columns & (pExpl%%columns)!=0) {
layouthelp2[(1:min(n.vars,rows)),1:(pExpl%%columns)] <-
matrix((1:(min(n.vars,rows)*(pExpl%%columns))),
nrow=min(n.vars,rows),ncol=(pExpl%%columns), byrow=TRUE )
uselayout2 <- TRUE
} else uselayout2 <- FALSE
}
if (!mfr) {
layouthelp <- t(layouthelp)
if (uselayout2) layouthelp2 <- t(layouthelp2)
}
if( write.plot!= "show") {
dev.set(which = dev.curr)
# Specify the window where to draw the plot.
}
if (! all(is.null(c(mfcol, mfrow))) | is.null(dev) | write.plot!="show") {
layout(layouthelp)
nonewwind <- FALSE
curr.layout <- layouthelp
} else {
nonewwind <- all.labels <- TRUE
}
winxlevel <- ceiling(pExpl/columns)
winylevel <- ceiling(n.vars/rows)
if(pExpl==1) winylevel <-1
if(n.vars==1) winxlevel <-1
displab1x <- (!all.labels & pExpl==1)
if(any(scale.lab == "ss")){
mxAll <- sapply( 1:n.vars, function(yi) max(mvabund.object.1[,yi],na.rm=TRUE))
wh.max <- which( abs( mxAll - max(mxAll, na.rm=TRUE)) < 1e-06 )[1]
mnAll <- sapply( 1:n.vars, function(yi)min(mvabund.object.1[,yi],na.rm=TRUE))
wh.min <- which( abs( mnAll - min(mnAll, na.rm=TRUE)) < 1e-06 )[1]
ylimleft <- mnAll[wh.min]-mnAll[wh.min]/10
ylimright <- mxAll[wh.max]+mxAll[wh.max]/10
if(is.null(yaxvalues[[wh.max]])) {
if(write.plot=="show" & is.null(dev)) {
do.call(dev.name, args=list(height=height,width=width))
} else {
do.call(dev.name, args=list())
}
# dev.set(which = dev.cur())
cat("\n \tPIPING TO 1st PLOT FORMULA FEATURE \n")
#DW changes, 4/5/15:making this initial plot invisible. Still shows as a blank page though!? Need to remove...
# do.call( "plotFormulafeature", c(list(mvabund.formula[[wh.max]],data=data, axes=axes, type=type[1],ylim=c(0,mxAll[wh.max]+mxAll[wh.max]/10), xvar.subset=1,ask=FALSE), dots ))
do.call( "plotFormulafeature", c(list(mvabund.formula[[wh.max]],data=data, axes=axes, type="n",xaxt="n",yaxt="n",bty="n",ylim=c(ylimleft,ylimright), xvar.subset=1), xlab="",ylab="", dots ))
# mtext(overall.main, outer = TRUE, cex = 0.9*par("cex.main"), col=par("col.main"), line = 1)
yaxTic <- axTicks(2)
yaxLab <- as.character(yaxTic)
# dev.off()
} else {
yaxTic <- yaxvalues[[wh.max]]
yaxLab <- yaxislabs[[wh.max]]
}
yaxvalues <- yaxislabs <- vector("list", length=n.vars)
for(yi in 1:n.vars) {
yaxvalues[[yi]] <- yaxTic
yaxislabs[[yi]] <- yaxLab
}
}
par(ask=aske)
# if( shift & type != "bx" )
# message("Overlapping points were shifted along the y-axis to make them visible.")
for (yi in 1:winylevel) {
if(any(c(pExpl, n.vars)==1) | nonewwind) {
rowsj <- 1:n.vars
} else {
rowsj <- (1:rows) + rows*(yi-1)
rowsj <- rowsj[rowsj<=n.vars]
}
for (xi in 1:winxlevel) {
if (!nonewwind) {
if (xi==winxlevel & uselayout2) {
layout(layouthelp2)
curr.layout <- layouthelp2
} else {
layout(layouthelp)
curr.layout <- layouthelp
}
}
if(any(c(pExpl, n.vars)==1) | nonewwind) {
colsj <- xvar.subset
} else {
colsj <- (1:columns) + columns*(xi-1)
colsj <- colsj[colsj<=pExpl]
colsj <- xvar.subset[colsj]
}
# Get max value for y axis.
max.y <- max(max(mvabund.object.1,na.rm=TRUE))
min.y <- min(min(mvabund.object.1,na.rm=TRUE))
for(yj in rowsj) {
if(mfr) {
if( all.labels | yj==rowsj[length(rowsj)] | (pExpl==1 & yj>n.vars-columns)) {
xlabyj <- xlab # [colsj]
} else {
xlabyj <- ""
}
if(displab1x & (yj%%columns !=1)) ylabj <- ""
else ylabj <- ylab[yj]
} else {
xlabyj <- xlab # [colsj]
if( all.labels | yj==rowsj[1] | (pExpl==1 & yj>n.vars-columns)) {
ylabj <- ylab[yj]
} else ylabj <- ""
}
if(any(scale.lab == "ss")) {
# ylimVal <- c(0,mxAll[wh.max]+mxAll[wh.max]/10)
ylimVal <- c(min.y, max.y)
} else {
# mx <- max(mvabund.object.1[,yj],na.rm=TRUE)
# ylimVal <- c(0,mx+mx/10)
ylimVal <- c(min.y, max.y)
}
cat("\n \tPIPING TO 2nd PLOT FORMULA FEATURE \n")
do.call( "plotFormulafeature", c(list(mvabund.formula[[yj]], data=data, ylab=ylabj, xlab=xlabyj, ask=aske, yaxis.ticks=yaxvalues[[yj]], yaxis.labs=yaxislabs[[yj]], fg=fg, cex.xaxis=par("cex.axis"), cex.yaxis=par("cex.axis"), axes=axes , col=col, pch=pch, type=type, las=las, main= main[yj], ylim=ylimVal, border=border, xvar.subset=colsj, cex.lab=par("cex.lab"),all.labels= if(mfr) all.labels else TRUE), dots ))
if (yj==rowsj[1] & overall.main!="" & par("oma")[3]>=1) {
mtext(overall.main, outer = TRUE, cex = 0.9*par("cex.main"), col=par("col.main"), line = 1)
}
}
}
}
if(xvar.select == TRUE & pExpl < NCOL(expl.data) & miss.xvarsubset ) {
# if(miss.xvarsubset)
tmp <- " \n(the variables with the highest average R^2)" # else {
# tmp <- " (user selected)" }
tmp2 <- paste(xlab[xvar.subset], collapse = ", ")
message("Only the independent variables ", tmp2 , " were included in the plot(s)", tmp, "." )
}
if(n.vars < p & miss.varsubset) {
# if(miss.varsubset)
tmp <- " \n(the variables with highest total abundance)" # else {
# tmp <- " (user selected)" }
tmp2 <- paste(colnames(mvabund.object.1), collapse = ", ")
message("Only the response variables ", tmp2 , " were included in the plot(s)", tmp, ".")
}
if(!is.null(subset.char)) message("Only the subset ", subset.char,
" of the cases was included in the plot(s) (user selected).")
pmfg <- par("mfg")
######### END plot #########
if(keep.window & !(pmfg[1]==rows & pmfg[2]==columns)) {
if( (rows-pmfg[1]+columns-pmfg[2])==1) {
# ie if only one free plot left.
par(mfrow=c(rows,columns),mfg=c(rows,columns))
return(invisible())
} else if(columns==1) {
par(mfrow=c(rows,columns),mfg=c(pmfg[1]+1,1))
return(invisible())
} else if(rows==1) {
par(mfrow=c(rows,columns),mfg=c(1, pmfg[2]+1))
return(invisible())
}
if(max(curr.layout)== curr.layout[pmfg[1], pmfg[2]]) {
# ie he current.layout is completely used.
if(curr.layout[2,1]-curr.layout[1,1]==1) {
wh.next <- which(curr.layout==0)
mfgrows <- wh.next[1]%%rows
if(mfgrows==0) mfgrows <- rows
par(mfrow=c(rows,columns),mfg=c(mfgrows,ceiling(wh.next[1]/rows)))
} else {
wh.next <- which(curr.layout==0)
mfgrows <- wh.next[1]%%rows
if(mfgrows==0) mfgrows <- rows
par(mfcol=c(rows,columns),mfg=c(mfgrows,ceiling(wh.next[1]/rows)))
}
} else {
# ie the curr layout has still at least one free place
plot(0, axes=FALSE, type="n", xlab="", ylab="")
# Check which is the next free place.
pmfg2 <- par("mfg")
if(pmfg2[1]>pmfg[1] | ( pmfg2[1]!= pmfg[1] & pmfg2[1]==1 ) ) {
par(mfrow=c(rows,columns),mfg=c(pmfg2[1],pmfg2[2]))
} else par(mfcol=c(rows,columns),mfg=c(pmfg2[1],pmfg2[2]))
}
}
}
Try the mvabund package in your browser
Any scripts or data that you put into this service are public.
mvabund documentation built on March 18, 2022, 7:25 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions default.plot.mvformula
################################################################################
## PLOT.mvformula: Plot functions for mvabund objects #
## (multivariate abundance data) #
################################################################################
default.plot.mvformula <- function( x,
y,
type="p",
main,
xlab,
ylab="Abundances",
col=if(type=="bx") "white" else "black",
fg="grey",
pch=1,
las=1,
write.plot="show",
filename="plot.mvabund",
n.vars=if(any(is.na(var.subset))) 12 else length(var.subset),
overall.main="",
data=parent.frame(),
var.subset=NA,
xvar.select=TRUE,
n.xvars=NA,
# n.xvars=if(any(is.na(xvar.subset))) min(3, sum(!is.interaction)) else length(xvar.subset),
xvar.subset=NA,
scale.lab="ss",
t.lab="t",
mfrow=NULL,
mfcol=NULL,
shift=TRUE,
border="black",
all.labels=FALSE,
keep.window=FALSE,
ask, ... ) {
#Kill Function if a BoxPlot has been passed
if (type=="bx") { stop("\nERROR: It is not ideal to use a boxplot for this data type\n") }
# If there is only one response variable this is changed to TRUE
# (overwritten independent of passed value)
dev <- dev.list()
dev.name <- getOption("device")
if(is.null(dev.name))
stop("Make sure that the 'device' option has a valid value, e.g. 'options(device = 'windows')'. Allowed values here are 'windows', 'win.graph', 'x11', 'X11'.")
# if(!(any(dev.name == c("windows", "win.graph", "x11", "X11")) ) )
# stop("Make sure that the 'device' option has a valid value, e.g. 'options(device = 'windows')'.
# Allowed values here are 'windows', 'win.graph', 'x11', 'X11'.")
m <- match.call(expand.dots = FALSE)
dots <- m$...
logWarn <- FALSE
if(!is.null(dots$log)){
# dots$log <- NULL
if(regexpr("y", dots$log ) !=-1) {
dots$log <- sub("y","", dots$log)
logWarn <- TRUE
}
}
usex <- TRUE
if (length(dots)>0){
# Delete arguments in ... that are defined lateron and cannot be used twice in
# the plot function.
deactive <- c("axes", "cex.axis")
deactivate <- (1:length(dots))[names(dots) %in% deactive ]
for (i in length(deactivate):1) {
dots[ deactivate[i] ] <- NULL #fixed the [[]] compile problem, could reduace functionallity.
}
dots <- lapply(dots, eval, parent.frame())
if(!is.null(dots$formula)) {
x <- dots$formula
dots$formula <- NULL
usex <- FALSE
}
subset.char <- dots$subset
} else subset.char <- NULL
miss.varsubset <- missing(var.subset)
# Change logical var.subset to numerical var.subset, if necessary. Note that
# NA values are logical as well, but should be excluded here.
if(!miss.varsubset){
if(is.logical(var.subset) & any(!is.na(var.subset)))
var.subset <- which(var.subset[!is.na(var.subset)])
}
miss.varsubset <- !is.numeric(var.subset) # If this function is called within
# another, the missing function could be tricked out.
var.subset <- as.vector(var.subset)
miss.xvarsubset <- missing(xvar.subset)
if(!miss.xvarsubset){
if(is.logical(xvar.subset) & any(!is.na(xvar.subset)))
var.subset <- which(xvar.subset[!is.na(xvar.subset)])
}
miss.xvarsubset <- any(is.na(xvar.subset)) || any(is.null(xvar.subset))
if ( (missing(x) & usex) ) stop("'formula' is missing or incorrect")
if(!inherits(x, "formula") ) {
# if no formula is passed a formula of supposedly passed
# matrices/vectors/data.frames will be tried to build.
# If there is only one mvabund object, it is assumed to be the response.
# If there are two mvabund objects, x (the first argument passed) is the
# response, !note that this is different in plot.mvabund
# The response cannot be a data.frame, when both objects are data.frames an
# error will be produced.
# If none is an mvabund object, x is the response, if it is not a data.frame.
if (missing(y)) { stop("A response and independent variables need to be provided for a formula.") }
if( (is.mvabund(y) & ! is.mvabund(x) ) | is.data.frame(x)) {
if(is.data.frame(y)) stop("a data.frame is not accepted as response for building the formula")
targs <- match.call(call = sys.call(which = 1), expand.dots = FALSE)
mvabund.formula <- x <- as.formula(paste(deparse(targs$y,width.cutoff = 500),"~",deparse(targs$x,width.cutoff = 500), sep=""))
# That might result in troubles if either targs$y or targs$x has the name "x"
# try to remove x (formula) earlier after giving it the name mvabund.formula
} else {
targs <- match.call(call = sys.call(which = 1), expand.dots = FALSE)
mvabund.formula <- x <- as.formula( paste(deparse(targs$x,width.cutoff = 500),"~", deparse(targs$y,width.cutoff = 500), sep="") )
}
# if(length(x) != 3) stop("'formula' is incorrect")
} else mvabund.formula <- x
resp <- attr(terms(mvabund.formula),"response")
term.labels <- attr(terms(mvabund.formula),"term.labels")
if(resp) {
if(length(mvabund.formula) != 3) stop("'formula' is incorrect")
mvabund.object.1<-as.matrix(model.response(model.frame( mvabund.formula, data=data)))
nam <- paste(deparse( mvabund.formula[[2]]), ".", sep="")
if(!is.null(colnames( mvabund.object.1 ))){
if(any(is.na(suppressWarnings(as.numeric(sub( nam ,"", colnames( mvabund.object.1 ))))))){
colnames( mvabund.object.1 ) <- sub( nam ,"", colnames( mvabund.object.1 ))
}
}
} else stop("'formula' has no response")
mvabund.object.1 <- mvabund(mvabund.object.1, neg=TRUE)
N <- NROW(mvabund.object.1) # number of sites
p <- NCOL(mvabund.object.1) # number of organism types
# Automatically correct n.vars if > p
n.vars <- min(n.vars,p)
if(length(term.labels)!=0){
foo.ext <- extend.x.formula(mvabund.formula,return.interaction=TRUE, extend.term=TRUE)
mvabund.formula <- foo.ext$formula
is.interaction <- foo.ext$is.interaction
# the length of this shows the number of x variables, including interaction
# term and factors, but factors only appear once and not like in model.matrix
# with (one column for every level) - 1
} else {
is.interaction <- c()
}
# Exclude Interactions.
if(any(is.interaction)){
mvabund.formula <- terms(mvabund.formula)[-which(is.interaction)]
print(2)
# terms.foo <- attr(model.frame(mvabund.formula), "terms")
# facs <- attr(terms.foo,"factors")[,is.interaction, drop = FALSE]
# datClass <- attr(terms.foo,"dataClasses")
# for(faci in 1:NCOL(facs)){
# if(any(datClass[facs[,faci]>0] == "factor"))
# stop("Formulas with interactions that include factors cannot yet be plotted")
# }
}
mf <- model.frame(mvabund.formula, data=data)
expl.data <- as.data.frame(mf[, -1])
# This is a matrix, ie factors are split to dummy variables
# not used for the actual plot, but for getting assign, etc
pExpl <- length(is.interaction) # NCOL(expl.data) + sum(is.interaction)
if(pExpl == 0) stop("The model is empty. There's nothing to plot.")
miss.xlab <- missing(xlab)
if (miss.xlab) xlab <- NULL
miss.main <- missing(main)
if (miss.main) {
nam <- paste(deparse( mvabund.formula[[2]]), ".", sep="")
# main <- colnames(as.matrix( model.frame( mvabund.formula)))[1:p]
main <- colnames(mvabund.object.1)
if(!is.null(main)){
if(any(is.na(suppressWarnings(as.numeric( sub(nam,"", main) ))))){
main <- colnames(mvabund.object.1) <- sub(nam,"", main)
} else colnames(mvabund.object.1) <- main
}
# main <- colnames(mvabund.object.1) <- colnames(as.matrix(
# model.frame( mvabund.formula)))[1:p]
n.mn <- nchar(main, type="chars")
if(any(n.mn > 14)) main[n.mn > 14] <- paste(substr(main[n.mn > 14], 1, 7),"...\n...",
substr(main[n.mn > 14], n.mn[n.mn > 14]-7, n.mn[n.mn > 14]) )
} else if (length(main)==p){
colnames(mvabund.object.1) <- main
} else if (!all(is.na(var.subset)) & length(main)==length(var.subset)){
colnames(mvabund.object.1)[var.subset]<- main
main <- colnames(mvabund.object.1)
} else if (length(main)==1) {
main <- rep(main,times =p)
} else stop("'main' must be a vector consisting of one or ", as.character(p), " characters" )
if (length(main)==1) { main <- rep(main,times= p) }
######### BEGIN edit var.subset, n.vars and mvabund.objects #########
# subset allows double variables
var.subset <- as.vector(var.subset)
var.subset.dim <- length(var.subset)
if (miss.varsubset) {
if(!any(length(ylab)==c(1,p))) stop("the length of 'ylab' is not correct")
sum.mvabund.object.1 <-
t(mvabund.object.1) %*% matrix(1,ncol=1,nrow=N)
# Find abundance ranks OF MVABUND.OBJECT.1.
var.subset<-order(sum.mvabund.object.1, decreasing = TRUE)
# Ensure no more than n.vars in var.subset.
if (n.vars < length(var.subset)) var.subset <- var.subset[1:n.vars]
var.subset.dim <- length(var.subset)
} else if ( p<max(var.subset) ) {
# Some dimension checks for subset, n.vars.
stop ("You have given an invalid 'var.subset'.")
} else if (n.vars!=var.subset.dim) {
n.vars <- var.subset.dim
}
mvabund.object.1 <- mvabund.object.1[,var.subset, drop=FALSE]
if(length(ylab)==p) {
ylab <- ylab[var.subset]
} else if (length(ylab)==1) ylab <- rep(ylab, times=n.vars)
main<- main[var.subset]
######### END edit var.subset, n.vars and mvabund.objects #########
mvabund.formulaUni <- formulaUnimva( mvabund.formula, var.subset=var.subset, split.x=TRUE, intercept=0 )
######### BEGIN edit xvar.subset, n.xvars etc #########
#DW changes, 30/10/14. 17/8: moved this outside of variable selection loop
#only define default n.xvars here, once is.interaction has been defined
if(is.na(n.xvars))
{
if(any(is.na(xvar.subset)))
n.xvars = min(3, sum(!is.interaction))
else
n.xvars = length(xvar.subset)
}
if(xvar.select & pExpl>1) #begin varaible selection (if required)
{
if(is.list(pch))
{
pch2 <- list()
for(i in 1:length(xvar.subset))
{
pch2[[xvar.subset[i] ]] <- pch[[i]]
}
pch <- pch2
}
if(is.list(col))
{
col2 <- list()
for(i in 1:length(xvar.subset))
{
col2[[xvar.subset[i]]] <- pch[[i]]
}
col<- col2
}
if(miss.xvarsubset)
{
# DW change, 17/8/15: changed to a warning if no n.xvars passed.
if(any(names(dots)==n.xvars)==FALSE)
warning(paste("No sQuote(n.xvars) passed so subset selection will be for",n.xvars,"variables"))
# Find the n.xvars independent variables with the highest average R^2.
xvar.subset <- best.r.sq(formula= mvabund.formula, # mvabund.formulaUni,
var.subset=var.subset, n.xvars= min(n.xvars, pExpl))
xvar.subset=xvar.subset$xs #to keep the indices only
if (is.null(xvar.subset))
{
xvar.subset <- 1:pExpl
warning("xvar.subset could not be found")
}
}
else if (length(xlab)==length(xvar.subset))
{
warning("The values of 'xlab' are taken to be corresponding to 'xvar.subset'")
lab <- xlab
xlab <- rep("",times=pExpl)
xlab[xvar.subset] <- lab
# It is required that xlab has a name for every column of x,
# if xvar.subset is passed, xlab of the same length is accepted.
}
else xvar.subset <- 1:pExpl
}
else xvar.subset <- 1:pExpl
pExpl <- length(xvar.subset)
######### END edit xvar.subset, n.xvars etc #########
########### factor plot #################
is.all.factor <- TRUE
for(i in xvar.subset) {
if (!is.factor(expl.data[,i])) is.all.factor <- FALSE
}
if(is.all.factor) {
# If the indep. var are all factors, draw a factor plot.
m[[1]]<- NULL
m$... <- NULL
m$data <- NULL
m$all.labels <- NULL
m$x <- NULL
m$y <- NULL
m$xvar.select <- m$n.xvars <- m$xvar.subset <- m$var.subset <-
m$nvars <- NULL
m <- lapply(m, eval, parent.frame())
lterm <- length(attr(terms(mvabund.formula),"term.labels"))
if(is.null(dots$legendtitle) & ( lterm==1 | lterm == pExpl)) {
m$legend.title <- attr(terms(mvabund.formula),"term.labels")[xvar.subset]
} else if(!is.null(dots$title)){
if(length(dots$title)>1) dots$title <- dots$title[xvar.subset]
}
if(shift) message("Overlapping points were shifted along the y-axis to make them visible.")
cat("\n \tPIPING to 1st MVFACTOR PLOT \n")
do.call( "default.plotMvaFactor", c(list(x=mvabund.object.1, y=expl.data[,xvar.subset, drop=FALSE],transformation="no"), m, dots))
return(invisible())
}
########### END factor plot #################
# print a warning if the y-axis was wanted to be log-transformed - this should
# be done inside the formula
if(logWarn) warning("logarithmic axis is only possible on the x axis in 'plot.mvformula'")
########## whether to write the plot ##################################
if (write.plot!="show") {
if (write.plot=="eps" | write.plot=="postscript") {
postscript(paste(filename,".eps", sep="") )
} else if (write.plot=="pdf") {
pdf(paste(filename,".pdf", sep="") )
} else if (write.plot=="jpeg" ) {
jpeg(paste(filename,".jpeg", sep=""))
} else if (write.plot=="bmp" ){
bmp(paste(filename,".bmp", sep=""))
}
# Specify the window where to draw the plot.
dev.curr <- dev.cur()
on.exit( dev.off(which = dev.curr) )
}
###########
mvabund.formula <- mvabund.formulaUni
###########
if (miss.xlab) xlab <- attr(terms(mvabund.formula[[1]]),"term.labels")
else
{
llab <- length(attr(terms(mvabund.formula[[1]]),"term.labels"))
# xlab needs to include names for every column of x.
# not only the ones chosen with x.var.subset
if (length(xlab)!= llab & length(xlab)>1)
stop("'names.ind' must be a vector consisting of ",
as.character(llab), " characters")
}
############################################################################
# multiple plots are produced using the formula #
# according to R's defaults in the the univariate case #
############################################################################
if ((dim(expl.data)[1])!=N ) {stop ("You have not given independent data with appropriate dimension.")}
######### BEGIN establish row and column sizes #########
if (is.null(mfcol))
{
#DW changes, 17-8-15: define mfrow here now that n.xvars is defined, and make it length one if one variable to plot
if(length(n.xvars[xvar.subset])==1)
mfrow=n.vars
else
mfrow=c(min(5,n.vars), min(3, n.xvars[xvar.select]))
# end change
perwindow <- mfrow
mfr <- TRUE
} else {
perwindow <- mfcol
mfr <- FALSE
}
if(length(perwindow)==1) {
if(any(c(pExpl, n.vars)==1)){
columns <-ceiling(sqrt(perwindow))
rows <-columns-1
if (rows*columns<perwindow) rows<-columns
} else if (( (pExpl<7) & n.vars<5 & (pExpl*n.vars)<=perwindow ) | (pExpl>1 & perwindow%%pExpl == 0) ) {
columns <- pExpl
rows <- ceiling(perwindow/columns)
} else {
tim <- ceiling((pExpl*n.vars)/perwindow)
if (pExpl > tim){
rows <- ceiling(pExpl/tim)
columns <- ceiling(perwindow/rows)
} else {
columns <- ceiling(n.vars/tim)
rows <- ceiling(perwindow/columns)
}
}
} else {
rows <- perwindow[1]
columns <- perwindow[2]
perwindow <- rows*columns
}
if(any(c(pExpl, n.vars)==1)) {
if (n.vars==1) all.labels <- TRUE
if(missing(ask)) {
if (perwindow < pExpl*n.vars) { ask <- TRUE }
else ask <- FALSE
}
} else {
if(missing(ask)) {
if (n.vars>columns | pExpl>rows) { ask <- TRUE }
else ask <- FALSE
}
}
opar <- par( "mfrow", "mfcol", "mar", "ask")
if((keep.window | all(opar$mfrow==c(rows,columns)))& write.plot=="show")
opar$mfrow <- opar$mfcol<- NULL
if(write.plot=="show") on.exit(par(opar))
if(write.plot=="show" & is.null(dev)) {
if (columns > (rows+1)){
width <- 16
height <- max(rows*width/columns,5)
} else {
height <- 11
width <- max(height*columns/rows,4)
}
dev.off()
do.call(dev.name, args=list(height=height,width=width))
}
######### END establish row and column sizes #########
if (all(par("mar")== c(5.1, 4.1, 4.1, 2.1)) & (rows*columns>1) ) {
if (overall.main=="") side3 <- 1.7 else side3 <- 1
# par(mar=c(3,4,side3,1)+0.1)
par(mar=c(1.5,1.75,side3,1)+0.1)
}
#Set outer margin of the image (in text lines), and margin line on which labels are drawn.
#DW change, 4/5/15
#par(oma=c(1,1,2,1))
par(oma=c(2,2,2,1),mgp=c(1.75,0.5,0))
######### BEGIN plot #########
t.lab <- substr(t.lab,1,1)
if(!t.lab %in% c("o", "t")) stop("You have passed an invalid 't.lab'")
if (t.lab=="o") {
# Obtain the transformation-axis ticklabels and tickmark values.
# at the moment transAxis only works with the original name (so 'x' must be used)
# and not if there are indices used, then the normal axis is drawn
trAxis <- try(transAxis(x))
if(class(trAxis)=="try-error"){
axes <- TRUE
yaxislabs <- yaxvalues <- NULL
} else {
if (trAxis$trans ) {
axes <- FALSE
opp <- par("mgp")
par(mgp=c(3,0.2,0))
if(write.plot=="show") on.exit(par(opp), add =TRUE)
} else axes <- TRUE
yaxislabs <- trAxis$yaxlabs[var.subset]
yaxvalues <- trAxis$yaxvalues[var.subset]
}
} else if (t.lab=="t") {
axes <- TRUE
yaxislabs <- yaxvalues <- NULL
} else stop("undefined value in 't.lab'")
if (write.plot=="show") { aske <- ask }
else { aske <- FALSE }
if (!mfr) {
rowsmfr <- rows
colsmfr <- rows <- columns
columns <- rowsmfr
}
if(any(c(pExpl, n.vars)==1)) {
layoutmat <- 1:min(pExpl*n.vars, perwindow)
} else {
layoutmat <- (1:(min(n.vars,rows)*min(pExpl,columns)))
layouthelp <- layouthelp2 <-matrix(0, ncol=columns,nrow=rows)
}
if(any(c(pExpl, n.vars)==1)) {
layouthelp <- matrix( c(layoutmat,
rep.int(0,times=max(c(0,rows*columns-max(layoutmat)))) ), ncol=columns,nrow=rows, byrow=TRUE )
uselayout2 <- FALSE
} else {
layouthelp[1:min(n.vars,rows),1:min(pExpl,columns)] <-
matrix(layoutmat, nrow=min(n.vars,rows),ncol=min(pExpl,columns), byrow=TRUE )
if(pExpl>columns & (pExpl%%columns)!=0) {
layouthelp2[(1:min(n.vars,rows)),1:(pExpl%%columns)] <-
matrix((1:(min(n.vars,rows)*(pExpl%%columns))),
nrow=min(n.vars,rows),ncol=(pExpl%%columns), byrow=TRUE )
uselayout2 <- TRUE
} else uselayout2 <- FALSE
}
if (!mfr) {
layouthelp <- t(layouthelp)
if (uselayout2) layouthelp2 <- t(layouthelp2)
}
if( write.plot!= "show") {
dev.set(which = dev.curr)
# Specify the window where to draw the plot.
}
if (! all(is.null(c(mfcol, mfrow))) | is.null(dev) | write.plot!="show") {
layout(layouthelp)
nonewwind <- FALSE
curr.layout <- layouthelp
} else {
nonewwind <- all.labels <- TRUE
}
winxlevel <- ceiling(pExpl/columns)
winylevel <- ceiling(n.vars/rows)
if(pExpl==1) winylevel <-1
if(n.vars==1) winxlevel <-1
displab1x <- (!all.labels & pExpl==1)
if(any(scale.lab == "ss")){
mxAll <- sapply( 1:n.vars, function(yi) max(mvabund.object.1[,yi],na.rm=TRUE))
wh.max <- which( abs( mxAll - max(mxAll, na.rm=TRUE)) < 1e-06 )[1]
mnAll <- sapply( 1:n.vars, function(yi)min(mvabund.object.1[,yi],na.rm=TRUE))
wh.min <- which( abs( mnAll - min(mnAll, na.rm=TRUE)) < 1e-06 )[1]
ylimleft <- mnAll[wh.min]-mnAll[wh.min]/10
ylimright <- mxAll[wh.max]+mxAll[wh.max]/10
if(is.null(yaxvalues[[wh.max]])) {
if(write.plot=="show" & is.null(dev)) {
do.call(dev.name, args=list(height=height,width=width))
} else {
do.call(dev.name, args=list())
}
# dev.set(which = dev.cur())
cat("\n \tPIPING TO 1st PLOT FORMULA FEATURE \n")
#DW changes, 4/5/15:making this initial plot invisible. Still shows as a blank page though!? Need to remove...
# do.call( "plotFormulafeature", c(list(mvabund.formula[[wh.max]],data=data, axes=axes, type=type[1],ylim=c(0,mxAll[wh.max]+mxAll[wh.max]/10), xvar.subset=1,ask=FALSE), dots ))
do.call( "plotFormulafeature", c(list(mvabund.formula[[wh.max]],data=data, axes=axes, type="n",xaxt="n",yaxt="n",bty="n",ylim=c(ylimleft,ylimright), xvar.subset=1), xlab="",ylab="", dots ))
# mtext(overall.main, outer = TRUE, cex = 0.9*par("cex.main"), col=par("col.main"), line = 1)
yaxTic <- axTicks(2)
yaxLab <- as.character(yaxTic)
# dev.off()
} else {
yaxTic <- yaxvalues[[wh.max]]
yaxLab <- yaxislabs[[wh.max]]
}
yaxvalues <- yaxislabs <- vector("list", length=n.vars)
for(yi in 1:n.vars) {
yaxvalues[[yi]] <- yaxTic
yaxislabs[[yi]] <- yaxLab
}
}
par(ask=aske)
# if( shift & type != "bx" )
# message("Overlapping points were shifted along the y-axis to make them visible.")
for (yi in 1:winylevel) {
if(any(c(pExpl, n.vars)==1) | nonewwind) {
rowsj <- 1:n.vars
} else {
rowsj <- (1:rows) + rows*(yi-1)
rowsj <- rowsj[rowsj<=n.vars]
}
for (xi in 1:winxlevel) {
if (!nonewwind) {
if (xi==winxlevel & uselayout2) {
layout(layouthelp2)
curr.layout <- layouthelp2
} else {
layout(layouthelp)
curr.layout <- layouthelp
}
}
if(any(c(pExpl, n.vars)==1) | nonewwind) {
colsj <- xvar.subset
} else {
colsj <- (1:columns) + columns*(xi-1)
colsj <- colsj[colsj<=pExpl]
colsj <- xvar.subset[colsj]
}
# Get max value for y axis.
max.y <- max(max(mvabund.object.1,na.rm=TRUE))
min.y <- min(min(mvabund.object.1,na.rm=TRUE))
for(yj in rowsj) {
if(mfr) {
if( all.labels | yj==rowsj[length(rowsj)] | (pExpl==1 & yj>n.vars-columns)) {
xlabyj <- xlab # [colsj]
} else {
xlabyj <- ""
}
if(displab1x & (yj%%columns !=1)) ylabj <- ""
else ylabj <- ylab[yj]
} else {
xlabyj <- xlab # [colsj]
if( all.labels | yj==rowsj[1] | (pExpl==1 & yj>n.vars-columns)) {
ylabj <- ylab[yj]
} else ylabj <- ""
}
if(any(scale.lab == "ss")) {
# ylimVal <- c(0,mxAll[wh.max]+mxAll[wh.max]/10)
ylimVal <- c(min.y, max.y)
} else {
# mx <- max(mvabund.object.1[,yj],na.rm=TRUE)
# ylimVal <- c(0,mx+mx/10)
ylimVal <- c(min.y, max.y)
}
cat("\n \tPIPING TO 2nd PLOT FORMULA FEATURE \n")
do.call( "plotFormulafeature", c(list(mvabund.formula[[yj]], data=data, ylab=ylabj, xlab=xlabyj, ask=aske, yaxis.ticks=yaxvalues[[yj]], yaxis.labs=yaxislabs[[yj]], fg=fg, cex.xaxis=par("cex.axis"), cex.yaxis=par("cex.axis"), axes=axes , col=col, pch=pch, type=type, las=las, main= main[yj], ylim=ylimVal, border=border, xvar.subset=colsj, cex.lab=par("cex.lab"),all.labels= if(mfr) all.labels else TRUE), dots ))
if (yj==rowsj[1] & overall.main!="" & par("oma")[3]>=1) {
mtext(overall.main, outer = TRUE, cex = 0.9*par("cex.main"), col=par("col.main"), line = 1)
}
}
}
}
if(xvar.select == TRUE & pExpl < NCOL(expl.data) & miss.xvarsubset ) {
# if(miss.xvarsubset)
tmp <- " \n(the variables with the highest average R^2)" # else {
# tmp <- " (user selected)" }
tmp2 <- paste(xlab[xvar.subset], collapse = ", ")
message("Only the independent variables ", tmp2 , " were included in the plot(s)", tmp, "." )
}
if(n.vars < p & miss.varsubset) {
# if(miss.varsubset)
tmp <- " \n(the variables with highest total abundance)" # else {
# tmp <- " (user selected)" }
tmp2 <- paste(colnames(mvabund.object.1), collapse = ", ")
message("Only the response variables ", tmp2 , " were included in the plot(s)", tmp, ".")
}
if(!is.null(subset.char)) message("Only the subset ", subset.char,
" of the cases was included in the plot(s) (user selected).")
pmfg <- par("mfg")
######### END plot #########
if(keep.window & !(pmfg[1]==rows & pmfg[2]==columns)) {
if( (rows-pmfg[1]+columns-pmfg[2])==1) {
# ie if only one free plot left.
par(mfrow=c(rows,columns),mfg=c(rows,columns))
return(invisible())
} else if(columns==1) {
par(mfrow=c(rows,columns),mfg=c(pmfg[1]+1,1))
return(invisible())
} else if(rows==1) {
par(mfrow=c(rows,columns),mfg=c(1, pmfg[2]+1))
return(invisible())
}
if(max(curr.layout)== curr.layout[pmfg[1], pmfg[2]]) {
# ie he current.layout is completely used.
if(curr.layout[2,1]-curr.layout[1,1]==1) {
wh.next <- which(curr.layout==0)
mfgrows <- wh.next[1]%%rows
if(mfgrows==0) mfgrows <- rows
par(mfrow=c(rows,columns),mfg=c(mfgrows,ceiling(wh.next[1]/rows)))
} else {
wh.next <- which(curr.layout==0)
mfgrows <- wh.next[1]%%rows
if(mfgrows==0) mfgrows <- rows
par(mfcol=c(rows,columns),mfg=c(mfgrows,ceiling(wh.next[1]/rows)))
}
} else {
# ie the curr layout has still at least one free place
plot(0, axes=FALSE, type="n", xlab="", ylab="")
# Check which is the next free place.
pmfg2 <- par("mfg")
if(pmfg2[1]>pmfg[1] | ( pmfg2[1]!= pmfg[1] & pmfg2[1]==1 ) ) {
par(mfrow=c(rows,columns),mfg=c(pmfg2[1],pmfg2[2]))
} else par(mfcol=c(rows,columns),mfg=c(pmfg2[1],pmfg2[2]))
}
}
}
Try the mvabund package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.