--- /dev/null
+#' @title Add a curve from a fitted linear model and a label to a plot.\r
+#'\r
+#' @description \code{stat_poly_eq} fits a polynomial and generates a label with an equation\r
+#' and/or coefficient of determination (R^2).\r
+#'\r
+#' @param mapping The aesthetic mapping, usually constructed with\r
+#' \code{\link[ggplot2]{aes}} or \code{\link[ggplot2]{aes_string}}. Only needs\r
+#' to be set at the layer level if you are overriding the plot defaults.\r
+#' @param data A layer specific dataset - only needed if you want to override\r
+#' the plot defaults.\r
+#' @param geom The geometric object to use display the data\r
+#' @param position The position adjustment to use for overlapping points on this\r
+#' layer\r
+#' @param show.legend logical. Should this layer be included in the legends?\r
+#' \code{NA}, the default, includes if any aesthetics are mapped. \code{FALSE}\r
+#' never includes, and \code{TRUE} always includes.\r
+#' @param inherit.aes If \code{FALSE}, overrides the default aesthetics, rather\r
+#' than combining with them. This is most useful for helper functions that\r
+#' define both data and aesthetics and shouldn't inherit behaviour from the\r
+#' default plot specification, e.g. \code{\link[ggplot2]{borders}}.\r
+#' @param ... other arguments passed on to \code{\link[ggplot2]{layer}}. This\r
+#' can include aesthetics whose values you want to set, not map. See\r
+#' \code{\link[ggplot2]{layer}} for more details.\r
+#' @param na.rm a logical indicating whether NA values should be stripped\r
+#' before the computation proceeds.\r
+#' @param formula a formula object\r
+#' @param eq.with.lhs logical indicating whether lhs of equation is to be\r
+#' included in label.\r
+#'\r
+#' @details This stat can be used to automatically annotate a plot with R^2,\r
+#' adjusted R^2 or the fitted model equation. It supports only linear models\r
+#' fitted with function \code{lm()}. The R^2 and adjusted R^2 annotations can be\r
+#' used with any linear model formula. The fitted equation label is correclty\r
+#' generated for polynomials or quasi-polynomials through the origin. Model\r
+#' formulas can use \code{poly()} or be defined algebraically with terms of\r
+#' powers of increasing magnitude with no missing intermediate terms, except\r
+#' possibly for the intercept indicated by "- 1" or "-1" in the formula. The\r
+#' validity of the \code{formula} is not checked in the current implementation,\r
+#' and for this reason the default aesthetics sets R^2 as label for the\r
+#' annotation. This stat only generates the label, the predicted values need\r
+#' to be sepearately added to the plot, so to make sure that the same model\r
+#' formula is used in all steps it is best to save the formula as an object\r
+#' and supply this object as argument to the different statistics.\r
+#'\r
+#' @section Computed variables:\r
+#' \describe{ \item{x}{x position for left edge}\r
+#' \item{y}{y position near upper edge}\r
+#' \item{eq.label}{equation for the\r
+#' fitted polynomial as a character string to be parsed}\r
+#' \item{rr.label}{\eqn{R^2} of the fitted model as a character string to be parsed}\r
+#' \item{adj.rr.label}{Adjusted \eqn{R^2} of the fitted model as a character string\r
+#' to be parsed}\r
+#' \item{hjust}{Set to zero to override the default of the "text" geom.}}\r
+#'\r
+#' @examples\r
+#' library(ggplot2)\r
+#' # generate artificial data\r
+#' set.seed(4321)\r
+#' x <- 1:100\r
+#' y <- (x + x^2 + x^3) + rnorm(length(x), mean = 0, sd = mean(x^3) / 4)\r
+#' my.data <- data.frame(x, y, group = c("A", "B"), y2 = y * c(0.5,2))\r
+#' # give a name to a formula\r
+#' formula <- y ~ poly(x, 3, raw = TRUE)\r
+#' # plot\r
+#' ggplot(my.data, aes(x, y)) +\r
+#' geom_point() +\r
+#' geom_smooth(method = "lm", formula = formula) +\r
+#' stat_poly_eq(formula = formula, parse = TRUE)\r
+#'\r
+#' @export\r
+#'\r
+stat_poly_eq <- function(mapping = NULL, data = NULL, geom = "text",\r
+ formula = NULL, eq.with.lhs = TRUE,\r
+ position = "identity",\r
+ na.rm = FALSE, show.legend = FALSE,\r
+ inherit.aes = TRUE, ...) {\r
+ ggplot2::layer(\r
+ stat = StatPolyEq, data = data, mapping = mapping, geom = geom,\r
+ position = position, show.legend = show.legend, inherit.aes = inherit.aes,\r
+ params = list(formula = formula, eq.with.lhs = eq.with.lhs,\r
+ na.rm = na.rm,\r
+ ...)\r
+ )\r
+}\r
+\r
+#' @rdname ggpmisc-ggproto\r
+#' @format NULL\r
+#' @usage NULL\r
+#' @export\r
+StatPolyEq <-\r
+ ggplot2::ggproto("StatPolyEq", ggplot2::Stat,\r
+ compute_group = function(data,\r
+ scales,\r
+ formula,\r
+ eq.with.lhs) {\r
+ mf <- lm(formula, data)\r
+ coefs <- coef(mf)\r
+ formula.rhs.chr <- as.character(formula)[3]\r
+ if (grepl("-1", formula.rhs.chr) || grepl("- 1", formula.rhs.chr)) {\r
+ coefs <- c(0, coefs)\r
+ }\r
+ rr <- summary(mf)$r.squared\r
+ adj.rr <- summary(mf)$adj.r.squared\r
+ eq.char <- as.character(signif(polynom::as.polynomial(coefs), 3))\r
+ if (eq.with.lhs) {\r
+ eq.char <- paste("italic(y)", eq.char, sep = "~`=`~")\r
+ }\r
+ rr.char <- format(rr, digits = 2)\r
+ adj.rr.char <- format(adj.rr, digits = 2)\r
+ data.frame(x = min(data$x),\r
+ y = max(data$y) - 0.1 * diff(range(data$y)),\r
+ eq.label = gsub("x", "~italic(x)", eq.char, fixed = TRUE),\r
+ rr.label = paste("italic(R)^2", rr.char, sep = "~`=`~"),\r
+ adj.rr.label = paste("italic(R)[adj]^2",\r
+ adj.rr.char, sep = "~`=`~"),\r
+ hjust = 0)\r
+ },\r
+ default_aes =\r
+ ggplot2::aes(label = ..rr.label.., hjust = ..hjust..),\r
+ required_aes = c("x", "y")\r
+ )\r
+\r
--- /dev/null
+to.org <- function(x){
+ sub("([0-9.-]+)([eE])\\+?(-?)0*([0-9]+)","$\\1\\\\times 10^{\\3\\4}$",
+ x)
+}
+to.latex <- function(x){
+ sub("([0-9]+)([eE])\\+?(-?)0*([0-9]+)","\\1\\\\ensuremath{\\\\times 10^{\\3\\4}}",
+ x)
+}
+cleanup.tolatex <- function(output) {
+ gsub("\\\\textrm\\{e\\}(-?)0*(\\d+)","\\\\ensuremath{\\\\times 10^{\\1\\2}}",output);
+}
+
+trimws <- function(x,left=TRUE,right=TRUE){
+ result <- x
+ if(left)
+ result <- gsub('^\\s+','',x,perl=TRUE)
+ if(right)
+ result <- gsub('\\s+$','',x,perl=TRUE)
+
+ return(result)
+}
+
+
+print.summary.glm.xtable <- function (x, digits = max(3, getOption("digits") - 3), symbolic.cor = x$symbolic.cor, signif.stars = getOption("show.signif.stars"), ...) {
+
+ cat("\\begin{verbatim}\n")
+ cat("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"),
+ "\n\n", sep = "")
+ cat("Deviance Residuals: \n")
+ if (x$df.residual > 5) {
+ x$deviance.resid <- quantile(x$deviance.resid, na.rm = TRUE)
+ names(x$deviance.resid) <- c("Min", "1Q", "Median", "3Q",
+ "Max")
+ }
+ xx <- zapsmall(x$deviance.resid, digits + 1)
+ cat("\\end{verbatim}\n")
+ print(table.to.latex(xx,digits=digits))
+ cat("\\begin{verbatim}\n")
+ if (length(x$aliased) == 0L) {
+ cat("\nNo Coefficients\n")
+ }
+ else {
+ df <- if ("df" %in% names(x))
+ x[["df"]]
+ else NULL
+ if (!is.null(df) && (nsingular <- df[3L] - df[1L]))
+ cat("\nCoefficients: (", nsingular, " not defined because of singularities)\n",
+ sep = "")
+ else cat("\nCoefficients:\n")
+ coefs <- x$coefficients
+ if (!is.null(aliased <- x$aliased) && any(aliased)) {
+ cn <- names(aliased)
+ coefs <- matrix(NA, length(aliased), 4L, dimnames = list(cn,
+ colnames(coefs)))
+ coefs[!aliased, ] <- x$coefficients
+ }
+ cat("\\end{verbatim}\n")
+ colnames(coefs) <- gsub("(Pr\\(>\\|z\\|\\))","$\\1$",colnames(coefs),perl=TRUE)
+ if (nrow(coefs) > 15) {
+ print(table.to.latex(coefs,longtable=TRUE))
+ } else {
+ print(table.to.latex(coefs))
+ }
+ cat("\\begin{verbatim}\n")
+ }
+ cat("\n(Dispersion parameter for ", x$family$family, " family taken to be ",
+ format(x$dispersion), ")\n\n", apply(cbind(paste(format(c("Null",
+ "Residual"), justify = "right"), "deviance:"), format(unlist(x[c("null.deviance",
+ "deviance")]), digits = max(5, digits + 1)), " on",
+ format(unlist(x[c("df.null", "df.residual")])), " degrees of freedom\n"),
+ 1L, paste, collapse = " "), sep = "")
+ if (nzchar(mess <- naprint(x$na.action)))
+ cat(" (", mess, ")\n", sep = "")
+ cat("AIC: ", format(x$aic, digits = max(4, digits + 1)),
+ "\n\n", "Number of Fisher Scoring iterations: ", x$iter,
+ "\n", sep = "")
+ correl <- x$correlation
+ if (!is.null(correl)) {
+ p <- NCOL(correl)
+ if (p > 1) {
+ cat("\nCorrelation of Coefficients:\n")
+ if (is.logical(symbolic.cor) && symbolic.cor) {
+ print(symnum(correl, abbr.colnames = NULL))
+ }
+ else {
+ correl <- format(round(correl, 2), nsmall = 2,
+ digits = digits)
+ correl[!lower.tri(correl)] <- ""
+ print(correl[-1, -p, drop = FALSE], quote = FALSE)
+ }
+ }
+ }
+ cat("\\end{verbatim}\n")
+ cat("\n")
+ invisible(x)
+}
+
+
+
+table.to.latex <- function(table,
+ digits=NULL,
+ format=NULL,
+ scientific=NA,
+ nsmall=0,
+ colspec="r",
+ useBooktabs=TRUE,
+ longtable=FALSE,
+ caption=NULL,
+ label=NULL,
+ rownames=TRUE,
+ centering=FALSE,
+ latex.pos="",
+ toprule=if(useBooktabs) "\\toprule" else "\\hline\\hline",
+ midrule=if(useBooktabs) "\\midrule" else "\\hline",
+ bottomrule=if(useBooktabs) "\\bottomrule" else "\\hline\\hline",
+ cols.as.is=FALSE,
+ ...) {
+ ans <- ""
+ header <- c(if(rownames){""}else{NULL},colnames(table));
+ .format.function <- function(x){
+ format(if(is.null(format)||is.null(scientific)) {x} else {
+ if(!is.na(suppressWarnings(as.numeric(x)))) as.numeric(x) else x},
+ digits=digits,nsmall=nsmall,scientific=scientific,...)
+ }
+ if(is.null(dim(table))){
+ header <- names(table)
+ rownames <- FALSE
+ }
+ res <- rbind(header,
+ cbind(if(rownames){rownames(table)}else{NULL},
+ if(is.null(dim(table))){t(sapply(table,.format.function))} else
+ {apply(table,1:2,.format.function)}
+ ))
+ res[,] <- gsub("_","\\\\_",res[,])
+ if (!missing(cols.as.is)) {
+ res.no.change <-
+ rbind(header,
+ cbind(if(rownames){rownames(table)}else{NULL},
+ if(is.null(dim(table))){
+ t(sapply(table,function(x){x}))
+ } else {
+ apply(table,1:2,function(x){x})
+ }
+ ))
+ res[,cols.as.is] <- res.no.change[,cols.as.is]
+ }
+ ignore.rows <- -1
+ coefrows <- 2:NROW(res)
+ if(rownames) {
+ coefrows <- 1:NROW(res)
+ ignore.rows <- TRUE
+ }
+## res[coefrows,ignore.rows] <- sub("(\\*+)","$^{\\1}$",
+## res[coefrows,ignore.rows])
+ if (NROW(res) > 1) {
+ res[coefrows,ignore.rows] <- sub("([0-9]+)([eE])\\+?(-?)0*([0-9]+)","\\1\\\\ensuremath{\\\\times 10^{\\3\\4}}",
+ res[coefrows,ignore.rows])
+ }
+ tabspec <- rep("l",ncol(res))
+ if (length(colspec)==length(tabspec) &&
+ length(colspec) > 1
+ ) {
+ tabspec <- colspec
+ }
+
+ tabbegin <- paste("\\begin{tabular}{",paste(tabspec,collapse=""),"}",sep="")
+ tabend <- "\\end{tabular}\n"
+ if (longtable) {
+ tabbegin <- paste("\\begin{longtable}{",paste(tabspec,collapse=""),"}",sep="")
+ tabend <- "\\end{longtable}\n"
+ if (!is.null(label)) {
+ tabend <- paste(sep="","\\label{",label,"}\\\\","\n",
+ tabend
+ )
+ }
+ if(!is.null(caption))
+ tabend <- paste(sep="","\\caption{",caption,"}\\\\","\n",
+ tabend
+ )
+ } else if (!is.null(label)) {
+ tabbegin <- paste0("\\begin{table}",latex.pos,"\n",
+ if (centering) {"\\centering\n"} else {""},
+ tabbegin)
+ if (!is.null(caption)) {
+ tabend <- paste0(tabend,"\\caption{",caption,"}\n")
+ }
+ tabend <- paste0(tabend,"\\label{",label,"}\n",
+ "\\end{table}","\n")
+ }
+ ans <- c(ans,tabbegin)
+ header <- NULL
+ if(length(toprule))
+ header <- c(header,toprule)
+ if (!rownames)
+ header <- c(header,paste("\\multicolumn{1}{c}{",gsub("\\_"," ",trimws(res[1,1])),"}",sep=""))
+ for(j in 2:ncol(res))
+ header <- c(header,paste(" & \\multicolumn{1}{c}{",gsub("\\_"," ",trimws(res[1,j])),"}",sep=""))
+ header <- c(header,"\\\\")
+ if(length(midrule))
+ header <- c(header,midrule)
+ if (longtable)
+ header <- c(header,"\\endhead")
+ ans <- c(ans,header)
+ if (NROW(res) > 1) {
+ for(i in 2:NROW(res)) {
+ ans <- c(ans,
+ paste(paste(res[i,],collapse=" & "),"\\\\"))
+ }
+ }
+ if(length(bottomrule))
+ ans <- c(ans,bottomrule)
+ ans <- c(ans,tabend)
+ structure(ans,class="Latex")
+
+}