Nothing
if(getRversion()>="2.15.1") {
# remove warnings due to ggplot2 syntax strangeness
utils::globalVariables(c("comp","y","label","angle","hjust"))
}
#' SCGLR generic plot
#' @export
#' @importFrom stats cor aggregate
#' @importFrom grid gpar circleGrob grid.newpage pushViewport viewport grid.layout
#' @method plot SCGLR
#' @description SCGLR generic plot
#' @param x an object from SCGLR class.
#' @param \dots optional arguments (see \link{customize}).
#' @param style named list of values used to customize the plot (see \link{customize})
#' @param plane a size-2 vector (or comma separated string) indicating which components are plotted (eg: c(1,2) or "1,2").
#' @return an object of class \code{\link{ggplot}}.
#' @examples \dontrun{
#' library(SCGLR)
#'
#' # load sample data
#' data(genus)
#'
#' # get variable names from dataset
#' n <- names(genus)
#' ny <- n[grep("^gen",n)] # Y <- names that begins with "gen"
#' nx <- n[-grep("^gen",n)] # X <- remaining names
#'
#' # remove "geology" and "surface" from nx
#' # as surface is offset and we want to use geology as additional covariate
#' nx <-nx[!nx%in%c("geology","surface")]
#'
#' # build multivariate formula
#' # we also add "lat*lon" as computed covariate
#' form <- multivariateFormula(ny,c(nx,"I(lat*lon)"),c("geology"))
#'
#' # define family
#' fam <- rep("poisson",length(ny))
#'
#' genus.scglr <- scglr(formula=form,data = genus,family=fam, K=4,
#' offset=genus$surface)
#'
#' summary(genus.scglr)
#'
#' barplot(genus.scglr)
#'
#' plot(genus.scglr)
#'
#' plot(genus.scglr, predictors=TRUE, factor=TRUE)
#'
#' pairs(genus.scglr)
#'
#' }
plot.SCGLR <- function(x, ..., style=getOption("plot.SCGLR"), plane=c(1, 2)) {
data <- x
if (class(data) != "SCGLR")
stop("This plot function need an SCGLR result")
if (dim(data$compr)[2] < 2)
stop("At least two axes are needed for this kind of plot!")
# customize from remaining arguments
dots <- list(...)
# merged with style argument
if (!is.null(style)) {
style[names(dots)] <- dots
dots <- style
}
allowed_cust <- c(
"title", "expand", "labels.offset", "labels.auto", "labels.size", "threshold",
"observations",
"observations.factor", "observations.size", "observations.color",
"observations.alpha",
"predictors",
"predictors.color", "predictors.arrows", "predictors.arrows.color",
"predictors.labels", "predictors.labels.color", "predictors.labels.size",
"predictors.labels.auto",
"predictors.alpha", "predictors.arrows.alpha", "predictors.labels.alpha",
"covariates",
"covariates.color", "covariates.arrows", "covariates.arrows.color",
"covariates.labels", "covariates.labels.color", "covariates.labels.size",
"covariates.labels.auto",
"covariates.alpha", "covariates.arrows.alpha", "covariates.labels.alpha",
"factor",
"factor.points", "factor.points.size", "factor.points.shape",
"factor.labels", "factor.labels.color", "factor.labels.size"
)
match_cust <- function(key) {
# build regexp foo.bar --> ^foo[^.]*\\.bar[^.]*$
key2 <- paste("^", gsub("\\.", "[^.]*\\\\.", key), "[^.]*$", sep="")
# find matching cust key
key2 <- grep(key2, allowed_cust, value=TRUE)
# if not found (length==0) or ambiguity (length >1) revert to original
if(length(key2) != 1)
key2 <- key
key2
}
# resolve abbreviated customization key names
names(dots) <- sapply(names(dots),match_cust)
# get custom value from key (if not found return default value)
cust <- function(key, def=NULL) {
if(is.null(key)||is.null(dots[[key]])) {
return(def)
} else {
return(dots[[key]])
}
}
# check if key has a value or is not false
has_cust <- function(key, def=NULL) {
key <- cust(key, def)
!is.null(key) && (!is.logical(key) || key)
}
labels.offset <- cust("labels.offset",0.01)
labels.auto <- cust("labels.auto",TRUE)
labels.size <- cust("labels.size",1)
# process plane
if(is.character(plane)) {
plane <- as.integer(trim(unlist(strsplit(plane,","))))
}
# sanity checking
if(length(plane) != 2 ) {
stop("Plane should have two components!")
}
if((min(plane)<1) || (max(plane)>ncol(data$compr))) {
stop("Invalid components for plane!")
}
# plan
axis_names <- colnames(data$compr)[plane]
# check factor
factor <- cust("factor",NULL)
if(!is.null(factor) && (is.character(factor) || factor)) {
if(is.logical(factor)) {
if(is.null(data$xFactors))
stop("No factor in data!")
factor <- names(data$xFactors)[1]
warning("No factor given, assuming first one! (",factor,")!")
} else {
if(!factor %in% names(data$xFactors))
stop("Invalid factor!")
}
}
# inertia
inertia <- data$inertia[plane]
# layers
covariates <- has_cust("covariates",TRUE)
predictors <- has_cust("predictors",FALSE)
observations <- has_cust("observations",FALSE)
# build base plot
p <- qplot((-1:1)*cust("expand",1.0), (-1:1)*cust("expand",1.0), geom="blank")+
coord_fixed()+
# thicker x unit arrow
xlab(paste(axis_names[1],"(",round(100*inertia[1],2),"%)")) +
geom_hline(yintercept=0)+
geom_segment(aes(x=-1.1,xend=1.1,y=0,yend=0),size=1,arrow=arrow(length=unit(0.02,"npc")))+
# thicker y unit arrow
ylab(paste(axis_names[2],"(",round(100*inertia[2],2),"%)")) +
geom_vline(xintercept=0)+
geom_segment(aes(y=-1.1,yend=1.1,x=0,xend=0),size=1,arrow=arrow(length=unit(0.02,"npc")))
# plot title
if(has_cust("title", FALSE)) {
p <- p + ggtitle(cust("title"))
} else {
if(observations) {
if(covariates||predictors) {
p <- p + ggtitle("Mixed individual and \ncorrelation plot\n")
} else {
p <- p + ggtitle("Individual plot\n")
}
} else {
if(covariates||predictors) {
p <- p + ggtitle("Correlation plot\n")
} else {
p <- p + ggtitle("SCGLR Plot\n")
}
}
}
# add unit circle
if(covariates||predictors)
p <- p + annotation_custom(circleGrob(r=0.5,gp=gpar(fill=NA)),-1,1,-1,1)
# add threshold circle
if((covariates||predictors)&&has_cust("threshold"))
p <- p + annotation_custom(circleGrob(r=0.5*cust("threshold"),gp=gpar(fill=NA,lty="dashed")),-1,1,-1,1)
# add observations
if(observations) {
obs <- as.data.frame(data$compr[,plane])
names(obs) <- c("x","y")
# colored according to factor ?
if(!is.null(factor)&&(cust("observations.factor",FALSE))) {
obs <- cbind(obs,data$xFactors[factor])
p <- p + geom_point(
aes_string(x="x",y="y",color=factor),
data=obs,
size=cust("observations.size",1),
alpha=cust("observations.alpha",1)
)
} else {
p <- p + geom_point(
aes(x=x, y=y),
data=obs,
size=cust("observations.size",1),
color=cust("observations.color","black"),
alpha=cust("observations.alpha",1)
)
}
}
# add linear predictor arrows
if(predictors) {
predictors <- cust("predictors")
co <- as.data.frame(cor(data$lin.pred, data$compr[,plane]))
names(co) <- c("x", "y")
co$norm <- sqrt(co$x^2+co$y^2)
co$label <- rownames(co)
co$color <- cust("predictors.color","red")
co$alpha <- cust("predictors.alpha",1)
if(cust("predictors.arrows",TRUE)) {
co$arrows.color <- cust("predictors.arrows.color",co$color)
co$arrows.alpha <- cust("predictors.arrows.alpha",co$alpha)
}
if(cust("predictors.labels",TRUE)) {
co$labels.color <- cust("predictors.labels.color",co$color)
co$labels.alpha <- cust("predictors.labels.alpha",co$alpha)
co$labels.size <- cust("predictors.labels.size",4*labels.size)
}
# adjust label position
co$angle <- atan2(co$y,co$x)*180/pi
co$hjust <- ifelse(abs(co$angle)>90,1,0)
if(cust("predictors.labels.auto",labels.auto)) {
co$angle <- ifelse(abs(co$angle)>90,co$angle+180,co$angle)
} else {
co$angle <- 0
}
# filter according to user's will
if(is.character(predictors)) {
co <- co[co$label %in% predictors,]
}
# filter according to given threshold
if(cust("threshold",FALSE)>0)
co <- co[co$norm>cust("threshold"),]
if(nrow(co)==0) {
warning("No predictors with correlation higher than threshold value ",cust("threshold"),"!")
} else {
# draw arrows ?
if(cust("predictors.arrows",TRUE)) {
p <- p + geom_segment(
aes(x=0,y=0,xend=x,yend=y),
data=co,
color=co$arrows.color,
alpha=co$arrows.alpha,
arrow=arrow(length=unit(0.02,"npc"))
)
} else {
# reset label position
co$hjust <- 0.5
co$angle <- 0
}
# draw labels ?
if(cust("predictors.labels",TRUE)) {
p <- p + geom_text(
aes(x=x*(1+labels.offset/norm),y=y*(1+labels.offset/norm),label=label,angle=angle,hjust=hjust),
data=co,
color=co$labels.color,
alpha=co$labels.alpha,
size=co$labels.size
)
}
}
}
# add co-variate arrows
if(covariates) {
covariates <- cust("covariates")
co <- as.data.frame(cor(data$xNumeric, data$compr[,plane]))
names(co) <- c("x", "y")
co$norm <- sqrt(co$x^2+co$y^2)
co$label <- rownames(co)
co$color <- cust("covariates.color","black")
co$alpha <- cust("covariates.alpha",1)
if(cust("covariates.arrows",TRUE)) {
co$arrows.color <- cust("covariates.arrows.color",co$color)
co$arrows.alpha <- cust("covariates.arrows.alpha",co$alpha)
}
if(cust("covariates.labels",TRUE)) {
co$labels.color <- cust("covariates.labels.color",co$color)
co$labels.alpha <- cust("covariates.labels.alpha",co$alpha)
co$labels.size <- cust("covariates.labels.size",4*labels.size)
}
# adjust label position
co$angle <- atan2(co$y,co$x)*180/pi
co$hjust <- ifelse(abs(co$angle)>90,1,0)
if(cust("covariates.labels.auto",labels.auto)) {
co$angle <- ifelse(abs(co$angle)>90,co$angle+180,co$angle)
} else {
co$angle <- 0
}
# filter according to users's will
if(is.character(covariates)) {
co <- co[co$label %in% covariates,]
}
# filter according to given threshold
if(cust("threshold",FALSE)>0)
co <- co[co$norm>cust("threshold"),]
if(nrow(co)==0) {
warning("No correlation higher than threshold value ",cust("threshold"),"!")
} else {
# draw arrows ?
if(cust("covariates.arrows",TRUE)) {
p <- p + geom_segment(
aes(x=0,y=0,xend=x,yend=y),
data=co,
color=co$arrows.color,
alpha=co$arrows.alpha,
arrow=arrow(length=unit(0.02,"npc"))
)
} else {
# reset label adjust
co$hjust <- 0.5
co$angle <- 0
}
# draw labels ?
if(cust("covariates.labels",TRUE)) {
p <- p + geom_text(
aes(x=x*(1+labels.offset/norm),y=y*(1+labels.offset/norm),label=label,angle=angle,hjust=hjust),
data=co,
color=co$labels.color,
size=co$labels.size,
alpha=co$labels.alpha
)
}
}
}
# add factors
if(!is.null(factor)) {
bary <- aggregate(data$compr[,plane],data$xFactors[factor],mean)
names(bary) <- c(factor,"x","y")
# draw points ?
if(cust("factor.points",TRUE)) {
p <- p + geom_point(
aes_string(x="x",y="y",color=factor),
data=bary,
size=cust("factor.points.size",4),
shape=cust("factor.points.shape",13)
)
}
# draw labels ?
if(cust("factor.labels",TRUE)) {
p <- p + geom_text(
aes_string(x="x",y="y",label=factor),
data=bary,
color=cust("factor.labels.color","black"),
size=cust("factor.labels.size",4*labels.size)
)
}
}
# return plot
p
}
#' @title Barplot of percent of overall X variance captured by component
#' @export
#' @method barplot SCGLR
#' @description Barplot of percent of overall X variance captured by component
#' @param height object of class 'SCGLR', usually a result of running \code{\link{scglr}}.
#' @param \dots optional arguments.
#' @return an object of class ggplot.
#' @seealso For barplot application see examples in \code{\link{plot.SCGLR}}.
#' @keywords internal
barplot.SCGLR <- function(height, ...) {
.Deprecated("screeplot")
screeplot.SCGLR(height, ...)
}
#' @title Screeplot of percent of overall X variance captured by component
#' @export
#' @method screeplot SCGLR
#' @description Screeplot of percent of overall X variance captured by component
#' @param x object of class 'SCGLR', usually a result of running \code{\link{scglr}}.
#' @param \dots optional arguments.
#' @return an object of class ggplot.
#' @seealso For screeplot application see examples in \code{\link{plot.SCGLR}}.
screeplot.SCGLR <- function(x, ...) {
inertia <- data.frame(inertia=x$inertia,comp=1:length(x$inertia))
ggplot(inertia)+
geom_col(aes(x=comp,y=inertia))+
scale_x_discrete(labels=names(x$inertia),limits=1:length(inertia$comp))+
labs(x="Components", y="Inertia", title="Inertia per component\n", ...)
}
#' @title Pairwise scglr plot on components
#' @export
#' @method pairs SCGLR
#' @description Pairwise scglr plot on components
#' @param x object of class 'SCGLR', usually a result of running \code{\link{scglr}}.
#' @param \dots optionally, further arguments forwarded to \code{\link{plot.SCGLR}}.
#' @param nrow number of rows of the grid layout.
#' @param ncol number of columns of the grid layout.
#' @param components vector of integers selecting components to plot (default is all components).
# @return an object of class ggplot.
#' @seealso For pairs application see examples in \code{\link{plot.SCGLR}}
pairs.SCGLR <- function(x, ..., nrow=NULL, ncol=NULL, components=NULL) {
prm <- list(...)
nr <- nrow
nc <- ncol
prm["x"] <- list(x)
prm[["components"]] <- NULL
# pairs of components
if(is.null(components)) {
ncomp <- ncol(x$compr)
} else {
ncomp <- components
}
# sanity check
if((min(ncomp)<1) || (max(ncomp)>ncol(x$compr))) {
stop("Invalid components for plane!")
}
# build pairs
cmp_pairs <- combn(ncomp, 2, simplify=FALSE)
# build plot list
one_plot <- function(cmp_pair) {
do.call("plot.SCGLR", c(prm, plane=list(cmp_pair),title=paste(cmp_pair,collapse = "/")))
}
plots <- lapply(cmp_pairs, one_plot)
# arrange them in a grid
if(is.null(nc) & is.null(nr)) {
nr <- as.integer(sqrt(length(plots)))
}
# if(require("gridExtra",quietly = TRUE)) {
# do.call("arrangeGrob", c(plots, nrow=nr, ncol=nc,main="toto"))
# } else {
do.call("arrange", c(plots, nrow=nr, ncol=nc))
# }
}
## equivalent du par pour les ggplot2
vp.layout <- function(x, y) viewport(layout.pos.row=x, layout.pos.col=y)
arrange <- function(..., nrow=NULL, ncol=NULL, as.table=FALSE) {
dots <- list(...)
n <- length(dots)
if(is.null(nrow) & is.null(ncol)) { nrow = floor(n/2) ; ncol = ceiling(n/nrow)}
if(is.null(nrow)) { nrow = ceiling(n/ncol)}
if(is.null(ncol)) { ncol = ceiling(n/nrow)}
## NOTE see n2mfrow in grDevices for possible alternative
grid.newpage()
pushViewport(viewport(layout=grid.layout(nrow,ncol) ) )
ii.p <- 1
for(ii.row in seq(1, nrow)){
ii.table.row <- ii.row
if(as.table) {ii.table.row <- nrow - ii.table.row + 1}
for(ii.col in seq(1, ncol)){
ii.table <- ii.p
if(ii.p > n) break
print(dots[[ii.table]], vp=vp.layout(ii.table.row, ii.col))
ii.p <- ii.p + 1
}
}
invisible(NULL)
}
#' @title Barplot of percent of overall X variance captured by component
#' @export
#' @method barplot SCGLRTHM
#' @description Barplot of percent of overall X variance captured by component by theme
#' @param height object of class 'SCGLRTHM', usually a result of running \code{\link{theme}}.
#' @param \dots optional arguments.
#' @return an object of class ggplot.
#' @keywords internal
barplot.SCGLRTHM <- function(height, ...) {
.Deprecated("screeplot.SCGLRTHM")
screeplot.SCGLRTHM(height, ...)
}
#' @title Screeplot of percent of overall X variance captured by component
#' @export
#' @method screeplot SCGLRTHM
#' @description Screeplot of percent of overall X variance captured by component by theme
#' @param x object of class 'SCGLRTHM', usually a result of running \code{\link{theme}}.
#' @param \dots optional arguments.
#' @return an object of class ggplot.
screeplot.SCGLRTHM <- function(x, ...) {
plots <- lapply(x$themes, function(t)
screeplot.SCGLR(t,...)+
labs(title=paste0("Theme ",t$label,": Inertia per component\n"))
)
do.call(arrange, plots)
}
#' SCGLRTHM generic plot
#' @export
#' @method plot SCGLRTHM
#' @description SCGLR generic plot for themes
#' @param x an object from SCGLRTHM class.
#' @param \dots see SCGLR plot method
#' @return an object of class \code{\link{ggplot}}.
plot.SCGLRTHM <- function(x, ...) {
plots <- lapply(x$themes, function(t)
plot.SCGLR(t,...)+
labs(title=paste0("Theme ",t$label,": Correlation plot\n"))
)
do.call(arrange, plots)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.