API tools faq
paste
Guest User

Untitled

a guest
Jun 3rd, 2015
233
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 29.32 KB | None | 0 0
raw download clone embed print report
  1. ### Title: Back to basics: High quality plots using base R graphics
  2. ### An interactive tutorial for the Davis R Users Group meeting on April 24, 2015
  3. ###
  4. ### Date created: 20150418
  5. ### Last updated: 20150423
  6. ###
  7. ### Author: Michael Koontz
  8. ### Email: mikoontz@gmail.com
  9. ### Twitter: @michaeljkoontz
  10. ###
  11. ### Purpose: Introduce basic to intermediate plotting capabilities of base R graphics
  12. ###
  13. ### Basic methods
  14. ### 1) Basic scatterplot and labeling a plot (Line 44)
  15. ### 2) Plotting groups in different ways on the same plot (Line 72)
  16. ### 3) Adding a legend (Line 120)
  17. ### 4) Adding a best fit line (Line 150)
  18. ### 5) Adding a 95% confidence interval (Line 150)
  19. ### 6) Shaded confidence intervals (Line 223)
  20. ### 7) Bar plots (Line 260)
  21. ### 8) Error bars (Line 274)
  22. ###
  23. ### Intermediate methods
  24. ### 1) Using other graphics devices like pdf() (Line 324)
  25. ### 2) Using par() for multipanel plots (Line 380)
  26. ### 3) Using par() for margin adjustments (Line 438)
  27. ### 4) Using axis() and mtext() (Line 484)
  28. ### 5) Pretty print from plotmath (Line 601)
  29.  
  30.  
  31. # We'll start with the very tractable 'trees' dataset, which is built into R. It describes the girth, height, and volume of 31 black cherry trees.
  32.  
  33. trees
  34. dim(trees)
  35. head(trees)
  36.  
  37. # Remember how we access the columns of a data.frame:
  38. trees$Girth
  39. trees$Volume
  40.  
  41. #Basic plot function takes x argument and a y argument
  42. #Default plot type is points, but you can change it to lines or both points and lines by adding the 'type' argument
  43.  
  44. plot(x=trees$Girth, y=trees$Volume)
  45. plot(x=trees$Girth, y=trees$Volume, type="l")
  46. plot(x=trees$Girth, y=trees$Volume, type="b")
  47.  
  48. # pch: 'plotting character' changes the type of point that is used (default is an open circle); remember pch=19!
  49. plot(x=trees$Girth, y=trees$Volume, pch=19)
  50.  
  51. # main: adds a title
  52. plot(x=trees$Girth, y=trees$Volume, pch=19, main="Girth vs. Volume for Black Cherry Trees")
  53.  
  54. # xlab: adds an x axis label
  55. plot(x=trees$Girth, y=trees$Volume, pch=19, main="Girth vs. Volume for Black Cherry Trees", xlab="Tree Girth (in)")
  56.  
  57. # ylab: adds a y axis label
  58. plot(x=trees$Girth, y=trees$Volume, pch=19, main="Girth vs. Volume for Black Cherry Trees", xlab="Tree Girth (in)", ylab="Tree Volume (cu ft)")
  59.  
  60. # las: rotates axis labels; las=1 makes them all parallel to reading direction
  61. plot(x=trees$Girth, y=trees$Volume, pch=19, main="Girth vs. Volume for Black Cherry Trees", xlab="Tree Girth (in)", ylab="Tree Volume (cu ft)", las=1)
  62.  
  63. # col: select a color for the plotting characters
  64. plot(x=trees$Girth, y=trees$Volume, pch=19, main="Girth vs. Volume for Black Cherry Trees", xlab="Tree Girth (in)", ylab="Tree Volume (cu ft)", las=1, col="blue")
  65.  
  66. # We can use the c() function to make a vector and have several colors, plotting characters, etc. per plot.
  67.  
  68. plot(x=trees$Girth, y=trees$Volume, pch=19, main="Girth vs. Volume for Black Cherry Trees", xlab="Tree Girth (in)", ylab="Tree Volume (cu ft)", las=1, col=c("black", "blue"))
  69.  
  70. plot(x=trees$Girth, y=trees$Volume, pch=c(1,19), main="Girth vs. Volume for Black Cherry Trees", xlab="Tree Girth (in)", ylab="Tree Volume (cu ft)", las=1, col="blue")
  71.  
  72. #------------
  73. # Plotting by group
  74. #------------
  75.  
  76. ### Those different colors and plotting characters that we just saw were arbitrary. The 2-element vector of colors or plotting characters just repeats for the whole data frame. What if we want to have more meaningful coloration, with a different color for each group?
  77.  
  78. ### We'll use the iris dataset to illustrate one way to do this. This dataframe describes the sepal length, sepal width, petal length, petal width, and species for 150 different irises.
  79.  
  80. # First look at the data:
  81. iris
  82. head(iris)
  83. dim(iris)
  84. str(iris)
  85.  
  86. # We can extend the idea of passing a vector of colors to the col= argument in the plot() function call.
  87.  
  88. # Let's cheat first, and see what the finished product will look like. First I define a new object with the three colors that I want to use.
  89. plot.colors <- c("violet", "purple", "blue")
  90. plot.colors
  91.  
  92. # Here's the cheating bit: I just looked at this dataframe and saw that there are exactly 50 observations for each species.
  93. iris
  94.  
  95. # I use the repeat function, rep() and the each= argument, to create a new vector with each element of plot.colors repeated 50 times in turn.
  96. color.vector <- rep(x=plot.colors, each=50)
  97. color.vector
  98.  
  99. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, col=color.vector)
  100.  
  101. # Notice the lengths of the x-vector, the y-vector, and the color vector are all the same.
  102. length(iris$Petal.Length)
  103. length(iris$Sepal.Length)
  104. length(color.vector)
  105.  
  106. # What if we want to automate the process? We can take advantage of the fact that the Species column is a factor.
  107. head(iris)
  108. iris$Species
  109. str(iris)
  110. as.numeric(iris$Species)
  111.  
  112. plot.colors <- c("violet", "purple", "blue")
  113.  
  114. color.vector <- plot.colors[iris$Species]
  115.  
  116. dev.off() # Just clearing the present plots
  117.  
  118. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, col=color.vector, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", las=1)
  119.  
  120. #-----------
  121. # Let's add a legend
  122. #-----------
  123.  
  124. # We use the legend() function to add a legend to an existing plot
  125.  
  126. legend("topleft", pch=19, col=plot.colors, legend=unique(iris$Species))
  127.  
  128. # You can customize the legend if you wish.
  129. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, col=color.vector, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", las=1)
  130.  
  131. # Here I pass a character vector to the legend= argument so that I can include the first letter of the species name
  132. # The bty="n" argument suppresses the border around the legend. (A personal preference)
  133. legend("topleft", pch=19, col=plot.colors, legend=c("I. setosa", "I. versicolor", "I. virginica"), bty="n")
  134.  
  135. # Italicize the labels in the legend using text.font=3
  136. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, col=color.vector, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", las=1)
  137.  
  138. legend("topleft", pch=19, col=plot.colors, legend=c("I. setosa", "I. versicolor", "I. virginica"), bty="n", text.font=3)
  139.  
  140.  
  141. #------------------
  142. # Test yourself.
  143. #------------------
  144.  
  145. #Using the ToothGrowth dataset built into R, plot the tooth length (the len column) as a function of the vitamin C dosage (the dose column). Use a different color for each method of administering the vitamin C (the supp column).
  146.  
  147. head(ToothGrowth)
  148.  
  149.  
  150. #------------------
  151. # Add a linear best fit line and confidence interval to a plot
  152. #------------------
  153.  
  154. # We'll use a simple linear regression for this, but the general recipe is the same every time.
  155. # The Recipe
  156. # 1) Estimate the parameters of the best fit line
  157. # 2) Make up your own x values that span the range of your data
  158. # 3) Get your y values by applying your mathematical model (e.g. a straight line) with the best fit parameters to your fabricated x values
  159. # 4) Plot these new y values against your fabricated x values.
  160.  
  161. # Save your model fit to an object. Here, we model Sepal.Length as a function of Petal.Length
  162. model1 <- lm(Sepal.Length ~ Petal.Length, data=iris)
  163.  
  164. # Now we have the parameter estimates for our y=ax+b line. The estimate for (Intercept) is b, and the estimate for Petal.Length is a.
  165. summary(model1)
  166.  
  167. # Make up our own x values; put them in a dataframe!
  168. range(iris$Petal.Length)
  169. xvals <- seq(from=1, to=7, by=0.1)
  170. xvals
  171. df <- data.frame(Petal.Length=xvals)
  172. df
  173.  
  174. # Take advantage of the predict() function, which returns a matrix with one row for each of your new x values, and 3 columns: 'fit' is the expected y value, 'lwr' is the lower 95% confidence interval, and 'upr' is the upper 95% confidence interval. Note this only works this seamlessly using the predict() function on an lm() model
  175.  
  176. CI <- predict(model1, newdata=df, interval="confidence")
  177. CI <- as.data.frame(CI) # Coerce the matrix to a dataframe, so we can access the column using the $ operator.
  178.  
  179. dim(df) # 61 x-values
  180. head(df) # Here are the first 6 of them
  181. dim(CI) # 61 records, 3 columns
  182. head(CI) # y-values, lower 95% CI bounds, and upper 95% CI bounds
  183.  
  184.  
  185. # Plot the actual data (in the iris dataframe)
  186. # This code copied from above
  187. plot.colors <- c("violet", "purple", "blue")
  188. color.vector <- plot.colors[iris$Species]
  189.  
  190. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", col=color.vector)
  191.  
  192. # Plot our best fit line. The x values are the Petal.Length column from the 'df' dataframe, and the y values are the 'fit' column from the CI dataframe.
  193. # Note that I use the lines() function, which just adds features to an existing plot.
  194. # The lwd= argument changes the line width
  195. lines(x=df$Petal.Length, y=CI$fit, lwd=2)
  196.  
  197. # Plot the confidence intervals
  198. # The lty= argument changes the line type. There are 6 different line types, and you can just put a number 1 through 6 if you'd like. Default is "solid" (aka 1)
  199. lines(x=df$Petal.Length, y=CI$lwr, lwd=2, lty="dashed", col="red")
  200. lines(x=df$Petal.Length, y=CI$upr, lwd=2, lty="dashed", col="red")
  201.  
  202. #----------------------
  203. # Going further... Prediction/forecast intervals
  204. #----------------------
  205.  
  206. forecast <- predict(model1, newdata=df, interval="prediction")
  207. forecast <- as.data.frame(forecast)
  208. head(forecast)
  209.  
  210. # Plot the actual data (in the iris dataframe)
  211. # This code copied from above
  212. plot.colors <- c("violet", "purple", "blue")
  213. color.vector <- plot.colors[iris$Species]
  214.  
  215. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", col=color.vector)
  216.  
  217. # New code using the forecast object
  218. lines(x=df$Petal.Length, y=forecast$fit, lwd=2)
  219. lines(x=df$Petal.Length, y=forecast$lwr, lwd=2, col="red", lty="dashed")
  220. lines(x=df$Petal.Length, y=forecast$upr, lwd=2, col="red", lty="dashed")
  221.  
  222.  
  223. #-------------
  224. # Shaded region of forecast interval
  225. #-------------
  226.  
  227. # This plot() call is copied from above
  228. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", col=color.vector)
  229.  
  230. lines(x=df$Petal.Length, y=forecast$fit, lwd=2)
  231.  
  232. # New code
  233. polygon.x <- c(df$Petal.Length, rev(df$Petal.Length))
  234. polygon.y <- c(forecast$lwr, rev(forecast$upr))
  235.  
  236. # By default, R won't fill in the polygon
  237. polygon(x=polygon.x, y=polygon.y)
  238.  
  239. # But we also may not want an opaque polygon
  240. polygon(x=polygon.x, y=polygon.y, col='darkgrey')
  241.  
  242. # The adjustcolor() function is nice for 'turning down' the opacity. It takes a color and the opacity level as arguments. I also use border=NA to suppress the border of the polygon
  243.  
  244. # First recreate the plot
  245. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", col=color.vector)
  246.  
  247. lines(x=df$Petal.Length, y=forecast$fit, lwd=2)
  248.  
  249. polygon(x=polygon.x, y=polygon.y, col=adjustcolor("black", alpha.f=0.4), border=NA)
  250.  
  251. # Add our legend back in
  252. legend("topleft", pch=19, col=plot.colors, legend=c("I. setosa", "I. versicolor", "I. virginica"), bty="n", text.font=3)
  253.  
  254. #--------------
  255. # Test yourself.
  256. #--------------
  257. # Layer the 95% confidence interval as a shaded region on top of the iris data, the best fit line for a Sepal.Length~Petal.Length model, and the forecast interval.
  258.  
  259.  
  260. #------------
  261. # Barplots
  262. #------------
  263.  
  264. model2 <- lm(Sepal.Length ~ Species, data=iris)
  265.  
  266. bar.heights <- predict(model2, newdata=data.frame(Species=c("setosa", "versicolor", "virginica")))
  267.  
  268. # The basic barplot function
  269. barplot(bar.heights)
  270.  
  271. # Let's add some flair
  272. barplot(bar.heights, names.arg=c("I. setosa", "I. versicolor", "I. virginica"), las=1, col=adjustcolor(plot.colors, alpha.f=0.5), main="Sepal length for 3 Irises", ylab="Sepal length (cm)")
  273.  
  274. #---------------
  275. # Error bars
  276. #---------------
  277. # Adding error bars to our barplot. These can be added to scatter plots in a similar way.
  278.  
  279. # We'll plot error bars representing 5 standard errors so you can see them more easily.
  280.  
  281. d <- summary(model2)
  282. CI <- 5 * coef(d)[ ,'Std. Error']
  283. lwr <- bar.heights - CI
  284. upr <- bar.heights + CI
  285.  
  286. # I used the ylim= argument to pass a 2-element numeric vector specifying the y extent of the barplot. I added some extra room on the top to account for error bars.
  287. # Importantly, assign the barplot to an object. I called it 'b' but you can call it whatever you like.
  288.  
  289. b <- barplot(bar.heights, names.arg=c("I. setosa", "I. versicolor", "I. virginica"), las=1, ylim=c(0,7.5), col=adjustcolor(plot.colors, alpha.f=0.5), main="Sepal length for 3 Irises", ylab="Sepal length (cm)")
  290.  
  291. # The object that you called your barplot is interpretted by R as the x values in the middle of each bar
  292. b
  293.  
  294. # We'll use the arrows() function to add arrows to an existing plot. With some modifications, our arrows will have an arrowhead at each end (code=3), and the 'arrowhead' will actually be perpendicular to the arrow shaft (angle=90)
  295. # Specify where each arrow starts (x0= and y=) and ends (x1= and y1=)
  296. arrows(x0=b, x1=b, y0=lwr, y1=upr, code=3, angle=90, length=0.1)
  297.  
  298.  
  299. #---------------
  300. # Test yourself
  301. #---------------
  302. # These data represent survivorship of plant seedlings in 4 different treatments: ambient, watered, heated + watered, and heated. Make a bar plot with their 95% confidence intervals. Note these are asymmetric (more uncertainty above the mean than below) like what might come from a logistic regression model.
  303.  
  304. prop <- c(0.18, 0.25, 0.13, 0.05)
  305. asympLCL <- c(0.14, 0.20, 0.11, 0.035)
  306. asympUCL <- c(0.24, 0.33, 0.18, 0.09)
  307.  
  308. #---------------
  309. # Test yourself. Error bars on scatter plots.
  310. #---------------
  311. # The randomly generated data below are measurements of the number of the number of angels who get their wings as a function of the number of bells that have been rung. There is some uncertainty in measuring wing acquisition (represented as the offset from the sampled mean). How would you add error bars to a scatter plot?
  312.  
  313. set.seed(13)
  314. n <- 20 # Number of experimental trials
  315. a <- 12
  316. b <- 1.5
  317. xvals <- round(runif(n)*50)
  318. yvals <- round(a + b*xvals + rnorm(n, sd=5))
  319. offset <- rpois(n, lambda=10)
  320. lwr <- yvals - offset
  321. upr <- yvals + offset
  322.  
  323.  
  324. #-----------------
  325. # Base R plotting skills: Other graphics devices
  326. #-----------------
  327.  
  328. # 1) Use the pdf() graphics device (or png() or postscript()) to get a permanent record of your plot in the exact final format.
  329.  
  330. # First set your working directory. A working directory is where R assumes all input and output files that it interfaces with should be. When you send your code to a collaborator, they can just change the working directory once instead of changing the filepath for every input (e.g. reading data into R) or output (e.g. making a plot)
  331.  
  332. # This is how I set mine:
  333. # 1) Make sure your script file is saved in a folder particular for the given project
  334. # 2) Run the file.choose() function
  335. # 3) Navigate to your script file and open it
  336. # 4) Copy the file path in the console up until the actual file name (i.e. just the folder path)
  337. # 5) Paste that folder into the setwd() function in quotes
  338. # 6) Delete the file.choose() command
  339.  
  340. file.choose()
  341.  
  342. pdf("iris plot.pdf")
  343. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, xlab="Petal length", ylab="Sepal length", col=plot.colors[iris$Species])
  344.  
  345. lines(x=df$Petal.Length, y=forecast$fit, lwd=2)
  346.  
  347. polygon.x <- c(df$Petal.Length, rev(df$Petal.Length))
  348. polygon.y <- c(forecast$lwr, rev(forecast$upr))
  349.  
  350. polygon(x=polygon.x, y=polygon.y, col=adjustcolor("black", alpha.f=0.4), border=NA)
  351.  
  352. # Add our legend back in
  353. legend("topleft", pch=19, col=plot.colors, legend=c("I. setosa", "I. versicolor", "I. virginica"), bty="n", text.font=3)
  354.  
  355. # Importantly, turn the device off at the end of your plotting block to complete the .pdf file.
  356. dev.off()
  357.  
  358.  
  359. #--------------
  360. # Figure for a paper
  361. #--------------
  362. # To make a figure for a paper that meets particular size and style guidelines, add some more arguments to the pdf() function call.
  363.  
  364. pdf("iris plot for paper.pdf", width=3.3, height=3.3, pointsize=9, family="Times")
  365.  
  366. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, xlab="Petal length", ylab="Sepal length", col=plot.colors[iris$Species])
  367. lines(x=df$Petal.Length, y=forecast$fit, lwd=2)
  368.  
  369. polygon.x <- c(df$Petal.Length, rev(df$Petal.Length))
  370. polygon.y <- c(forecast$lwr, rev(forecast$upr))
  371.  
  372. polygon(x=polygon.x, y=polygon.y, col=adjustcolor("black", alpha.f=0.4), border=NA)
  373.  
  374. # Add our legend back in
  375. legend("topleft", pch=19, col=plot.colors, legend=c("I. setosa", "I. versicolor", "I. virginica"), bty="n", text.font=3)
  376.  
  377. # Importantly, turn the device off at the end of your plotting block to complete the .pdf file.
  378. dev.off()
  379.  
  380. #-----------------
  381. # Base R plotting skills: Managing par()
  382. #-----------------
  383.  
  384. # 2) You can change some fundamental plotting parameters by using par() before a plot.
  385.  
  386. # This is how you can:
  387. # Make multipanel plots
  388. # Give your plots more room in the inner margins
  389. # Give your plots some room in the outer margins
  390.  
  391. #---------------
  392. # Multi-panel plots
  393. #---------------
  394.  
  395. setosa <- subset(iris, subset=Species=="setosa")
  396. versicolor <- subset(iris, subset=Species=="versicolor")
  397. virginica <- subset(iris, subset=Species=="virginica")
  398.  
  399. # Using par(mfrow=c(number of rows, number of columns))
  400. par(mfrow=c(2,2))
  401.  
  402. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab="Sepal Length", xlab="Petal Length")
  403.  
  404. plot(x=setosa$Petal.Length, y=setosa$Sepal.Length, pch=19, las=1, col=plot.colors[1], main="I. setosa", ylab="Sepal Length", xlab="Petal Length")
  405.  
  406. plot(x=versicolor$Petal.Length, y=versicolor$Sepal.Length, pch=19, las=1, col=plot.colors[2], main="I. versicolor", ylab="Sepal Length", xlab="Petal Length")
  407.  
  408. plot(x=virginica$Petal.Length, y=virginica$Sepal.Length, pch=19, las=1, col=plot.colors[3], main="I. virginica", ylab="Sepal Length", xlab="Petal Length")
  409.  
  410. dev.off()
  411.  
  412. # Using layout(matrix())
  413. plot.matrix <- matrix(c(1,1,2,3), nrow=2, ncol=2)
  414. plot.matrix
  415.  
  416. layout(mat=plot.matrix)
  417.  
  418. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab="Sepal Length", xlab="Petal Length")
  419.  
  420. plot(x=versicolor$Petal.Length, y=versicolor$Sepal.Length, pch=19, las=1, col=plot.colors[2], main="I. versicolor", ylab="Sepal Length", xlab="Petal Length")
  421.  
  422. plot(x=virginica$Petal.Length, y=virginica$Sepal.Length, pch=19, las=1, col=plot.colors[3], main="I. virginica", ylab="Sepal Length", xlab="Petal Length")
  423.  
  424. # Or we can adjust the widths of the plots (the heights, too). The widths= argument should be the same length as the number of columns while the heights= argument should be the same length as the number of rows. The values represent multipliers.
  425.  
  426. layout(mat=plot.matrix, widths=c(2,1))
  427. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab="Sepal Length", xlab="Petal Length")
  428.  
  429. plot(x=versicolor$Petal.Length, y=versicolor$Sepal.Length, pch=19, las=1, col=plot.colors[2], main="I. versicolor", ylab="Sepal Length", xlab="Petal Length")
  430.  
  431. plot(x=virginica$Petal.Length, y=virginica$Sepal.Length, pch=19, las=1, col=plot.colors[3], main="I. virginica", ylab="Sepal Length", xlab="Petal Length")
  432.  
  433. dev.off()
  434.  
  435. #----------------
  436. # Specifying margins
  437. #----------------
  438. # We want to make sure not to waste white space, so we can specify the plotting margins using par(mar=c(bottom, left, top, right)). This can be important if your labels aren't fitting on the plot.
  439. # The default is par(mar=c(5.1, 4.1, 4.1, 2.1)).
  440.  
  441. # Let's say I want to make the text bigger because the plot is for a presentation. We can use the cex (character expander) family of arguments in the plot() call.
  442. # cex= makes the plotted points larger or smaller
  443. # cex.main= makes the title larger
  444. # cex.lab= makes the axis labels larger
  445. # cex.axis= makes the numbers of the axis larger
  446.  
  447. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab="Sepal Length", xlab="Petal Length")
  448.  
  449. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab="Sepal Length", xlab="Petal Length", cex=2, cex.main=4, cex.lab=3, cex.axis=1.5)
  450.  
  451. # A solution
  452. par(mar=c(5,6,5,2))
  453. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab="Sepal Length", xlab="Petal Length", cex=2, cex.main=4, cex.lab=3, cex.axis=1.5)
  454.  
  455. # Note that the par(mgp=c(labels, tick marks, axis line)) method will adjust the location of those components away from the edge of the plot.
  456. # Default is par(mgp=c(3,1,0))
  457.  
  458. par(mar=c(6,7,5,2), mgp=c(4, 1.5, 0))
  459. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab="Sepal Length", xlab="Petal Length", cex=2, cex.main=4, cex.lab=3, cex.axis=1.5)
  460.  
  461. #--------------
  462. # Specifying outer margins
  463. #--------------
  464. # Outer margins are especially great for multi-panel plots. Let's recreate the one that we had earlier and add some outer margins with par(oma=c(bottom, left, top, right))
  465. # The default is par(oma=c(0,0,0,0))
  466.  
  467. dev.off()
  468. par(mfrow=c(2,2), oma=c(5,5,5,0), mar=c(2,2,4,2))
  469.  
  470. # Note here that I use the font.main= argument in each of these plot calls to make the title in normal font for the first panel and then italicized for the other panels (it is bold by default; equivalent to a font.main=2)
  471.  
  472. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab=NA, xlab=NA, cex.main=2, font.main=1)
  473.  
  474. plot(x=setosa$Petal.Length, y=setosa$Sepal.Length, pch=19, las=1, col=plot.colors[1], main="I. setosa", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  475.  
  476. plot(x=versicolor$Petal.Length, y=versicolor$Sepal.Length, pch=19, las=1, col=plot.colors[2], main="I. versicolor", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  477.  
  478. plot(x=virginica$Petal.Length, y=virginica$Sepal.Length, pch=19, las=1, col=plot.colors[3], main="I. virginica", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  479.  
  480.  
  481. #-----------------
  482. # Base R plotting skills: Using mtext() and axis()
  483. #-----------------
  484. # The two add on functions that I find myself using most often are mtext(), which adds text to a margin, and axis(), which draws in a user-specified axis.
  485.  
  486. # mtext() lets me:
  487. # 1) Make axis labels that are shared between more than one panel
  488. # 2) rotate the y-axis label so it is easier to read for presentations
  489.  
  490. # axis() lets me:
  491. # 1) Specify exactly how many tick marks and where they'll be
  492. # 2) Specify exactly how those tick marks will be labelled
  493.  
  494. # We actually set up the previous 4-panel plot nicely to see the functionality of mtext() and axis(). Let's recreate that plot.
  495.  
  496. par(mfrow=c(2,2), oma=c(5,6,5,0), mar=c(2,2,3,2))
  497.  
  498. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab=NA, xlab=NA, cex.main=2, font.main=1)
  499.  
  500. plot(x=setosa$Petal.Length, y=setosa$Sepal.Length, pch=19, las=1, col=plot.colors[1], main="I. setosa", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  501.  
  502. plot(x=versicolor$Petal.Length, y=versicolor$Sepal.Length, pch=19, las=1, col=plot.colors[2], main="I. versicolor", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  503.  
  504. plot(x=virginica$Petal.Length, y=virginica$Sepal.Length, pch=19, las=1, col=plot.colors[3], main="I. virginica", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  505.  
  506. # By default, an mtext() call will put an axis label on the specified side (1=bottom, 2=left, 3=top, 4=right) of the most recently created panel. It also puts it right at the margin. Use the line= argument to move it further from the margin however far you choose. These line units are the same units as those used by par(mar=c()).
  507.  
  508. mtext(side=1, text="Petal Length", line=3)
  509.  
  510. # We can specify outer=TRUE to put this text on the outer margins. This works only because we used par(oma=c()) to make the outer margins greater than 0 on some sides.
  511.  
  512. # Recreate the plot
  513. par(mfrow=c(2,2), oma=c(5,6,5,0), mar=c(2,2,3,2))
  514.  
  515. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab=NA, xlab=NA, cex.main=2, font.main=1)
  516.  
  517. plot(x=setosa$Petal.Length, y=setosa$Sepal.Length, pch=19, las=1, col=plot.colors[1], main="I. setosa", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  518.  
  519. plot(x=versicolor$Petal.Length, y=versicolor$Sepal.Length, pch=19, las=1, col=plot.colors[2], main="I. versicolor", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  520.  
  521. plot(x=virginica$Petal.Length, y=virginica$Sepal.Length, pch=19, las=1, col=plot.colors[3], main="I. virginica", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  522.  
  523. # We can just use cex= to make the text bigger here (rather than cex.lab=) because we're already in the mtext() function
  524. mtext(side=1, text="Petal Length", outer=TRUE, line=3, cex=2)
  525.  
  526. # Do the same for the y axis
  527. mtext(side=2, text="Sepal Length", outer=TRUE, line=3, cex=2)
  528.  
  529. # And the same for an overall title
  530. mtext(side=3, text="Sepal length vs. petal length for Irises", outer=TRUE, line=2, cex=2.1)
  531.  
  532. # We can also use the las= argument to rotate axis text using mtext(). We can make the y-axis label in the reading direction this way.
  533.  
  534. # Recreate the plot. Add some extra outer margin space to the left.
  535.  
  536. par(mfrow=c(2,2), oma=c(5,9,5,0), mar=c(2,2,3,2))
  537.  
  538. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="All Irises", ylab=NA, xlab=NA, cex.main=2, font.main=1)
  539.  
  540. plot(x=setosa$Petal.Length, y=setosa$Sepal.Length, pch=19, las=1, col=plot.colors[1], main="I. setosa", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  541.  
  542. plot(x=versicolor$Petal.Length, y=versicolor$Sepal.Length, pch=19, las=1, col=plot.colors[2], main="I. versicolor", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  543.  
  544. plot(x=virginica$Petal.Length, y=virginica$Sepal.Length, pch=19, las=1, col=plot.colors[3], main="I. virginica", ylab=NA, xlab=NA, cex.main=2, font.main=3)
  545.  
  546. mtext(side=1, text="Petal Length", outer=TRUE, line=3, cex=2)
  547.  
  548. # Use las=1 here to rotate the text. Use \n to insert a carriage return.
  549. mtext(side=2, text="Sepal\nLength", outer=TRUE, line=1.5, cex=2, las=1)
  550.  
  551. # I added an \n in the middle of the character string to insert a carriage return. I also changed the line= argument from 2 to 0.
  552. mtext(side=3, text="Sepal length vs. petal length\nfor Irises", outer=TRUE, line=0, cex=2.1)
  553.  
  554. #-----------------
  555. # Test yourself
  556. #-----------------
  557. # Add a best fit line and shaded 95% confidence interval to each panel of a 4-panel Iris plot similar to the one we just created.
  558.  
  559.  
  560. #------------------
  561. # Using axis()
  562. #------------------
  563. # Put the axis tick marks where you want them, and label them how you like
  564.  
  565. # First reset the plotting device.
  566. dev.off()
  567.  
  568. # Here's the original plot
  569. plot.colors <- c("violet", "purple", "blue")
  570. color.vector <- plot.colors[iris$Species]
  571.  
  572. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", col=color.vector)
  573.  
  574. # In the plot() call, suppress the x-axis with xaxt="n" and the y-axis with yaxt="n". Then we can draw them in how we'd like.
  575.  
  576. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", col=color.vector, yaxt="n")
  577.  
  578. # Use axis() in a similar way to using mtext(). The side= argument specifies where you want the axis (1=bottom, 2=left, 3=top, 4=right)
  579.  
  580. axis(side=2, at=c(4.5, 6, 7.5), las=1)
  581.  
  582. # Or change how the tick marks are labeled
  583. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", col=color.vector, yaxt="n")
  584.  
  585. axis(side=2, at=c(4.5, 5, 6, 7.5), labels=c("4.5", "Critical\nPoint", "6", "7.5"), las=1)
  586.  
  587. # Or meaningful points
  588. s.tick <- mean(setosa$Sepal.Length)
  589. ve.tick <- mean(versicolor$Sepal.Length)
  590. vi.tick <- mean(virginica$Sepal.Length)
  591.  
  592. ticks <- c(s.tick, ve.tick, vi.tick)
  593. ticks
  594.  
  595. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", col=color.vector, yaxt="n")
  596.  
  597. axis(side=2, at=ticks, labels=round(ticks,digits=2), las=1)
  598.  
  599. #-----------------
  600. # Base R plotting skills: Using Pretty Print
  601. #-----------------
  602. # R can render some nice typesetting for axis labels, titles, etc.
  603.  
  604. # Here's the breakdown (I'll use the main= argument for illustration):
  605. # 1) Want only text? Just use a text string.
  606. # ex: main="text here"
  607. # 2) Want evaluated statements and text? Use paste().
  608. # ex: main=paste("Mean sepal length =", mean(iris$Sepal.Length))
  609. # 3) Want special symbols? Use expression().
  610. # ex: main=expression(alpha)
  611. # 4) Want special symbols and text? Use expression(paste())
  612. # ex: main=expression(paste(mu, "Sepal length="))
  613. # 5) Want text, symbols, and evaluated statements? Use bquote(). A tilda (~) goes between each component. Wrap statements you want evaluated in .()
  614. # ex: main=bquote(mu[sepal ~ length] ~ "=" .(mean(iris$Sepal.Length)))
  615.  
  616. # See ?plotmath for all of the fancy typesetting you can use.
  617.  
  618. # Original plot with just text as title
  619. plot(x=iris$Petal.Length, y=iris$Sepal.Length, pch=19, las=1, main="Iris sepal length vs. petal length", xlab="Petal length", ylab="Sepal length", col=color.vector)
  620.  
  621. # Evaluated statement AND text using paste()
  622. plot( x=iris$Petal.Length,
  623. y=iris$Sepal.Length,
  624. pch=19,
  625. las=1,
  626. main=paste("Mean sepal length =", mean(iris$Sepal.Length)),
  627. xlab="Petal length",
  628. ylab="Sepal length",
  629. col=color.vector)
  630.  
  631. # Just a special symbol? Use expression()
  632. plot( x=iris$Petal.Length,
  633. y=iris$Sepal.Length,
  634. pch=19,
  635. las=1,
  636. main=expression(mu),
  637. xlab="Petal length",
  638. ylab="Sepal length",
  639. col=color.vector)
  640.  
  641. # Special symbol plus text? Use expression(paste())
  642. plot( x=iris$Petal.Length,
  643. y=iris$Sepal.Length,
  644. pch=19,
  645. las=1,
  646. main=expression(paste(mu["sepal length"])),
  647. xlab="Petal length",
  648. ylab="Sepal length",
  649. col=color.vector)
  650.  
  651. # Special symbol plus text plus evaluated statements? Use bquote(). Put a ~ in between each component. Equal signs go in quotes. Wrap statements you want evaluated in .()
  652. plot( x=iris$Petal.Length,
  653. y=iris$Sepal.Length,
  654. pch=19,
  655. las=1,
  656. main=bquote(mu[sepal ~ length] ~ "=" ~ .(mean(iris$Sepal.Length))),
  657. xlab="Petal length",
  658. ylab="Sepal length",
  659. col=color.vector)
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!