SDS300CR_Lab_2_template

Author

Julia Lewis-Borgia

Part 0: Read in packages we need

Make a code chunk below and read in any packages you may need.

library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.4
✔ ggplot2   3.5.0     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(palmerpenguins)
library(plotrix)
data(package = 'palmerpenguins')

Part 1: Reshape data and test the Intermediate Disturbance Hypothesis

  1. Read in this data from github
ocean_point <- read.csv("https://raw.githubusercontent.com/jbaumann3/BIOL234_Biostats_MHC/main/Ocean_Point_Data%202018.csv")
head(ocean_point)
  Number  trans_description year tidal_ht_m tidal_ht_cat_num bare_rock calothr
1      0 a_cobble_protected 2017      -0.13                1         0       0
2      0 a_cobble_protected 2017      -0.13                1         5       0
3      0 a_cobble_protected 2017      -0.13                1         5       0
4      0 a_cobble_protected 2018       0.38                2        10      12
5      0 a_cobble_protected 2018       0.38                2         0      93
6      0 a_cobble_protected 2018       0.78                2         0      65
  semibal mytilus Halichondria Spongomorpha chaetomo ulva_lac ascophyl lamin_di
1       0       0            0            0        0        0        0        0
2       0       0            0            0        0        5        0        0
3       0       0            0            0        0        0        0        0
4       8       0            0            0        0        0       10        0
5       7       0            0            0        0        0        0        0
6      30       0            0            0        0        0        0        0
  lamin_sa fucus_spir fucus_vesic fucus_dist fucus_sp chondrus mastocar
1        0          0           0          0        0        0       40
2        0          0           5          0        0        0       20
3        0          0           5          0        0        0       20
4        0          0          60          0        0        0        0
5        0          0           0          0        0        0        0
6        0          0           5          0        0        0        0
  petrocel corralina ceramium lithothamn Enteromorpha tectura litt_sax litt_lit
1        0         2        0         58            0       0        0       64
2        0         5        0         60            0       0        0        0
3        0        15        0         55            0       0        0        0
4        0         0        0          0            0       0        0      112
5        0         0        0          0            0       0        0      112
6        0         0        0          0            0       0        0        0
  litt_obt nucella carcinus Cancer.crabs
1        0       0        0            0
2        0       0        0            0
3       16       0        0            0
4       80       0        0            0
5       16       0        0            0
6        0       0        0            0
  1. Today we are interested in how much algal diversity is in a place given the number of grazers present. In short- is there a relationship between grazers (snail) and algae? Bare rock, carcinus, and cancer crabs columns can be dropped as we are not using them. Do this below
#drop columns we don't need!
part_ocean_point <- ocean_point %>% select(-one_of('bare_rock', 'carcinus', 'Cancer.crabs'))
head(part_ocean_point)
  Number  trans_description year tidal_ht_m tidal_ht_cat_num calothr semibal
1      0 a_cobble_protected 2017      -0.13                1       0       0
2      0 a_cobble_protected 2017      -0.13                1       0       0
3      0 a_cobble_protected 2017      -0.13                1       0       0
4      0 a_cobble_protected 2018       0.38                2      12       8
5      0 a_cobble_protected 2018       0.38                2      93       7
6      0 a_cobble_protected 2018       0.78                2      65      30
  mytilus Halichondria Spongomorpha chaetomo ulva_lac ascophyl lamin_di
1       0            0            0        0        0        0        0
2       0            0            0        0        5        0        0
3       0            0            0        0        0        0        0
4       0            0            0        0        0       10        0
5       0            0            0        0        0        0        0
6       0            0            0        0        0        0        0
  lamin_sa fucus_spir fucus_vesic fucus_dist fucus_sp chondrus mastocar
1        0          0           0          0        0        0       40
2        0          0           5          0        0        0       20
3        0          0           5          0        0        0       20
4        0          0          60          0        0        0        0
5        0          0           0          0        0        0        0
6        0          0           5          0        0        0        0
  petrocel corralina ceramium lithothamn Enteromorpha tectura litt_sax litt_lit
1        0         2        0         58            0       0        0       64
2        0         5        0         60            0       0        0        0
3        0        15        0         55            0       0        0        0
4        0         0        0          0            0       0        0      112
5        0         0        0          0            0       0        0      112
6        0         0        0          0            0       0        0        0
  litt_obt nucella
1        0       0
2        0       0
3       16       0
4       80       0
5       16       0
6        0       0
  1. The data are survey data from the rocky intertidal in coastal Maine. Students worked at a transect (trans_description) in different years (year), they measured how high they were above (or below) the low tide line (tidal_ht), those numerical tide heights were converted into a ranked list to make bins (tidal_height_cat_number, 1-8 or so), and then students recorded either the percent of 2D area that was covered by a thing (if it was something that didn’t move) OR how many individuals of a thing they saw (if a crab or snail or something). This abundance or cover data is in columnar form by species (or category) in the remaining columns. This works nicely when we are recording data on a clipboard in the field. BUT, R wants all of those names in 1 column and all of those numbers in a second column. We can use pivot_longer to convert this dataframe to long format.

Pivot the data to long format - one column for names and one for values

part_ocean_point_long <-part_ocean_point %>% pivot_longer(cols=calothr:nucella,
                    names_to='species',
                    values_to='value')
show(part_ocean_point_long) 
# A tibble: 4,030 × 7
   Number trans_description   year tidal_ht_m tidal_ht_cat_num species     value
    <int> <chr>              <int>      <dbl>            <int> <chr>       <dbl>
 1      0 a_cobble_protected  2017      -0.13                1 calothr         0
 2      0 a_cobble_protected  2017      -0.13                1 semibal         0
 3      0 a_cobble_protected  2017      -0.13                1 mytilus         0
 4      0 a_cobble_protected  2017      -0.13                1 Halichondr…     0
 5      0 a_cobble_protected  2017      -0.13                1 Spongomorp…     0
 6      0 a_cobble_protected  2017      -0.13                1 chaetomo        0
 7      0 a_cobble_protected  2017      -0.13                1 ulva_lac        0
 8      0 a_cobble_protected  2017      -0.13                1 ascophyl        0
 9      0 a_cobble_protected  2017      -0.13                1 lamin_di        0
10      0 a_cobble_protected  2017      -0.13                1 lamin_sa        0
# ℹ 4,020 more rows
  1. Make a new column to tell us if value = %cover or abundance. The grazes (animals) are counted in abundance and the algaes (everything else) are % cover. You may also want to make an algae vs grazer column. Do this below
algaes_list <- list("calothr", "Halichondria", "Spongomorpha", "chaetomo", "ulva_lac", "ascophyl", "lamin_di", "lamin_sa", "fucus_spir", "fucus_vesic", "fucus_dist", "fucus_sp", "chondrus", "mastocar", "petrocel", "corralina", "ceramium", "lithothamn","Enteromorpha")

part_ocean_point_long$alg_or_grz <- with(part_ocean_point_long, ifelse(species %in% algaes_list, "algae", "graze"))

part_ocean_point_long$pct_cov_or_abun <- with(part_ocean_point_long, ifelse(alg_or_grz == "algae", "percent cover", "abundance"))

col_order <- c('Number', 'trans_description', 'year',
               'species', 'value', 'pct_cov_or_abun', 'alg_or_grz', 'tidal_ht_m', 'tidal_ht_cat_num')

show(part_ocean_point_long)
# A tibble: 4,030 × 9
   Number trans_description   year tidal_ht_m tidal_ht_cat_num species     value
    <int> <chr>              <int>      <dbl>            <int> <chr>       <dbl>
 1      0 a_cobble_protected  2017      -0.13                1 calothr         0
 2      0 a_cobble_protected  2017      -0.13                1 semibal         0
 3      0 a_cobble_protected  2017      -0.13                1 mytilus         0
 4      0 a_cobble_protected  2017      -0.13                1 Halichondr…     0
 5      0 a_cobble_protected  2017      -0.13                1 Spongomorp…     0
 6      0 a_cobble_protected  2017      -0.13                1 chaetomo        0
 7      0 a_cobble_protected  2017      -0.13                1 ulva_lac        0
 8      0 a_cobble_protected  2017      -0.13                1 ascophyl        0
 9      0 a_cobble_protected  2017      -0.13                1 lamin_di        0
10      0 a_cobble_protected  2017      -0.13                1 lamin_sa        0
# ℹ 4,020 more rows
# ℹ 2 more variables: alg_or_grz <chr>, pct_cov_or_abun <chr>
  1. We are interested in a visual that assess the intermediate disturbance hypothesis, which states that there is some middle range of disturbance (grazing in this case) that leads to the most DIVERSE algal ecosystem. Calculate total grazer count and total algal diversity for each tidal height cat num.
algae_div <- part_ocean_point_long %>% 
  filter(alg_or_grz == "algae", value > 0)  %>%
  select(tidal_ht_cat_num, value, species) %>%
  group_by(tidal_ht_cat_num) %>% 
  summarize(count = n_distinct(species))
names(algae_div)[2] <- "Algae_Diversity"

graze_part <- part_ocean_point_long %>% 
  group_by(tidal_ht_cat_num) %>%
  summarize(across(alg_or_grz, list(Graze_Count = ~ sum(. == "graze"))))
names(graze_part)[2] <- "Graze_Count"

final_count = merge(graze_part, algae_div, on='tidal_ht_cat_num', all.y = TRUE)
final_count[is.na(final_count)] <- 0
final_count
   tidal_ht_cat_num Graze_Count Algae_Diversity
1                 1         147              10
2                 2          91               8
3                 3          77               7
4                 4         245               8
5                 5          84               4
6                 6          84               1
7                 7          84               1
8                 8          84               1
9                 9         105               1
10               10          84               1
  1. Graph to visually assess the IDH– grazer abundance on the X and algal diversity on the y! What do trends do you see? Do you see any evidence for the IDH? Why or why not?
idh_plot <- ggplot(final_count, aes(x = Graze_Count, y = Algae_Diversity)) + 
  geom_point(size = 3) + 
  labs(title = "Grazer Abundace vs Algal Diversity", x = "Grazer Abundance", y = "Algal Diversity")
  
idh_plot

There seems to be some support for the IDH as the highest amount of algal diversity occurs when grazer abundance seems to be in the middle of the range (at about 150) while the data points with a smaller or larger grazer count have a lower amount of algal diversity.

Part 2: dplyr/tidyverse pipelines and ANOVA

  1. Read in this coral cover data from github.

In the following chunk, starting with the raw data from github, build a dplyr/tidyverse pipeline that calculates mean coral cover by reef type. Don’t forget to calculate a relevant and useful error term for plotting

coral_cover <- read.csv("https://raw.githubusercontent.com/jbaumann3/BIOL234_Biostats_MHC/main/coralcover.csv")
head(coral_cover)
  site       type lat transect diver cc_percent
1    1  Back Reef   3        1     4  5.8435758
2    1  Back Reef   3        2     4  0.9505263
3    1  Back Reef   3        3     4  5.2423389
4    1  Back Reef   3        4     5  5.0040475
5    1  Back Reef   3        5     5  5.8954916
6    2 Patch Reef   3        1     4  5.2826190
means <- coral_cover %>%
  group_by(type) %>% 
summarize(mean = mean(cc_percent, na.rm=TRUE), st_error = std.error(cc_percent, na.rm=TRUE))
means
# A tibble: 3 × 3
  type        mean st_error
  <chr>      <dbl>    <dbl>
1 Back Reef  11.9     1.77 
2 Nearshore   3.16    0.531
3 Patch Reef 12.7     2.11 

2.Make a visual that allows us to assess the following Null Hypothesis: There is no difference in coral cover by reef type. You should make a graph that shows relevant means, error bars, and raw data

cc_by_reef_type2 <- ggplot(means, aes(x = type, y = mean, color = type)) +
  geom_point(color = "red") + 
    geom_errorbar(aes(x=type, ymin=mean-st_error, ymax=mean+st_error), width=0.4, alpha=0.9) + 
  geom_jitter(data = coral_cover, aes(x = type, y = cc_percent), alpha = 0.4) + 
  labs(title = "Coral Cover Percentage by Reef Type", x = "Reef Type", y = "Coral Cover Percentage") 
cc_by_reef_type2

  1. Run the correct statistical test to accompany the plot above. This should be a One-Way ANOVA + TukeyHSD test. Run these below AND interpret them. What do you see evidence for? Do we reject or fail to reject the null?
anova_test <- aov(cc_percent ~ type, data = coral_cover)
summary(anova_test)
            Df Sum Sq Mean Sq F value  Pr(>F)   
type         2   1172   585.9   6.683 0.00215 **
Residuals   74   6487    87.7                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tukey <- TukeyHSD(anova_test)
tukey
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = cc_percent ~ type, data = coral_cover)

$type
                           diff        lwr       upr     p adj
Nearshore-Back Reef  -8.7600964 -15.479776 -2.040416 0.0072335
Patch Reef-Back Reef  0.8217506  -5.010031  6.653533 0.9393628
Patch Reef-Nearshore  9.5818470   2.905196 16.258498 0.0027941

We first see from the anova test that we have a high F-value of 6.683 and a low p-value of about 0.00215. Additionally, looking at the p-values between specific coral reef types, we see that patch reef and back reef have no statistically significant difference, with a very high p-value and 0 contained in the confidence interval. However, nearshore is significantly different to back reef as is patch reef to nearshore. Both of these pairs have very low p-values and don’t contain 0 in their confidence intervals. Thus, based on the p-value from our anova test, we can reject the null hypothesis. However, it should be noted that the patch reefs and back reefs do not have the same apparent differences to each other as other reef types have.

  1. Load the penguins data from the palmerpenguins package. Do a two-way ANOVA to assess the effects of species and binary assigned sex at birth (aka the ‘sex’ column here) on beak length. This should be an interactive two-way ANOVA (Y~X*Z). Make a graph and run the stats. Assess the result using the following H0: “There is no effect of sex or species and there is no interaction between sex and species on penguin beak length”
#Graph goes here
means2 <- penguins %>%
  group_by(species, sex) %>% 
  summarize(mean = mean(bill_length_mm, na.rm=TRUE), st_error = std.error(bill_length_mm, na.rm=TRUE))
`summarise()` has grouped output by 'species'. You can override using the
`.groups` argument.
means2 <- na.omit(means2)
means2
# A tibble: 6 × 4
# Groups:   species [3]
  species   sex     mean st_error
  <fct>     <fct>  <dbl>    <dbl>
1 Adelie    female  37.3    0.237
2 Adelie    male    40.4    0.267
3 Chinstrap female  46.6    0.533
4 Chinstrap male    51.1    0.268
5 Gentoo    female  45.6    0.269
6 Gentoo    male    49.5    0.348
anov2_plot <- ggplot(means2, aes(x = species, y = mean, color = sex), na.rm = TRUE) + 
  geom_point(position = position_dodge(0.5))+ 
  geom_errorbar(aes(x=species, ymin=mean-st_error, ymax=mean+st_error), alpha=0.9, position = position_dodge(0.5)) + 
  labs(title="Effect of Penguin Species and Sex on Bill Length", y ="Bill Length (mm)") +
  scale_fill_brewer(palette="Set3") 
anov2_plot

What do we see in the graph?

We first see that there’s a big difference between the bill lengths of Adelie penguins compared to both chinstrap and gentoo penguins, on average. There is also a large difference in bill length between male and female penguins of each species type.

#stats go here!
anova2_test <- aov(bill_length_mm ~ species*sex, data = penguins)
summary(anova2_test)
             Df Sum Sq Mean Sq F value Pr(>F)    
species       2   7015    3508 654.189 <2e-16 ***
sex           1   1136    1136 211.807 <2e-16 ***
species:sex   2     24      12   2.284  0.103    
Residuals   327   1753       5                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
11 observations deleted due to missingness
TukeyHSD(anova2_test)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = bill_length_mm ~ species * sex, data = penguins)

$species
                      diff       lwr        upr     p adj
Chinstrap-Adelie 10.009851  9.209424 10.8102782 0.0000000
Gentoo-Adelie     8.744095  8.070779  9.4174106 0.0000000
Gentoo-Chinstrap -1.265756 -2.094535 -0.4369773 0.0010847

$sex
                diff      lwr      upr p adj
male-female 3.693368 3.194089 4.192646     0

$`species:sex`
                                     diff       lwr        upr     p adj
Chinstrap:female-Adelie:female   9.315995  7.937732 10.6942588 0.0000000
Gentoo:female-Adelie:female      8.306259  7.138639  9.4738788 0.0000000
Adelie:male-Adelie:female        3.132877  2.034137  4.2316165 0.0000000
Chinstrap:male-Adelie:female    13.836583 12.458320 15.2148470 0.0000000
Gentoo:male-Adelie:female       12.216236 11.064727 13.3677453 0.0000000
Gentoo:female-Chinstrap:female  -1.009736 -2.443514  0.4240412 0.3338130
Adelie:male-Chinstrap:female    -6.183118 -7.561382 -4.8048548 0.0000000
Chinstrap:male-Chinstrap:female  4.520588  2.910622  6.1305547 0.0000000
Gentoo:male-Chinstrap:female     2.900241  1.479553  4.3209291 0.0000002
Adelie:male-Gentoo:female       -5.173382 -6.341002 -4.0057622 0.0000000
Chinstrap:male-Gentoo:female     5.530325  4.096547  6.9641020 0.0000000
Gentoo:male-Gentoo:female        3.909977  2.692570  5.1273846 0.0000000
Chinstrap:male-Adelie:male      10.703707  9.325443 12.0819703 0.0000000
Gentoo:male-Adelie:male          9.083360  7.931851 10.2348686 0.0000000
Gentoo:male-Chinstrap:male      -1.620347 -3.041035 -0.1996591 0.0149963

Looking at stats and graph, what do we see?

Species has a very low p-value, indicating that the species type is associated with varying bill length. Sex also has a very low p-value, which means that the sex of the penguin is also associated with varying bill length. However, the p-value of the interaction between species and sex is too high, meaning that the connection between species and sex may not be significant to bill length.

Then looking at the Tukey test, the only interaction that doesn’t have a significant impact is the Gentoo:female-Chinstrap:female. Every other interaction has a significant impact on bill length.

Thus, based on the graph and statistical results we can reject the null hypothesis that species has no effect on bill length and similarly reject the nulls that there sex has no effect on bill length and that there is no interaction between sex and species on penguin beak length.

Part 3: Time series data!

Generally speaking, R has a tough time dealing with time formatting. This is a common issue in programming. There are special data formats called Datetime that we can use to work through this.

  1. Read in this time series data from github
temp_data <- read.csv("https://raw.githubusercontent.com/jbaumann3/Belize-RT-Baumann-et-al-2021/main/temperature%20and%20light/RT_in_situ_temp.csv")
head(temp_data)
        Date Month Day Year Time..GMT.   temp       site
1 12/11/2017    12  11 2017 3:00:00 PM 26.524 False Caye
2 12/11/2017    12  11 2017 3:15:00 PM 26.549 False Caye
3 12/11/2017    12  11 2017 3:30:00 PM 26.598 False Caye
4 12/11/2017    12  11 2017 3:45:00 PM 26.671 False Caye
5 12/11/2017    12  11 2017 4:00:00 PM 26.744 False Caye
6 12/11/2017    12  11 2017 4:15:00 PM 26.769 False Caye

These data are from an experiment I did recently on coral reefs in Belize. You should see temperature data from 2 sites across a range of dates.

  1. Make a new column the combines the date and time into 1 column, name that column datetime
 temp_data$datetime <- paste(temp_data$Date, temp_data$Time..GMT.)
head(temp_data)
        Date Month Day Year Time..GMT.   temp       site              datetime
1 12/11/2017    12  11 2017 3:00:00 PM 26.524 False Caye 12/11/2017 3:00:00 PM
2 12/11/2017    12  11 2017 3:15:00 PM 26.549 False Caye 12/11/2017 3:15:00 PM
3 12/11/2017    12  11 2017 3:30:00 PM 26.598 False Caye 12/11/2017 3:30:00 PM
4 12/11/2017    12  11 2017 3:45:00 PM 26.671 False Caye 12/11/2017 3:45:00 PM
5 12/11/2017    12  11 2017 4:00:00 PM 26.744 False Caye 12/11/2017 4:00:00 PM
6 12/11/2017    12  11 2017 4:15:00 PM 26.769 False Caye 12/11/2017 4:15:00 PM
  1. Using functions from the lubridate package (a lifesaver!), convert that column to datetime format
temp_data$datetime <- mdy_hms(temp_data$datetime)
head(temp_data)
        Date Month Day Year Time..GMT.   temp       site            datetime
1 12/11/2017    12  11 2017 3:00:00 PM 26.524 False Caye 2017-12-11 15:00:00
2 12/11/2017    12  11 2017 3:15:00 PM 26.549 False Caye 2017-12-11 15:15:00
3 12/11/2017    12  11 2017 3:30:00 PM 26.598 False Caye 2017-12-11 15:30:00
4 12/11/2017    12  11 2017 3:45:00 PM 26.671 False Caye 2017-12-11 15:45:00
5 12/11/2017    12  11 2017 4:00:00 PM 26.744 False Caye 2017-12-11 16:00:00
6 12/11/2017    12  11 2017 4:15:00 PM 26.769 False Caye 2017-12-11 16:15:00
  1. Make a time series plot of temperature that shows temp (or light) on y and datetime on x. Color by site
time_series <- ggplot(temp_data, aes(x = datetime, y=temp, color = site)) +   
  geom_line() + 
  geom_point() +
  labs(x = "datetime", y = "temperature")
time_series