Environmental Risk

Part 2

Planetary Boundaries

Coral Bleaching

Longterm Monitoring Data

Collected by CRIOBE since 2005
Point intercept transects
3 transects per site

Longterm Monitoring Data

Task 2.1

Explore the coral_cover set:

A) Which years are covered?

B) Are there any gaps?

C) How many sites?

D) Same number of transects per site and year?

year	date	site	transect	percent
2005	2005-02-16	Haapiti	1	32
2005	2005-02-16	Haapiti	2	18
2005	2005-02-16	Haapiti	3	36
2005	2005-02-16	Taotaha	1	32
2005	2005-02-16	Taotaha	2	42
2005	2005-02-16	Taotaha	3	38
2005	2005-02-17	Tetaiuo	1	44
2005	2005-02-17	Tetaiuo	2	32
2005	2005-02-17	Tetaiuo	3	44
2005	2005-02-18	Entre 2 baies	1	70
2005	2005-02-18	Entre 2 baies	2	50

Longterm Monitoring Data

Task 2.2

Make a similar plot

The points are the mean percent values per site
The line is the mean of these mean values
The shaded area in this mean ± the standard deviation (sd())

Bleaching event 2019

Bleaching event 2019 - Sites

Bleaching event 2019

Task 2.3

Focus on data from 2019 until 2023
Explore differences between sites
How can you include the additional information?
- Lines on top of each other to show variability
- Different panels for each site to see what is happening where
- Other ideas?

Bleaching event 2019

Why is the change in cover different between sites?

Degree Heating Weeks

Measure for accumulated heat stress
Satellite derived data
Based on Sea Surface Temperature (SST)

Calculation

Calculate difference between SST and bleaching threshold: long-term mean (MMM) + 1°C
Sum up all differences > 0 for the last 12 weeks (84 days)
Divide by 7 (=> unit is °C weeks)
Continue with next day

\[ \textrm{DHW} = \frac{1}{7}\sum_{i}^{84} BT\textrm{, where }BT \geq 1 \]

Degree Heating Weeks

Bleaching 2019

Strongest bleaching event in 30 years

Analysis: Impact of DHW on coral cover

Plan

Calculate the change in coral cover in 2019 (before bleaching) and 2020 (after bleaching)
Add the DWH_max data for each site
Run a regression on the change in coral cover and DWH_max
Visualize the results

1. Change in cover

Aim

Calculate the change in coral cover in 2019 (before bleaching) and 2020 (after bleaching)

Considerations

The transect were not done at exactly the same spot => Transect 1 in 2020 not exactly at same spot as in 2019

Consequence

Calculate the mean cover per site and year before calculating the difference

1. Change in cover

Aim

Calculate the change in coral cover in 2019 (before bleaching) and 2020 (after bleaching)

Considerations

The coral cover in 2019 differed between the sites, leading to different “start values”. This affects how much the cover can be reduced.

Consequence

Calculate the relative change in cover

dat_change_coral_cover <- dat_coral_cover %>% 
  filter(year >= 2019, year <= 2020) %>% 
  group_by(site, year) %>% 
  summarise(percent = mean(percent)) %>% 
  pivot_wider(names_from = year, 
              values_from = percent) %>% 
  mutate(change = `2020` - `2019`) %>% 
  select(-`2020`, -`2019`)

1. Change in cover

Aim

Calculate the change in coral cover in 2019 (before bleaching) and 2020 (after bleaching)

Considerations

The coral cover in 2019 differed between the sites, leading to different “start values”. This affects how much the cover can be reduced.

Consequence

Calculate the relative change in cover

dat_change_coral_cover <- dat_coral_cover %>% 
  filter(year >= 2019, year <= 2020) %>% 
  group_by(site, year) %>% 
  summarise(percent = mean(percent)) %>% 
  pivot_wider(names_from = year, 
              values_from = percent) %>% 
  # calculate relative change
  mutate(rel_change = 100*(`2019` - `2020`)/`2019`) %>% 
  select(-`2020`, -`2019`)

2. DWH_max

Task 2.4

Use or left_join() to combine dat_change_coral_cover and the DWH_max data based on site.

dat_change_coral_cover

site	rel_change
Afareaitu	38.98305
Aroa	41.58416
Entre 2 baies	69.74790
Gendron	75.58140
Haapiti	64.61538
Maatea	26.53061
Motu Ahi	38.88889

DWH_max data

site	max_dhw
Haapiti	6.212857
Taotaha	6.190429
Tetaiuo	6.229857
Gendron	6.264714
Tiahura	6.240286
Entre 2 baies	6.211000
Pihaena	6.164286

2. DWH_max

3. Regression

Type of linear model
Tries to explain how one variable (e.g. height) influences another (e.g. weight)
Can be used to
- predict e.g. weight only with height
- assess importance (significant or not)

3. Regression

Finds the best fitting line through a point cloud

\[ y_i = a + b x_i + \epsilon_i \]

\(a\) is the y-intercept, \(b\) is the slope
\(y_i\) is the expected y value at \(x_i\)
The data does not fit perfectly to the line, each actual data point can vary by \(\epsilon\)

3. Regression

Example

Here, \(a\) = 3.41 and \(b\) = 0.38

What is the expected weight for a height of 166 cm?

3. Regression

Example

Here, \(a\) = 3.41 and \(b\) = 0.38
Expected weight for a height of 166 cm is 66.5 kg

3. Regression

Example

Here, \(a\) = 3.41 and \(b\) = 0.38
Expected weight for a height of 166 cm is 66.5 kg
In the actual data, the corresponding weight is 64.4 kg

3. Regression

Example

Here, \(a\) = 3.41 and \(b\) = 0.38
Expected weight for a height of 166 cm is 66.5 kg
In the actual data, the corresponding weight is 64.4 kg
At this point, the error \(\epsilon_{166~cm}\) is -2.1 kg

3. Regression

Residuals

These errors are called residuals
Over all the data points, the residuals are expected to follow a normal distribution with mean 0 (to allow positive and negative values) and standard deviation (\(\sigma\)):

\(\epsilon \sim normal(0, \sigma)\)

The standard deviation (\(\sigma\)) is calculated by the linear model function

3. Regression

Residuals

3. Regression

In R, it is done with the lm() function:

m_weight <- lm(weight ~ height, 
            data = dat_weight)

3. Regression

In R, it is done with the lm() function:

m_weight <- lm(weight ~ height, 
            data = dat_weight)

summary(m_weight)


Call:
lm(formula = weight ~ height, data = dat_weight)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0941 -1.7353 -0.0545  1.0740  4.8862 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.41513    7.22304   0.473    0.638    
height       0.38023    0.04304   8.835 1.24e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.019 on 48 degrees of freedom
Multiple R-squared:  0.6192,    Adjusted R-squared:  0.6113 
F-statistic: 78.06 on 1 and 48 DF,  p-value: 1.242e-11

3. Regression

In R, it is done with the lm() function:

m_weight <- lm(weight ~ height, 
            data = dat_weight)

summary(m_weight)


Call:
lm(formula = weight ~ height, data = dat_weight)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0941 -1.7353 -0.0545  1.0740  4.8862 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.41513    7.22304   0.473    0.638    
height       0.38023    0.04304   8.835 1.24e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.019 on 48 degrees of freedom
Multiple R-squared:  0.6192,    Adjusted R-squared:  0.6113 
F-statistic: 78.06 on 1 and 48 DF,  p-value: 1.242e-11

Translation

The intercept is 3.41513 and not significantly different to 0 (we don’t care)
The slope is 0.38023 and significantly different to 0
There is a significant effect of height on weight
The Adjusted R² is 0.6113, meaning that our model explains about 61% of the variance in the data

3. Regression

Task 2.5

Back to the coral cover data set

Analyse the impact of max_dhw on rel_change

Is this effect significant?
What is the expected change in coral cover per unit DHW?

3. Regression


Call:
lm(formula = rel_change ~ max_dhw, data = dat_change_coral_cover)

Residuals:
     Min       1Q   Median       3Q      Max 
-21.2426  -7.4508  -0.2689   8.6999  15.8141 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -501.16     109.01  -4.597 0.000768 ***
max_dhw        90.52      17.97   5.036 0.000380 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11.84 on 11 degrees of freedom
Multiple R-squared:  0.6975,    Adjusted R-squared:   0.67 
F-statistic: 25.36 on 1 and 11 DF,  p-value: 0.0003802

Meaning

The intensity of the heat stress, measured as the maximum degree heating weeks in 2019 (DHW_max), significantly impacted the decline in coral cover (t(11) = 5.036, p < .001, adj. R² = 0.67). The expected relative decline in coral cover increased by 90.5% ± 18.0% (SE) per unit of DHW_max.

4. Visualisation

Weight data set

Layers

Raw Data

4. Visualisation

Weight data set

Layers

Raw Data
Linear equation:

\(\Delta\) Relative cover = -501.2 + 90.5 * DHW_max

4. Visualisation

Weight data set

Layers

Raw Data
Linear equation:

\(\Delta\) Relative cover = -501.2 + 90.5 * DHW_max

Standard error of predictions

4. Visualisation

Use model to predict data

Steps

Create a new data frame with a sequence of x values for which you want to predict y values.

Code

ndat_weight <-  data.frame(height = seq(min(dat_weight$height),
                                        max(dat_weight$height),
                                        length = 100))

4. Visualisation

Use model to predict data

Steps

Create a new data frame with a sequence of x values for which you want to predict y values.
Predict y values for x values in new data frame

Code

pred_weight <- predict(m_weight,
                       newdata = ndat_weight,
                       se.fit = T) %>% 
  bind_cols(ndat_weight) # add height to predictions

4. Visualisation

Use model to predict data

Steps

Create a new data frame with a sequence of x values for which you want to predict y values.
Predict y values for x values in new data frame
Plot model as line, SE as ribbon, and raw data in background

Code

pred_weight %>% 
  # only height is the same in both data sets
  ggplot(aes(x = height))+
  # Plot SE
  geom_ribbon(aes(ymin = fit - se.fit,
                  ymax = fit + se.fit),
              fill = "grey")+
  #Plot model line
  geom_line(aes(y = fit), 
            col = "darkred", linewidth = 1.2)+
  # Plot raw data
  geom_point(data = dat_weight,
             aes(y = weight), size = 4)+
  # Formatting
  labs(x = "Height (cm)", y = "Weight (kg)")+
  theme_light(base_size = 28)+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

4. Visualisation

Task 2.6

Plot the coral cover data and model in a similar way

Thermotolerance

Many studies show that the thermotlerance differs between genera and species!

Thermotolerance

Download this RStudio project
Double-click on bleaching_anaylsis.Rproj
Have a look at scripts/1_regression_coralcover_dhwmax.qmd. It contains the current analysis for coral cover.
Open scripts/2_regression_coralcover_genus_dhwmax.qmd. I prepared the analysis for the three main coral genera. I indicated the analysis steps in the scripts. Continue on your own

Thermotolerance

Tasks

For each genus, assess the impact of DHW_max on the relative coral decline with a linear model
Visualize the results
Interpret the results. Does the impact of heat stress differs between genera? What does this mean?

Thermotolerance

Statistical Rethinking

Data analysis with Bayesian methods

Course

Accompanying course with slides and videos

Environmental Risk

Planetary Boundaries

Coral Bleaching

Coral Bleaching

Longterm Monitoring Data

Longterm Monitoring Data

Longterm Monitoring Data

Bleaching event 2019

Bleaching event 2019

Bleaching event 2019

Bleaching event 2019 - Sites

Bleaching event 2019

Bleaching event 2019

Bleaching event 2019

Degree Heating Weeks

Degree Heating Weeks

Degree Heating Weeks

Degree Heating Weeks

Degree Heating Weeks

Degree Heating Weeks

Bleaching 2019

Analysis: Impact of DHW on coral cover

Plan

1. Change in cover

1. Change in cover

1. Change in cover

2. DWHmax

2. DWHmax

3. Regression

3. Regression

3. Regression

3. Regression

3. Regression

3. Regression

3. Regression

3. Regression

3. Regression

3. Regression

3. Regression

Translation

3. Regression

3. Regression

Meaning

4. Visualisation

Weight data set

Layers

4. Visualisation

Weight data set

Layers

4. Visualisation

Weight data set

Layers

4. Visualisation

Use model to predict data

4. Visualisation

Use model to predict data

4. Visualisation

Use model to predict data

4. Visualisation

Thermotolerance

Thermotolerance

Thermotolerance

Thermotolerance

Read more

2. DWH_max

2. DWH_max