/
supp_figure.Rmd
128 lines (104 loc) · 4.68 KB
/
supp_figure.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
```{r }
# Creating a figure describing the relationship between bay laurel density and symptomatic leaf count from 2004-2012
# Load libraries
library(plyr)
library(dplyr)
library(ggplot2)
library(grid)
library(tidyr)
library(knitr)
setwd("~/GitHub/ggplot_figures")
# Read in data ####
bay_dens <- readRDS("bay_dens.rds")
str(bay_dens)
summary(bay_dens)
bay_slc <- readRDS("bay_slc.rds")
str(bay_slc)
summary(bay_slc)
# Tidy data ####
bay_dens <- bay_dens %>% gather(year, density, umca_ct_04:umca_ct_11) %>%
mutate(year = extract_numeric(year), year = as.factor(year+2000))
bay_slc <- bay_slc %>% gather(year, slc, cum_slc_04:cum_slc_11) %>%
mutate(year = extract_numeric(year), year = as.factor(year+2000))
bay_vars1 <- inner_join(bay_dens,bay_slc, by = c("plotid","year"))
str(bay_vars1)
head(bay_vars1)
#saveRDS(bay_vars1, file = "bay_variables.rds")
#bay_vars1 <- readRDS("bay_variables.rds")
# Check variable distributions ####
par(mfrow=c(2,1))
hist(bay_vars1$density, main = "Bay laurel density")
hist(bay_vars1$slc, main = "Bay laurel symptomatic leaf count")
```
These show some pretty severe right skew, so I tried a log transformation:
```{r Log-transformed histograms, fig.width=6, fig.height=9, echo=TRUE}
par(mfrow=c(2,1))
hist(log(bay_vars1$density), main = "Bay laurel density")
hist(log(bay_vars1$slc), main = "Bay laurel symptomatic leaf count")
```
```{r Version 1: Plot leaf count vs density faceted by year}
# Plot data ####
df1 <- bay_vars1 %>%
group_by(plotid, year)
p1 <- qplot(log(density+1), log(slc+1), data = df1, color = year, facets = year~., xlab = "Bay Laurel Density", ylab = "Symptomatic Leaf Count", main = "ln(Symptomatic Leaf Count) vs ln(Bay Laurel Density)")
# Add smooth line to get a better idea of the magnitude of change ####
p1 + geom_smooth(aes(group=year), method="lm", size = 0.5)
# Calculate the regression equation and r-square for each linear model ####
reg_eqn <- function(df) {
y = log(df$slc+1)
x = log(df$density+1)
m <- lm(y ~ x, data = df)
a <- coef(m)[1]
b <- coef(m)[2]
r2 <- summary(m)$r.squared
expr <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(a, digits = 2),
b = format(b, digits = 2),
r2 = format(r2, digits = 3)))
c(lab = as.character(as.expression(expr)))
}
eqns <- ddply(bay_vars1, "year", reg_eqn)
# Add the equations to the plot ####
p1 + geom_smooth(aes(group=year), method="lm", size = 0.5) +
geom_text(aes(x = 0, y = 3, label = lab), data = eqns, parse = T, hjust = 0, vjust = -1.5) +
theme_bw() +
theme(legend.position="none") +
theme(panel.margin = unit(0.7,"lines"))
```
```{r Final version ggplot code}
# Finalized plotting code ####
# Hat-tip to Noam Ross' blog post here <http://www.noamross.net/blog/2013/11/20/formatting-plots-for-pubs.html>
p1a <- ggplot(data = bay_vars1 %>% group_by(year),
aes(x = log(density+1), y = log(slc+1))) +
geom_point(size=3, alpha = 0.3) +
geom_smooth(method="lm", size = 0.5, color = "black") +
geom_text(aes(x = 0, y = 3, label = lab), data = eqns, parse = T,
hjust = 0, vjust = -1.2) +
facet_grid(year~.)
p1a + theme_bw() +
theme(legend.position = "none") +
theme(panel.margin = unit(0.7, "lines")) +
theme(panel.grid=element_blank()) +
labs(x = "Bay Laurel Density", y = "Bay Laurel Inoculum Load") +
theme(text = element_text(size = 16)) +
theme(axis.title.x = element_text(size=16)) +
theme(axis.title.y = element_text(size=16, angle=90)) +
theme(axis.ticks = element_line(size = 1.1)) +
theme(axis.ticks.length = unit(.2, "cm"))
ggsave(filename = "haas_supp_update1.pdf")
ggsave(filename = "haas_supp_update1.tiff", width = 6, height = 8, units = "in", dpi = 300)
ggsave(filename = "haas_supp_update1.png")
```
A couple of alternative ways to export the plot
```{r exporting plot}
# You can also do this from the plotting device
#`ggsave` defaults to saving the last plot that you displayed, and for a default size uses the size of the current graphics device. It also guesses the type of graphics device from the extension. This means the only argument you need to supply is the filename
ggsave("supp_figure/slc_vs_dens.svg", width = 6, height = 9)
#This is a `base` approach for exporting the final figure to PDF
pdf(file = "supp_figure/slc_vs_dens.pdf", width = 6, height = 9, useDingbats=F)
p1a + theme_bw() +
theme(legend.position = "none") +
theme(panel.margin = unit(0.7, "lines")) +
labs(x = "Log bay laurel density", y = "Log bay laurel symptomatic leaf count")
dev.off()
```