Lab_3_Assignment

Author

Chii Kojima

#load packages
library(tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.4.0      ✔ purrr   1.0.1 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.3.0      ✔ stringr 1.5.0 
✔ readr   2.1.3      ✔ forcats 0.5.2 
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(lubridate)


Attaching package: 'lubridate'

The following objects are masked from 'package:base':

    date, intersect, setdiff, union

library(patchwork)

Essentials:

1.) Read in https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-05-04/water.csv check the format of the date column and change it using lubridate so it is correct

#read in the water data
water<-read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-05-04/water.csv')

Rows: 473293 Columns: 13
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (9): report_date, status_id, water_source, water_tech, facility_type, co...
dbl (4): row_id, lat_deg, lon_deg, install_year

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

str(water)

spc_tbl_ [473,293 × 13] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ row_id       : num [1:473293] 3957 33512 35125 37760 38118 ...
 $ lat_deg      : num [1:473293] 8.073 7.374 0.773 0.781 0.779 ...
 $ lon_deg      : num [1:473293] 38.6 40.5 34.9 35 35 ...
 $ report_date  : chr [1:473293] "04/06/2017" "08/04/2020" "03/18/2015" "03/18/2015" ...
 $ status_id    : chr [1:473293] "y" "y" "y" "y" ...
 $ water_source : chr [1:473293] NA "Protected Spring" "Protected Shallow Well" "Borehole" ...
 $ water_tech   : chr [1:473293] NA NA NA NA ...
 $ facility_type: chr [1:473293] NA "Improved" "Improved" "Improved" ...
 $ country_name : chr [1:473293] "Ethiopia" "Ethiopia" "Kenya" "Kenya" ...
 $ install_year : num [1:473293] NA 2019 NA NA NA ...
 $ installer    : chr [1:473293] "Private-CRS" "WaterAid" NA NA ...
 $ pay          : chr [1:473293] NA NA NA NA ...
 $ status       : chr [1:473293] NA NA NA NA ...
 - attr(*, "spec")=
  .. cols(
  ..   row_id = col_double(),
  ..   lat_deg = col_double(),
  ..   lon_deg = col_double(),
  ..   report_date = col_character(),
  ..   status_id = col_character(),
  ..   water_source = col_character(),
  ..   water_tech = col_character(),
  ..   facility_type = col_character(),
  ..   country_name = col_character(),
  ..   install_year = col_double(),
  ..   installer = col_character(),
  ..   pay = col_character(),
  ..   status = col_character()
  .. )
 - attr(*, "problems")=<externalptr>

water$report_date<-mdy_hms(water$report_date)#Converts the data column into data/time object

Warning: All formats failed to parse. No formats found.

str(water) #Data is no longer a factor but is now a POSIXct

spc_tbl_ [473,293 × 13] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ row_id       : num [1:473293] 3957 33512 35125 37760 38118 ...
 $ lat_deg      : num [1:473293] 8.073 7.374 0.773 0.781 0.779 ...
 $ lon_deg      : num [1:473293] 38.6 40.5 34.9 35 35 ...
 $ report_date  : POSIXct[1:473293], format: NA NA ...
 $ status_id    : chr [1:473293] "y" "y" "y" "y" ...
 $ water_source : chr [1:473293] NA "Protected Spring" "Protected Shallow Well" "Borehole" ...
 $ water_tech   : chr [1:473293] NA NA NA NA ...
 $ facility_type: chr [1:473293] NA "Improved" "Improved" "Improved" ...
 $ country_name : chr [1:473293] "Ethiopia" "Ethiopia" "Kenya" "Kenya" ...
 $ install_year : num [1:473293] NA 2019 NA NA NA ...
 $ installer    : chr [1:473293] "Private-CRS" "WaterAid" NA NA ...
 $ pay          : chr [1:473293] NA NA NA NA ...
 $ status       : chr [1:473293] NA NA NA NA ...
 - attr(*, "spec")=
  .. cols(
  ..   row_id = col_double(),
  ..   lat_deg = col_double(),
  ..   lon_deg = col_double(),
  ..   report_date = col_character(),
  ..   status_id = col_character(),
  ..   water_source = col_character(),
  ..   water_tech = col_character(),
  ..   facility_type = col_character(),
  ..   country_name = col_character(),
  ..   install_year = col_double(),
  ..   installer = col_character(),
  ..   pay = col_character(),
  ..   status = col_character()
  .. )
 - attr(*, "problems")=<externalptr>

2 Read in some data from Justin’s PhD!

bzcoral<-read.csv('https://raw.githubusercontent.com/jbaumann3/BIOL234_Biostats_MHC/main/coralcover.csv')
head(bzcoral)

  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

3 Look at the data and generate a hypothesis to test (X, Y, and/or Z has no effect on coral cover (cc_percent), for example). cc_percent is the only numerical variable we care about here! Everything else is categorical. Write your hypothesis in BOLD below.

Null hypothesis: The percentage of Coral Coverage is not different in different sites

Alternate hypothesis: The percentage of Coral Coverage is different in different sites

4 Filter our the columns that you are not using for your hypothesis test

#Select and remove columns/rows
bzcoral2<-select(bzcoral,site, cc_percent)
head(bzcoral2)

  site cc_percent
1    1  5.8435758
2    1  0.9505263
3    1  5.2423389
4    1  5.0040475
5    1  5.8954916
6    2  5.2826190

5 Using the pipe, %>%, group your data and calcualte the mean(s) you need for your visual hypothesis test

bzcoral3<-bzcoral %>%
  select(site, cc_percent) 

  head(bzcoral3)

  site cc_percent
1    1  5.8435758
2    1  0.9505263
3    1  5.2423389
4    1  5.0040475
5    1  5.8954916
6    2  5.2826190

6 Graph your results! Means + errorbars required :) Make a nice, easy to see graph with clear labels and text

meanbzcoral<-bzcoral3%>%
  group_by(site)%>% #group = site
  summarize(mean=mean(cc_percent),sd=sd(cc_percent),n=n(),se=sd/sqrt(n))

meanbzcoral

# A tibble: 13 × 5
    site  mean     sd     n    se
   <int> <dbl>  <dbl> <int> <dbl>
 1     1  4.59  2.07      5 0.925
 2     2 11.7   8.16      6 3.33 
 3     3  2.69  2.94      6 1.20 
 4     4 17.7   7.28      6 2.97 
 5     5  7.88  3.09      6 1.26 
 6     6  3.99  2.22      6 0.907
 7     7 25.1   4.92      6 2.01 
 8     8 28.6  12.3       6 5.03 
 9     9  2.79  1.56      6 0.639
10    11 14.2   7.35      6 3.00 
11    12  6.86  3.57      6 1.46 
12    13  4.08  2.54      6 1.04 
13    14  1.33  0.800     6 0.327

meanbzcoral$site=as.factor(meanbzcoral$site) #Make site as factor

ggplot(data=meanbzcoral,aes(x=site,y=mean,color=site))+
  geom_point()+
  geom_errorbar(data=meanbzcoral,aes(x=site,ymin=mean-se,ymax=mean+se))

7 Assess your hypothesis! What does your graph show (Note: We did not do stats, so please do not say ‘significant’)

Since the mean percentage of coral coverage is different in different sites, we reject our null hypothesis. Therefore the percentage of coral coverage is different in different sites

Depth

1: Read in intertidal transect data. View it, identify the columns that contain species/cover items and pivot from wide to long format!

op<-read.csv('https://raw.githubusercontent.com/jbaumann3/BIOL234_Biostats_MHC/main/Ocean_Point_Data%202018.csv')
head(op)

  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

longop<-op %>%
  pivot_longer(bare_rock:Cancer.crabs,names_to='species',values_to='number')

longop

# A tibble: 4,495 × 7
   Number trans_description   year tidal_ht_m tidal_ht_cat_num species    number
    <int> <chr>              <int>      <dbl>            <int> <chr>       <dbl>
 1      0 a_cobble_protected  2017      -0.13                1 bare_rock       0
 2      0 a_cobble_protected  2017      -0.13                1 calothr         0
 3      0 a_cobble_protected  2017      -0.13                1 semibal         0
 4      0 a_cobble_protected  2017      -0.13                1 mytilus         0
 5      0 a_cobble_protected  2017      -0.13                1 Halichond…      0
 6      0 a_cobble_protected  2017      -0.13                1 Spongomor…      0
 7      0 a_cobble_protected  2017      -0.13                1 chaetomo        0
 8      0 a_cobble_protected  2017      -0.13                1 ulva_lac        0
 9      0 a_cobble_protected  2017      -0.13                1 ascophyl        0
10      0 a_cobble_protected  2017      -0.13                1 lamin_di        0
# … with 4,485 more rows

2 Filter data such that we retain only ‘species’ that are animals. These include: Carinbus, Cancer Crabs, nucella, litt_obt, litt_lit, litt_sax. Everything else is either rock or algae. Please do this filter step AFTER you pivot to long format. Note: It is easier to do this when you are in wide format, but this is good practice! Ask questions if you need help :). Hint: The ‘|’ key is how you say ‘or’ in code :)

#Select and remove columns/rows
opspecies<-filter(longop,species=="semibal"|species=="Cancer.crabs"|species=="nucella"|species=="carcinus"|species=="litt_lit"|species=="litt_sax"|species=="litt_obt")

opspecies

# A tibble: 1,085 × 7
   Number trans_description   year tidal_ht_m tidal_ht_cat_num species    number
    <int> <chr>              <int>      <dbl>            <int> <chr>       <dbl>
 1      0 a_cobble_protected  2017      -0.13                1 semibal         0
 2      0 a_cobble_protected  2017      -0.13                1 litt_sax        0
 3      0 a_cobble_protected  2017      -0.13                1 litt_lit       64
 4      0 a_cobble_protected  2017      -0.13                1 litt_obt        0
 5      0 a_cobble_protected  2017      -0.13                1 nucella         0
 6      0 a_cobble_protected  2017      -0.13                1 carcinus        0
 7      0 a_cobble_protected  2017      -0.13                1 Cancer.cr…      0
 8      0 a_cobble_protected  2017      -0.13                1 semibal         0
 9      0 a_cobble_protected  2017      -0.13                1 litt_sax        0
10      0 a_cobble_protected  2017      -0.13                1 litt_lit        0
# … with 1,075 more rows

2 Using the same data–> rename the factors in the trans_description column based on wave exposure. a_cobble_protected is low, c_flat_profile and d_semi_expos are both moderate, e_exposed is high. Hint: use ifelse(). Justin can help!

opspecies$wave<-ifelse(opspecies$trans_description=="a_cobble_protected","low",
                       ifelse(opspecies$trans_description==" e_exposed","high","moderate"))
opspecies

# A tibble: 1,085 × 8
   Number trans_description   year tidal_ht_m tidal_ht_ca…¹ species number wave 
    <int> <chr>              <int>      <dbl>         <int> <chr>    <dbl> <chr>
 1      0 a_cobble_protected  2017      -0.13             1 semibal      0 low  
 2      0 a_cobble_protected  2017      -0.13             1 litt_s…      0 low  
 3      0 a_cobble_protected  2017      -0.13             1 litt_l…     64 low  
 4      0 a_cobble_protected  2017      -0.13             1 litt_o…      0 low  
 5      0 a_cobble_protected  2017      -0.13             1 nucella      0 low  
 6      0 a_cobble_protected  2017      -0.13             1 carcin…      0 low  
 7      0 a_cobble_protected  2017      -0.13             1 Cancer…      0 low  
 8      0 a_cobble_protected  2017      -0.13             1 semibal      0 low  
 9      0 a_cobble_protected  2017      -0.13             1 litt_s…      0 low  
10      0 a_cobble_protected  2017      -0.13             1 litt_l…      0 low  
# … with 1,075 more rows, and abbreviated variable name ¹tidal_ht_cat_num

3 Using the same data, make a simplified tidal height column. We want tidal height cat <4 to be ‘low’, >7 to be m’oderate’high’, and in between to be ‘intermediate’. Hint: You can use if_else for this as well!

opspecies$tide<-ifelse(opspecies$tidal_ht_cat_num<4,"low",
                       ifelse(opspecies$tidal_ht_cat_num>7,"high", "moderate"))
opspecies

# A tibble: 1,085 × 9
   Number trans_description   year tidal_ht_m tidal…¹ species number wave  tide 
    <int> <chr>              <int>      <dbl>   <int> <chr>    <dbl> <chr> <chr>
 1      0 a_cobble_protected  2017      -0.13       1 semibal      0 low   low  
 2      0 a_cobble_protected  2017      -0.13       1 litt_s…      0 low   low  
 3      0 a_cobble_protected  2017      -0.13       1 litt_l…     64 low   low  
 4      0 a_cobble_protected  2017      -0.13       1 litt_o…      0 low   low  
 5      0 a_cobble_protected  2017      -0.13       1 nucella      0 low   low  
 6      0 a_cobble_protected  2017      -0.13       1 carcin…      0 low   low  
 7      0 a_cobble_protected  2017      -0.13       1 Cancer…      0 low   low  
 8      0 a_cobble_protected  2017      -0.13       1 semibal      0 low   low  
 9      0 a_cobble_protected  2017      -0.13       1 litt_s…      0 low   low  
10      0 a_cobble_protected  2017      -0.13       1 litt_l…      0 low   low  
# … with 1,075 more rows, and abbreviated variable name ¹tidal_ht_cat_num

4 Using the same dataframe, build 2 new dataframes, one calculating mean and error (standard error) for semibal (barnacle) abundance and one doing the same for nucella (whelk) abundance by tidal height and wave exposure (trans_description)

semibal<-filter(opspecies,species=='semibal') %>%
  group_by(wave,tide)%>%
  summarize(mean=mean(number),sd=sd(number),n=n(),se=sd/sqrt(n))

`summarise()` has grouped output by 'wave'. You can override using the
`.groups` argument.

semibal

# A tibble: 6 × 6
# Groups:   wave [2]
  wave     tide      mean    sd     n    se
  <chr>    <chr>    <dbl> <dbl> <int> <dbl>
1 low      high      0     0        9 0    
2 low      low       6.11  9.53     9 3.18 
3 low      moderate  2.24  5.24    17 1.27 
4 moderate high      0     0       30 0    
5 moderate low       1.97  4.63    36 0.771
6 moderate moderate 20.4  25.6     54 3.49

nucella<-filter(opspecies,species=='nucella') %>%
　group_by(wave,tide)%>%
  summarize(mean=mean(number),sd=sd(number),n=n(),se=sd/sqrt(n))

`summarise()` has grouped output by 'wave'. You can override using the
`.groups` argument.

nucella

# A tibble: 6 × 6
# Groups:   wave [2]
  wave     tide      mean    sd     n    se
  <chr>    <chr>    <dbl> <dbl> <int> <dbl>
1 low      high      0     0        9  0   
2 low      low       0     0        9  0   
3 low      moderate  0     0       17  0   
4 moderate high      0     0       30  0   
5 moderate low       2.67  8.11    36  1.35
6 moderate moderate  2.96 10.8     54  1.47

5 Plot barnacle and snail abundance + error across tidal heights & wave exposures. You can plot them both on the same graphs (merging your dataframes could help) or you can make multiple graphs and patchwork them together.

semibal

# A tibble: 6 × 6
# Groups:   wave [2]
  wave     tide      mean    sd     n    se
  <chr>    <chr>    <dbl> <dbl> <int> <dbl>
1 low      high      0     0        9 0    
2 low      low       6.11  9.53     9 3.18 
3 low      moderate  2.24  5.24    17 1.27 
4 moderate high      0     0       30 0    
5 moderate low       1.97  4.63    36 0.771
6 moderate moderate 20.4  25.6     54 3.49

p1=ggplot(data=semibal,aes(x=wave,y=mean,color=tide))+
  geom_point()+
  geom_errorbar(data=semibal,aes(x=wave,ymin=mean-se,ymax=mean+se))+
  labs(title="Semibal")

p2=ggplot(data=nucella,aes(x=wave,y=mean,color=tide))+
  geom_point()+
  geom_errorbar(data=nucella,aes(x=wave,ymin=mean-se,ymax=mean+se))+
  labs(title="Nucella")

p1+p2

6 Write a short interpretative statement. We didn’t run any stats, so avoid the word ‘significant.’ How do snail and whelk abundances appear to vary by tidal height and wave exposure? I do not need you to tell me anything about the ecology of the system, but feel free to do so if you’d like :) I just need you to interpret your graphs. ( nucella (whelk))

For snail, in low wave condition there are more abundance in the cobination of low tide, but in moderate wave the most abundance was found in the combination of moderate wave. For whelk, the moderate wave and moderate tide seems the best condition for them to survive.