library(tidyverse)── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.5.0 ✔ purrr 0.3.4
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.0 ✔ stringr 1.4.1
✔ readr 2.1.2 ✔ forcats 0.5.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()library(ggplot2)
library(broom)Today we will use data from a study I published in 2021. The data are on github. Read this in as a csv and use it for the lab!
data1 <- read.csv("https://raw.githubusercontent.com/jbaumann3/Belize-RT-Baumann-et-al-2021/main/rt_master_data.csv")
head(data1) timepoint bw_date ID species colony coralID colony_ID home_reef
1 T0 12/11/2017 PSTR 1 1 PSTR P1 1 NSP1 NS
2 T0 12/11/2017 PSTR 1 2 PSTR P1 2 NSP1 NS
3 T0 12/11/2017 PSTR 1 3 PSTR P1 3 NSP1 NS
4 T0 12/11/2017 PSTR 1 4 PSTR P1 4 NSP1 NS
5 T0 12/11/2017 PSTR 1 5 PSTR P1 5 NSP1 NS
6 T0 12/11/2017 PSTR 1 6 PSTR P1 6 NSP1 NS
transplant_reef collection_site transplant_site transplant pale.or.bleached
1 NS False Caye False Caye NS Native N
2 NS False Caye False Caye NS Native N
3 NS False Caye False Caye NS Native N
4 NS False Caye False Caye NS Native N
5 NS False Caye False Caye NS Native N
6 NS False Caye False Caye NS Native N
partial.mortality mortality parent.colony.length parent.colony.width
1 N alive 19 14
2 N alive 19 14
3 N alive 19 14
4 N alive 19 14
5 N alive 19 14
6 N alive 19 14
parent.colony.height parent.colony.area parent.colony.volume bw
1 8 266 2128 16.99472
2 8 266 2128 28.50501
3 8 266 2128 16.50844
4 8 266 2128 13.15963
5 8 266 2128 14.14723
6 8 266 2128 22.20343
bw_sd bw_se surface.area sa_sd sa_se chla symb_den symb_sd
1 0.014142136 0.010 21.85694 0.028311047 0.016345390 NA NA NA
2 0.028284271 0.020 20.92004 0.266249295 0.153719102 NA NA NA
3 0.007071068 0.005 17.12328 0.192485165 0.111131362 NA NA NA
4 0.021213203 0.015 12.47212 0.115079269 0.066441047 NA NA NA
5 0.000000000 0.000 11.49451 0.029319816 0.016927804 NA NA NA
6 0.007071068 0.005 13.84078 0.004490647 0.002592676 NA NA NA
symb_se carbs_sa carbs_dw lipids_sa lipids_dw protein_dw protein_sa
1 NA NA NA NA NA NA NA
2 NA NA NA NA NA NA NA
3 NA NA NA NA NA NA NA
4 NA NA NA NA NA NA NA
5 NA NA NA NA NA NA NA
6 NA NA NA NA NA NA NAWe are interested in coral surface area and we have a column called surface area (with the space). R kind of hates spaces. We CAN just call the variable surface area all the time and it will work. OR we could rename the column or make a new column. Pick your favorite solution–> Do something to the dataframe that eliminates the need for us to use backticks (`) aroundsurface area`. In other words– remove that space somehow.
#used read.csvThe data above are relatively complex. There is A LOT going on, as is typical of large datasets. For now, what you need to know is that these are coral data from Belize. We did an experiment in which corals were collected, cut into smaller fragments, and placed into one of 4 treatments (in the ocean). Corals were collected from a Nearshore Reef (False Caye) or Offshore Reef (Silk Caye), they were then either left in their home environment (nearshore native / offshore native) OR transplanted to the nearshore (nearshore transplant) or offshore (offshore transplant). MANY columns are not relevant for us, but BW, which stands for buoyant weight, represents the mass of each coral and the column ‘surface area’ represents the surface area of each coral.
H0: There is no relationship between coral surface area and coral buoyant weight in SSID corals (you will need to filter so only have 1 species of coral)
data2 <- data1 %>% select("species", "bw", "surface.area")
data2 <- data2 %>% filter(species == "PSTR")
head(data2) species bw surface.area
1 PSTR 16.99472 21.85694
2 PSTR 28.50501 20.92004
3 PSTR 16.50844 17.12328
4 PSTR 13.15963 12.47212
5 PSTR 14.14723 11.49451
6 PSTR 22.20343 13.84078stats <- lm(data = data2, bw ~ surface.area)
summary(stats)
Call:
lm(formula = bw ~ surface.area, data = data2)
Residuals:
Min 1Q Median 3Q Max
-20.5937 -4.1412 -0.5371 3.0975 29.7465
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.21829 0.62256 6.776 3.4e-11 ***
surface.area 0.84963 0.02665 31.877 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.349 on 514 degrees of freedom
(15 observations deleted due to missingness)
Multiple R-squared: 0.6641, Adjusted R-squared: 0.6634
F-statistic: 1016 on 1 and 514 DF, p-value: < 2.2e-16plot <- ggplot(data2, aes(x = surface.area, y = bw)) +
geom_point() +
geom_smooth(method='lm', se=FALSE) +
labs(title = "Coral Surface Area vs Coral Buoyant Weight", x = "surface area", y = "buoyant weight")
plot`geom_smooth()` using formula = 'y ~ x'Warning: Removed 15 rows containing non-finite outside the scale range
(`stat_smooth()`).Warning: Removed 15 rows containing missing values or values outside the scale range
(`geom_point()`).
#there is a relationshipDo we reject H0? Why or why not- cite your graph and table
We reject the null hypothesis. First based on the graph we see a clear linear relationship between surface area and buoyant weight and then looking at the p-value on the table, it’s below 0.05, so we can reject the null and state that there is a relationship surface area and buoyant weight.
H0: There is no effect of surface area or species on buoyant weight
additive <- lm(bw ~ surface.area + species, data = data1)
tidy(additive)# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 4.67 0.606 7.72 2.96e- 14
2 surface.area 0.828 0.0248 33.3 1.46e-163
3 speciesSSID 2.89 0.464 6.23 6.92e- 10add_graph<-additive %>%
augment() %>%
ggplot(aes(x=surface.area, y=bw, color=species))+
geom_point()+
geom_line(aes(y=.fitted))+
theme_classic()
add_graph
Do we reject or fail to reject H0? What is significant and what is not? How do we know if there is a species effect? Use both the graph and the table to explain.
my response:
Based on the very small p-values (far under 0.05) shown in the table, we see that it is likely that the surface area has a result on buoyancy weight and also that species likely has an effect on buoyancy weight. This is displayed on the graph because it’s clear that a greater surface area results more often in a greater buoyancy weight. Similarly, looking at species there is an observable difference in their weights based off of the two species’ lines of best fit. Therefore, we reject the null hypothesis.
my response
stats <- lm(data = data1, bw ~ surface.area * species)
summary(stats)
Call:
lm(formula = bw ~ surface.area * species, data = data1)
Residuals:
Min 1Q Median 3Q Max
-24.150 -4.871 -0.570 3.781 32.856
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.21829 0.69708 6.051 2.04e-09 ***
surface.area 0.84963 0.02984 28.470 < 2e-16 ***
speciesSSID 4.17376 1.07855 3.870 0.000116 ***
surface.area:speciesSSID -0.07092 0.05383 -1.318 0.187925
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.109 on 983 degrees of freedom
(77 observations deleted due to missingness)
Multiple R-squared: 0.5311, Adjusted R-squared: 0.5297
F-statistic: 371.2 on 3 and 983 DF, p-value: < 2.2e-16plot <- ggplot(data1, aes(x = surface.area, y = bw, color = species)) +
geom_point() +
geom_smooth(method='lm', se=FALSE) +
labs(title = "Coral Surface Area and Coral Species vs Coral Buoyant Weight", x = "surface area", y = "buoyant weight", color = "Coral Species:")
plot`geom_smooth()` using formula = 'y ~ x'Warning: Removed 77 rows containing non-finite outside the scale range
(`stat_smooth()`).Warning: Removed 77 rows containing missing values or values outside the scale range
(`geom_point()`).
Interactive mode: We know based on the table’s p-values that surface area has an effect on buoyancy and likewise species also has an effect on buoyancy. However, the interaction between surface area and species doesn’t have a small enough p-value, so there’s no evidence that this interaction has an effect on buoyancy. Thus we fail to reject the null hypothesis. This is reflected in the graph. We see that the different species as surface area increases has a similar affect on buoyant weight.
The difference between the two is that an interactive multiple regression model accounts for the interaction between variables x and z itself having an effect on y while the additive model tests only how x is affecting y and how z is affecting why. There is not a significant interaction between the explanatory variables in part(b) based on two high a p-value for the interaction. As stated above, the graph displays how the species and surface area together has a similar affect on buoyant weight.
my response
anova <- aov(bw ~ surface.area * species, data = data1)
summary(anova) Df Sum Sq Mean Sq F value Pr(>F)
surface.area 1 54224 54224 1073.041 < 2e-16 ***
species 1 1962 1962 38.836 6.83e-10 ***
surface.area:species 1 88 88 1.736 0.188
Residuals 983 49674 51
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
77 observations deleted due to missingnessI noticed here that while the p-values of from the lm() and the aov() are relatively similar, they are slightly different, though in this case it doesn’t change the effect on the null hypothesis. We also have the F-value to work with on aov() which we didn’t have on lm(). The TukeyHSD test didn’t work here because one of the independent variables, Species, is not numeric.
Your goal is to recreate this graph:
Below I lay out a step by step process that can help you get to this point. You are welcome to recreate this graph however you’d like though. Make sure you answer the questions in part C after you make the graph!
#code for graph abovegraph <- ggplot(data = data1, aes(x = timepoint, y = bw)) +
geom_point() +
facet_wrap(~species) +
labs(x = "Time Point", y = "Buoyant Weight (g)") +
theme_classic() +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
strip.text.x = element_text(size = 14))
graphWarning: Removed 35 rows containing missing values or values outside the scale range
(`geom_point()`).
Do you notice any trends? Write a response in your .qmd
my response:
graph <- ggplot(data = data1, aes(x = timepoint, y = bw, color = transplant)) +
geom_point() +
geom_jitter(width = 0.15)+
scale_color_manual(values=c("#003366", "#CC0000", "#006633", "#330066", "#006666")) +
facet_wrap(~species) +
labs(x = "Time Point", y = "Buoyant Weight (g)") +
theme_classic() +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
strip.text.x = element_text(size = 14))
graphWarning: Removed 35 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 35 rows containing missing values or values outside the scale range
(`geom_point()`).
Now what do you see? my response:
df2 <- data1 %>%
group_by(timepoint, transplant) %>%
summarise_at(.vars = vars(bw, bw_se, bw_sd), .funs = mean, na.rm = TRUE)
head(df2)# A tibble: 6 × 5
# Groups: timepoint [2]
timepoint transplant bw bw_se bw_sd
<chr> <chr> <dbl> <dbl> <dbl>
1 T0 control NaN NaN NaN
2 T0 NS Native 20.6 0.0164 0.0233
3 T0 NS Transplant 14.2 0.00958 0.0136
4 T0 OS Native 14.2 0.00978 0.0138
5 T0 OS Transplant 21.4 0.00906 0.0130
6 T1 NS Native 22.9 0.00389 0.00552graph <- ggplot(data = data1, aes(x = timepoint, y = bw, color = transplant)) +
geom_point(alpha = 0.15) +
geom_jitter(alpha = 0.15, width = 0.25) +
geom_point(data = df2, aes(x = timepoint, y = bw), size = 3) +
geom_jitter(data = df2, alpha = 0.15, width = 0.25) +
geom_line(data = df2) +
geom_errorbar(data = df2, aes(x = timepoint, ymin = bw-(bw_sd), ymax = bw+(bw_sd)), width=.2, position=position_dodge(0.05)) +
scale_color_manual(values=c("#003366", "red", "chartreuse4", "#330066", "cyan4")) +
facet_wrap(~species) +
labs(x = "Time Point", y = "Buoyant Weight (g)") +
theme_classic() +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
strip.text.x = element_text(size = 14))
graphWarning: Removed 35 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 35 rows containing missing values or values outside the scale range
(`geom_point()`).Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_line()`).`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
I added error bars and the means for each transplant. My issue now is that I can’t figure out how to connect them. Whenever, I try to use geom_line(), the different line doesn’t connect across different time points.
What do you see now? Are there any noticeable differences between transplant treatments? Write your final thoughts here:
my response:
Yes there are differences between transplant treatments. NS Native and OS Native grow at the same rate my with NS Native has a higher buoyant weight the whole time. NS Transplant and OS Transplant switch buyant weight as time point increases.