## Making a graphic *function* using R's grid graphics, ## advanced example with gridlines replacing y axis ## (Lineplot of a regression with CI) ## ## HELPER FUNCTIONS: ## mmest.fit() ## lighten() ## myGridPlot() ## ## Data and example from Iversen & Soskice 2003 ## ## Chris Adolph ## 1 February 2018 # Helper functions: MM-estimator fitting mmest.fit <- function(y, x, ci, quiet=TRUE) { require(MASS) y <- y[order(x)] x <- x[order(x)] result <- rlm(y~x, method="MM") if (!quiet) print(result) fit <- list(x=x) fit$y <- result$fitted.values fit$lower <- fit$upper <- NULL if (length(na.omit(ci))>0) for (i in 1:length(ci)) { pred <- predict(result, interval="confidence", level=ci[i]) fit$lower <- cbind(fit$lower, pred[,2]) fit$upper <- cbind(fit$upper, pred[,3]) } fit } # Here's an effort at a color lightener that could use work lighten <- function(col, pct=0.75, alpha=1){ if (abs(pct)>1) { warning("invalid pct in lighten(); must be between 0 and 1.") pcol <- col2rgb(col)/255 } else { col <- col2rgb(col)/255 if (pct>0) { pcol <- col + pct*(1-col) } else { pcol <- col*pct } } rgb(pcol[1], pcol[2], pcol[3], alpha) } myGridPlot <- function(x, y, limits, xat, yat, logX=FALSE, logY=FALSE, col="black", mainTitle="Plot title", xaxisTitle="x", yaxisTitle="y", file="gridplot") { ## Check limits and logging if (logX) { if (any(limits[1:2]<=0)) { stop("Can't log x-axis with non-positive limits") } else { limits[1:2] <- log(limits[1:2]) } x <- log(x) xat <- log(xat) } if (logY) { if (any(limits[3:4]<=0)) { stop("Can't log y-axis with non-positive limits") } else { limits[3:4] <- log(limits[3:4]) } y <- log(y) yat <- log(yat) } ## Open a new pdf device pdf(paste0(file, ".pdf"), width=5, height=5) ## Set up layout widths wd <- unit.c(unit(6,"strheight", yaxisTitle), unit(5,"null")) ## Set up layout heights ht <- unit.c(unit(4,"strheight", mainTitle), unit(5,"null"), unit(6,"strheight", xaxisTitle)) ## Create layout overlay <- grid.layout(nrow=3, ncol=2, widths=wd, heights=ht, respect=TRUE) ## Display layout ##grid.show.layout(overlay) ## "Push" layout to create initial viewport pushViewport(viewport(layout=overlay)) ## Push the main title viewport pushViewport(viewport(layout.pos.col=2, layout.pos.row=1, xscale=c(0,1), yscale=c(0,1), gp=gpar(fontsize=12), name="maintitle", clip="on" ) ) ## Note the use of a grid primitive grid.text(mainTitle, x=unit(0.5,"npc"), # Why NPC? y=unit(0.5,"npc"), gp=gpar(fontface="bold") ) ## Go back up to the top level Viewport upViewport(1) ## Push the Y-axis title viewport pushViewport(viewport(layout.pos.col=1, layout.pos.row=2, xscale=c(0, 1), yscale=c(0, 1), gp=gpar(fontsize=12), name="ytitle", clip="on" )) ## Plot the actual y-axis title grid.text(yaxisTitle, x=unit(0.15, "npc"), y=unit(0.5, "npc"), rot=90 ) ## Go back up to the top level Viewport upViewport(1) ## Push the X-axis title viewport pushViewport(viewport(layout.pos.col=2, layout.pos.row=3, xscale=c(0, 1), yscale=c(0, 1), gp=gpar(fontsize=12), name="xtitle", clip="on" )) ## Plot the actual x axis title grid.text(xaxisTitle, x=unit(0.5, "npc"), y=unit(0.25, "npc") ) ## Go back up to the top level Viewport upViewport(1) ## Push the main plot Viewport. Note the scales pushViewport(viewport(layout.pos.col=2, layout.pos.row=2, xscale=c(limits[1], limits[2]), yscale=c(limits[3], limits[4]), gp=gpar(fontsize=12), name="mainplot", clip="on" ) ) ## get the fit from the data fit <- mmest.fit(y=y, x=x, ci=0.95) ## Make the x-coord of a CI polygon xpoly <- c(fit$x, rev(fit$x), fit$x[1]) ## Make the y-coord of a CI polygon ypoly <- c(fit$lower, rev(fit$upper), fit$lower[1]) ## Choose a pastel color for the polygon lightercol <- lighten(col) ## Plot the polygon first grid.polygon(x=xpoly, y=ypoly, gp=gpar(col=lightercol, border=FALSE, fill=lightercol), default.units="native") ## Go back up to the top level Viewport upViewport(1) ## Push the main plot Viewport. Note the clip setting pushViewport(viewport(layout.pos.col=2, layout.pos.row=2, xscale=c(limits[1], limits[2]), yscale=c(limits[3], limits[4]), gp=gpar(fontsize=12), name="mainplot", clip="off" ) ) ## Add light horizontal lines at ticks ## (Best to do after polygons & before other features) for (i in 1:length(yat)) grid.lines(x=unit( c(0,1), "npc"), y=unit( rep(yat[i], 2), "native"), gp=gpar(col="gray80", lwd=0.6)) ## Go back up to the top level Viewport upViewport(1) ## Push the main plot Viewport. Note the scales pushViewport(viewport(layout.pos.col=2, layout.pos.row=2, xscale=c(limits[1], limits[2]), yscale=c(limits[3], limits[4]), gp=gpar(fontsize=12), name="mainplot", clip="on" ) ) ## Then plot the regression line grid.lines(x=fit$x, y=fit$y, gp=gpar(lty="solid", col=col), default.units="native") ## If you want a box around the graph, ## now's the time to plot it ##grid.rect(gp=gpar(linejoin="round")) ## Return to layout level upViewport(1) ## Wait! We haven't made axes! ## Axes plot facing out of the Viewport pushViewport(viewport(layout.pos.col=2, layout.pos.row=2, xscale=c(limits[1], limits[2]), yscale=c(limits[3], limits[4]), gp=gpar(fontsize=12), name="mainplot", clip="off" )) ## Make an axis as a "grob" (graphical object), but don't plot yet yaxis <- yaxisGrob(at = yat, # Where to put ticks label = TRUE, # Only takes logical values! main = TRUE # Top (TRUE) or bottom (FALSE) ) ## Edit the y-axis to show unlogged units if (logY) { yaxis <- editGrob(yaxis, gPath("labels"), label=exp(yat) ) } ## Edit the y-axis to remove the major line and tick lines yaxis <- removeGrob(yaxis, gPath("major")) yaxis <- removeGrob(yaxis, gPath("ticks")) ## Edit the y-axis to make the labels snug yaxis <- editGrob(yaxis, gPath("labels"), x=unit(0, "npc") - unit(0.5, "char") ## the left edge of the plot - a small spacer ) ## Draw the y-axis grob grid.draw(yaxis) ## Make an x-axis as a "grob" xaxis <- xaxisGrob(at = xat, label = TRUE, main = TRUE ) ## Edit the x-axis to show unlogged units if (logX) { xaxis <- editGrob(xaxis, gPath("labels"), label=exp(xat) ) } ## Edit the x-axis to fit the plot xaxis <- editGrob(xaxis, gPath("major"), x=unit(c(0, 1), "npc") ) ## Draw the x-axis grob grid.draw(xaxis) ## Return to layout level upViewport(1) ## Save the graph! dev.off() }