R/utils.R

In VFP: Variance Function Program

# TODO: Add comment
# 
# Author: schueta6
###############################################################################

#' Determine Monotony of Vector.
#' 
#' @param x		(numeric) vector
#' 
#' @return 1 if monotonically increasing, -1 if monotonically decreasing, otherwise,
#'         a vector of 1 or -1 indicating increasing or decreasing status
#' 
#' @author Andre Schuetzenmeister \email{andre.schuetzenmeister@@roche.com}
#' 
#' @examples 
#' x1 <- seq(-1, 1, .1)
#' x2 <- seq(1, -1, -.1)
#' y <- x1^2 / 2 + x1/3 - 5
#' get_mon(x1)
#' get_mon(x2)
#' get_mon(y)

get_mon <- function(x) {
	if(all(x == cummax(x)))	{					# monotonically increasing
		return(1)
	} else if(all(-1 * x == cummax(-1 * x))) {	# monotonically decreasing
		return(-1) 
	} else {
		return(sign(diff(x)))
	}
}


#' Select the Best Fitting Model.
#' 
#' In case the user has not specified a specific model out of multiple fitted
#' models the one with the lowes AIC value will be selected. If multiple models
#' have the same AIC the one with the lowes complexity will be chosen. 
#' 
#' @param object		(object) of class "VFP"
#' @param model.no		(integer) specifying a fitted model stored in 'object'
#' @param cpx			(integer) vector specifying the order of complexity
#' 						used to select the less complex model when multiple, 
#'                      identical AIC values occur, taken from the fitted object
#' 						or c(1, 2, 3, 9, 5, 4, 7, 6, 8) is used if not specified
#' 
#' @return (integer) index of selected model in 'object$Models' or 'object$AIC'
#' 
#' @author Andre Schuetzenmeister \email{andre.schuetzenmeister@@roche.com}
#' 
#' @examples 
#' \donttest{
#' library(VCA)
#' data(VCAdata1)
#' lst <- anovaVCA(y~(device+lot)/day/run, VCAdata1, by="sample")
#' mat <- get_mat(lst)		# automatically selects "total"
#' res <- fit_vfp(model.no=1:9, Data=mat)
#' get_model(res)
#' }

get_model <- function(object, model.no = NULL, cpx=NULL) {
	
	stopifnot(is(object, "VFP"))
	stopifnot(model.no %in% 1:9)
	
	if(is.null(cpx)) {
		if(is.null(object$complexity)) {
			cpx <- c(1, 2, 3, 9, 5, 4, 7, 6, 8)
		} else {
			cpx <- object$complexity
		}
	}
	
	cpx <- paste0("Model_", cpx)
	
	if (is.null(model.no)) {				# automatically determine best fitting model
		AIC <- object$AIC
		nam <- names(AIC)
		if("Model_10" %in% nam)
			AIC <- AIC[-which(nam == "Model_10")]
		idx <- which(AIC == min(AIC, na.rm=TRUE))
		if(length(idx) > 1) {
			mods <- nam[idx]
			bmod <- cpx[which(cpx %in% mods)[1]]
			num  <- which(names(AIC) == bmod)
			attr(num, "model") <- bmod
		} else {
			num <- idx	
			attr(num, "model") <- nam[idx]
		}
		
		AIC <- NULL
	} else {
		models <- sub("Model_", "", names(object$RSS))
		if(!model.no %in% models)
			stop("Specified model ", paste0("'", model.no, "'"),
					" is not among fitted models: ", 
					paste(models, collapse=", "),"!")
		else {
			num <- which(models == model.no)
			attr(num, "model") <- paste0("Model_", model.no)
		}
	}
	num
}


#' Condition-Handling Without Losing Information.
#' 
#' Function is intented to wrap expressions provided and catching
#' all potentially useful information generated by the wrapped expression, i.e.
#' errors, warnings, and messages.
#' 
#' @param expr		(expression) for which exception handling should be provided
#' @param file		(character) string specifying a file to which all captured output
#' 					shall be written
#' 
#' @return 	(list) with element "result", "status" (0 = no warnings, no errors), 1 = warnings
#' 			were caught, 2 = errors were caught no result generated, "warnings", "errors", 
#' 			"messages"
#' 
#' @author Andre Schuetzenmeister \email{andre.schuetzenmeister@@roche.com}
#' 
#' @aliases conditionHandler
#' 
#' @examples 
#' condition_handler(warning("This is a warning!"))
#' f <- function(expr){warning("This a warning!"); eval(expr)}
#' condition_handler(f(1/2))
#' condition_handler(stop("This is an error!"))
#' condition_handler(1/"a")

condition_handler <- function(expr,  file=NULL) 
{
	Warnings <- Errors <- Messages <- NULL
	status   <- 0
	# handle warnings
	WHandler <- function(w) 
	{
		status <<- 1
		Warnings <<- c(Warnings, w$message)#list(w))
		invokeRestart("muffleWarning")		
	}
	# handle messages
	MHandler <- function(m)
	{
		Messages <<- c(Messages, m$message)#list(m))
		invokeRestart("muffleMessage")		
	}
	# handle errors
	if(!is.null(file) && file.exists(file))
	{
		capture.output(
				result <- try(withCallingHandlers(expr, warning = WHandler, message=MHandler), silent=TRUE),
				file=file, append=TRUE)
	}
	else
	{
		result <- try(withCallingHandlers(expr, warning = WHandler, message=MHandler), silent=TRUE)
	}
	
	if(is(result[[1]], "try-error") || is(result, "try-error"))
	{		
		Errors <- attr(result, "condition")$message
		result <- NA
		status <- 2
	}
	
	if(is.null(result))
	{
		status <- 2
		#	Errors <- "A model could not be fitted. Please check 'messages' and 'warnings'!"
		Errors <- paste("A model could not be fitted:", Warnings, Messages, sep=" | ")
	}
	
	res <- list(result = result, status=status, warnings = Warnings, errors = Errors, messages = Messages)
	
	res
} 


#' Transform list of VCA-object into VFP-matrix required for fitting.
#'
#' @param obj		(list) of VCA-objects
#' @param vc		(integer, character) either an integer specifying a variance component
#' 					or the name of a variance component; can also be a vector of integers
#' 					specifying a continuous sequence of variance components always including
#' 					'error' (repeatability)
#' 
#' @author Andre Schuetzenmeister \email{andre.schuetzenmeister@@roche.com}
#' 
#' @aliases getMat.VCA
#' 
#' @examples 
#' \donttest{
#' library(VCA)
#' data(VCAdata1)
#' lst <- anovaVCA(y~(device+lot)/day/run, VCAdata1, by="sample")
#' get_mat(lst)  # automatically selects 'total'
#' # pooled version of intermediate precision (error+run+day)
#' get_mat(lst, 4:6)
#' # only repeatability ('error')
#' get_mat(lst, "error")
#' # use remlVCA instead
#' lst2 <- remlVCA(y~(device+lot)/day/run, VCAdata1, by="sample")
#' get_mat(lst2)
#' }

get_mat <- function(obj, vc=1)
{
	stopifnot(is(obj, "list"))
	stopifnot(all(sapply(obj, class) == "VCA"))
	stopifnot(is.numeric(vc) || is.character(vc))
	tab <- obj[[1]]$aov.tab
	if(is.numeric(vc))
	{
		stopifnot(all(vc %in% 1:nrow(tab)))
		nm <- vc
		vc <- rownames(tab)[vc]
	}
	else
	{
		stopifnot(all(vc %in% rownames(tab)))
		nm <- (1:nrow(tab))[which(rownames(tab) %in% vc)]	
	}
	
	if(length(nm) > 1)
	{
		if(max(nm) != nrow(tab) || !(all(diff(nm) == 1)))
			stop("When specifying a sequence of variance components, it must include 'error' and be continuous!")
	}
	if(length(vc) > 1)		
	{
		if(obj[[1]]$EstMethod == "ANOVA") {
			tmp <- lapply(obj, stepwiseVCA)			# perform all combinations of VCA if ANOVA was used
			idx <- suppressWarnings(which(unlist(lapply(tmp[[1]], function(x) all(rownames(x$aov.tab) == c("total", vc))))))
			for(i in 1:length(tmp))
				tmp[[i]] <- tmp[[i]][[idx]]
		} else {									# objects generated by REML-estimation
			tmp <- lapply(obj, getIP.remlVCA, vc=vc[1])
		}
		
		obj <- tmp
		vc  <- "total"
	}
	
	mat 	<- t(sapply(obj, function(x) c(Mean=x$Mean, x$aov.tab[vc, c("DF", "VC")])))
	rn		<- rownames(mat)
	cn		<- colnames(mat)
	mat		<- matrix(as.numeric(mat), ncol=3, dimnames=list(rn, cn))
	mat		<- mat[order(unlist(mat[,"Mean"])),]
	mat		<- as.data.frame(mat)
	mat$SD 	<- sqrt(mat$VC)
	mat$CV 	<- mat$SD*100/mat$Mean
	mat
}



#' Adapted Version of Function 'signif'
#' 
#' This function adapts base-function \code{\link{signif}}
#' by always returning integer values in case the number of
#' requested significant digits is less than the the number of
#' digits in front of the decimal separator. In case of -Inf, Inf,
#' NA, and NaN the same value will be returned.
#' 
#' @param x			(numeric) value to be rounded to the desired number
#' 					of significant digits
#' @param digits	(integer) number of significant digits
#' @param as.char	(logical) TRUE = return character string, otherwise,
#' 					a numeric value will be returned
#' @param ...		additional parameters
#' 
#' @return 	number with 'digits' significant digits, if 'as.char=TRUE' "character" objects will be
#' 			returned otherwise objects of mode "numeric"
#' 
#' @aliases Signif
#' 
#' @author Andre Schuetzenmeister \email{andre.schuetzenmeister@@roche.com}
#' 
#' @examples 
#' identical(signif2(1.23456, 4), "1.235")
#' identical(signif2(-1.23456, 4), "-1.235")
#' identical(signif2(0.123456, 5), "0.12346")
#' identical(signif2(-12.5678, digits=2), "-13")
#' identical(signif2(0.490021, digits=4), "0.4900") 
#' identical(signif2(c(1.203, 4.56), digits=3, as.char=FALSE), c(1.20, 4.56))
#' identical(signif2(c(1.203, 4.56), digits=3, as.char=TRUE), c("1.20", "4.56"))
#' identical(signif2(0.003404, 3), "0.00340")
#' identical(signif2(0.00034004, 4), "0.0003400")
#' identical(signif2(12345.67, 4), "12346")
#' identical(signif2(12345.67, 3), "12346")
#' identical(signif2(123456, 3), "123456")

signif2 <- function(x, digits=4, as.char=TRUE, ...)
{

	call 	<- match.call()
	stopifnot(is.numeric(x))
	if(length(x) > 1)
		return(sapply(x, signif2, digits=digits, as.char=as.char))
	if(x %in% c(-Inf, NA, NaN, Inf))
		return(x)
	sgn		<- sign(x)
	x		<- abs(x)
	# number of digits left of comma
	Ndbc	<- nchar(substr(as.character(x), 1, regexpr("\\.", as.character(x))-1))
	if(Ndbc == 0) {
		Ndbc <- nchar(x)
	} else if(Ndbc==1 && substr(x, 1, 1) == "0") {
		Ndbc <- 0	
	}
	x 		<- signif(x, ifelse(Ndbc > digits, Ndbc, digits))
	NcX		<- nchar(x)
	comma	<- grepl("\\.", x)
	if(comma) {
		if(Ndbc > 0)
			NcX <- NcX - 1
		else { 							
			x2  <- sub("0.", "", x)
			Nz  <- regexpr("^[0]+", x2)
			if(Nz > -1){					# non-significant zeros right of comma
				NcX <- NcX - (2 + attr(Nz, "match.length"))
			} else {
				NcX <- NcX - 2
			}
		}
	}
	
	if(NcX < digits) {
		x <- paste0(x, ifelse(comma, "", "."), paste(rep(0, digits-NcX), collapse=""))
	}
	if(sgn == -1) {
		x <- paste0("-", x)
	}
	if(as.char) {
		x <- as.character(x)
	} else {
		x <- as.numeric(x)
	}
	x
}



#' Add a Grid to an Existing Plot.
#' 
#' It is possible to use automatically determined grid lines (\code{x=NULL, y=NULL}) or specifying the number 
#' of cells \code{x=3, y=4} as done by \code{grid}. Additionally, x- and y-locations of grid-lines can be specified,
#' e.g. \code{x=1:10, y=seq(0,10,2)}.
#' 
#' @param x         (integer, numeric) single integer specifies number of cells, numeric vector specifies vertical grid-lines
#' @param y         (integer, numeric) single integer specifies number of cells, numeric vector specifies horizontal grid-lines
#' @param col       (character) color of grid-lines
#' @param lwd       (integer) line width of grid-lines
#' @param lty       (integer) line type of grid-lines
#' 
#' @aliases addGrid
#' 
#' @author Andre Schuetzenmeister \email{andre.schuetzenmeister@@roche.com}

add_grid <- function(x=NULL, y=NULL, col="lightgray", lwd=1L, lty=3L)
{
	if(all(is.null(c(x,y))) || all(length(c(x,y))<2))               # call grid function
		grid(nx=x, ny=y, col=col, lwd=lwd, lty=lty)
	else
	{
		if(length(x) == 0)                                          # NULL
			xticks <- axTicks(side=1)
		else if(length(x) == 1)
		{
			U <- par("usr")
			xticks <- seq.int(U[1L], U[2L], length.out = x + 1)
		}
		else
			xticks <- x
		
		if(length(y) == 0)                                          # NULL
			yticks <- axTicks(side=2)
		else if(length(y) == 1)
		{
			U <- par("usr")
			yticks <- seq.int(U[3L], U[4L], length.out = y + 1)
		}
		else
			yticks <- y
		
		abline(v=xticks, col=col, lwd=lwd, lty=lty)
		abline(h=yticks, col=col, lwd=lwd, lty=lty)
	}
}





#' Convert Color-Name or RGB-Code to Possibly Semi-Transparent RGB-code.
#' 
#' Function takes the name of a color and converts it into the rgb space. Parameter "alpha" allows
#' to specify the transparency within [0,1], 0 meaning completey transparent and 1 meaning completey
#' opaque. If an RGB-code is provided and alpha != 1, the RGB-code of the transparency adapted color 
#' will be returned.
#' 
#' @param col (character) name of the color to be converted/transformed into RGB-space (code). Only
#'               those colors can be used which are part of the set returned by function colors(). Defaults
#'               to "black".
#' @param alpha (numeric) value specifying the transparency to be used, 0 = completely transparent, 
#'               1 = opaque.
#' 
#' @return RGB-code
#' 
#' @author Andre Schuetzenmeister \email{andre.schuetzenmeister@@roche.com}
#' 
#' @aliases as.rgb
#' 
#' @examples 
#' # convert character string representing a color to RGB-code
#' # using alpha-channel of .25 (75\% transparent)
#' as_rgb("red", alpha=.25)
#' 
#' # same thing now using the RGB-code of red (alpha=1, i.e. as_rgb("red"))
#' as_rgb("#FF0000FF", alpha=.25)

as_rgb <- function(col="black", alpha=1)
{
	if(length(col) > 1 && (length(alpha) == 1 || length(alpha) < length(col)))         # unclear which alpha to use or only one alpha specified
	{
		if(length(alpha) < length(col) && length(alpha) > 1)
			warning("Multiple (but too few) 'alpha' specified! Only use 'alpha[1]' for each color!")
		return(sapply(col, as_rgb, alpha=alpha[1]))
	}
	if(length(col) > 1 && length(col) <= length(alpha))                                 # process each color separately
	{
		res <- character()
		for(i in 1:length(col))
			res <- c(res, as_rgb(col[i], alpha[i]))
		return(res)
	}
	if( col %in% colors() )
		return( rgb(t(col2rgb(col))/255, alpha=alpha) )
	else
	{
		col <- sub("#", "", col)
		R <- as.numeric(paste("0x", substr(col, 1,2), sep=""))
		G <- as.numeric(paste("0x", substr(col, 3,4), sep=""))
		B <- as.numeric(paste("0x", substr(col, 5,6), sep=""))
		return( rgb(R/255, G/255, B/255, alpha=alpha, maxColorValue=1) )
	}        
}



#' Add Legend to Right Margin.
#' 
#' This function accepts all parameters applicable in and forwards them to function \code{\link{legend}}.
#' There will be only made some modifications to the X-coordinate ensuring that the legend is plotted in
#' the right margin of the graphic device. Make sure that you have reserved sufficient space in the right
#' margin, e.g. 'plot.VFP(....., mar=c(4,5,4,10))'.
#' 
#' @param x			(character, numeric) either one of the character strings "center","bottomright", "bottom", "bottomleft", 
#' 					"left", "topleft", "top", "topright", "right" or a numeric values specifying the X-coordinate in user
#' 					coordinates
#' @param y			(numeric) value specifying the Y-coordiante in user coordinates, only used in case 'x' is numeric
#' @param offset	(numeric) value in [0, 0.5] specifying the offset as fraction in regard to width of the right margin
#' @param ...		all parameters applicable in function \code{\link{legend}}
#' 
#' @author Andre Schuetzenmeister \email{andre.schuetzenmeister@@roche.com}
#' 
#' @aliases legend.rm
#' 
#' @examples 
#' \donttest{
#' library(VCA)
#' data(VCAdata1)
#' # perform VCA-anaylsis
#' lst <- anovaVCA(y~(device+lot)/day/run, VCAdata1, by="sample")
#' # transform list of VCA-objects into required matrix
#' mat <- get_mat(lst)		# automatically selects "total"
#' mat
#' 
#' # fit all 9 models batch-wise
#' res <- fit_vfp(model.no=1:9, Data=mat)
#' 
#' plot(res, mar=c(5.1, 4.1, 4.1,15), Crit=NULL)
#' 
#' legend_rm(cex=1.25, text.font=10,
#' 	 legend=c(
#'     paste0("AIC:    ", signif(as.numeric(res$AIC["Model_6"]), 4)),
#' 	   paste0("Dev:    ", signif(as.numeric(res$Deviance["Model_6"]), 4)),
#'     paste0("RSS:    ", signif(as.numeric(res$RSS["Model_6"]),4))))

#' }

legend_rm <- function(	x=c("center","bottomright", "bottom", "bottomleft", 
				"left", "topleft", "top", "topright", "right"), 
		y=NULL, offset=.05, ...)
{
	stopifnot(	is.numeric(x) || is.character(x) )
	if(is.character(x))
	{
		x <- match.arg(x[1], choices=c("center","bottomright", "bottom", "bottomleft", 
						"left", "topleft", "top", "topright", "right"))
	}else
	{
		stopifnot(is.numeric(y))
	}
	
	par(xpd=TRUE)
	args <- list(...)
	
	USR  <- par("usr")
	PLT  <- par("plt")
	FIG  <- par("fig")
	
	wrm <- FIG[2] - PLT[2]						# width right margin
	hrm	<- PLT[4] - PLT[3]
	
	if(is.character(x))
	{
		X.orig	<- x
		xjust 	<- 0.5							# defaults to center
		x 		<- PLT[2] + 0.5 * wrm
		yjust   <- 0.5
		y		<- PLT[3] + 0.5 * hrm
		
		if(grepl("left", X.orig))
		{
			xjust 	<- 0
			x 		<- PLT[2] + offset * wrm
		}		
		if(grepl("right", X.orig))
		{
			xjust 	<- 1
			x 		<- PLT[2] + (1-offset) * wrm
		}
		if(grepl("top", X.orig))
		{
			yjust 	<- 1
			y 		<- PLT[3] + (1-offset) * hrm
		}
		if(grepl("bottom", X.orig))
		{
			yjust 	<- 0
			y 		<- PLT[3] + offset * hrm
		}
	}
	
	x <- grconvertX(x, from="nic", to="user")
	y <- grconvertY(y, from="nic", to="user")
	
	args$x 		<- x
	args$y 		<- y
	args$xjust 	<- xjust
	args$yjust 	<- yjust
	
	do.call(legend, args)
	par(xpd=FALSE)
}