Climate change and temperature anomalies

weather <- 
  read_csv("https://data.giss.nasa.gov/gistemp/tabledata_v4/NH.Ts+dSST.csv", 
           skip = 1, 
           na = "***")
# SELECT YEAR, MONTHS AND PIVOT LONGER
weather %>% 
  select(1:13) %>% 
  pivot_longer(2:13, names_to = "month", values_to = "delta") -> tidyweather

Plotting Information

tidyweather <- tidyweather %>%
  mutate(date = ymd(paste(as.character(Year), month, "1")),
         month = month(date, label=TRUE),
         Year = year(date))

ggplot(tidyweather, aes(x=date, y = delta))+
  geom_point()+
  geom_smooth(color="red") +
  theme_bw() +
  labs (
    title = "Weather Anomalies",
    x = "Year",
    y = "Delta"
  )

ggplot(tidyweather, aes(x=date, y = delta))+
  geom_point()+
  geom_smooth(color="red") +
  facet_wrap(~month) +
  theme_bw() +
  labs (
    title = "Weather Anomalies",
    x = "Year",
    y = "Delta"
  )

There is no apparent difference between the effects of increasing temperature in monthly data. January and February seem to represent the strongest change, but it does not vary significantly from other months.

comparison <- tidyweather %>% 
  filter(Year>= 1881) %>%     #remove years prior to 1881
  #create new variable 'interval', and assign values based on criteria below:
  mutate(interval = case_when(
    Year %in% c(1881:1920) ~ "1881-1920",
    Year %in% c(1921:1950) ~ "1921-1950",
    Year %in% c(1951:1980) ~ "1951-1980",
    Year %in% c(1981:2010) ~ "1981-2010",
    TRUE ~ "2011-present"
  ))
comparison %>% 
  ggplot(aes(x = delta, fill = interval, color = interval)) +
  geom_density(alpha = 0.3) +
  labs (
    title = "Weather Delta by Decade",
    x = "Delta",
    y = "Density",
    fill = "Decade",
    color = "Decade"
  )

#creating yearly averages
average_annual_anomaly <- tidyweather %>% 
  group_by(Year) %>%   #grouping data by Year
  # creating summaries for mean delta 
  # use `na.rm=TRUE` to eliminate NA (not available) values 
  summarise(mean_delta = mean(delta, na.rm=TRUE)) 

average_annual_anomaly %>% 
  ggplot(aes(x=Year, y=mean_delta)) +
  geom_point(aes(color=mean_delta>0)) + # DRAWING POINTS ABOVE ZERO A DIFFERENT COLOUR
  geom_smooth(method = "loess", color="black") +
  theme_bw() + 
  labs(
    title = "Weather Delta average by Year",
    x = "Year",
    y = "Mean Delta", 
    color = "Above Zero"
  )

Confidence Interval for delta

CI using formula

formula_ci <- comparison %>% 

  # choose the interval 2011-present
  filter(Year >= 2011) %>% 
  group_by(Year) %>% 
  summarise(
    mean_delta = mean(delta),
    sd_delta = sd(delta),
    count = n(),
    # We're choosing a 95% confidence interval:
    t_critical = qt(0.975, count-1),
    se_delta = sd(delta/sqrt(count)),
    margin_of_error = t_critical*se_delta,
    delta_low  = mean_delta - margin_of_error,
    delta_high = mean_delta + margin_of_error
  )
  
  # calculate summary statistics for temperature deviation (delta) 
  # calculate mean, SD, count, SE, lower/upper 95% CI
  # what dplyr verb will you use? 

#print out formula_CI
formula_ci
## # A tibble: 12 × 9
##     Year mean_delta sd_delta count t_critical se_delta margin_…¹ delta…² delta…³
##    <dbl>      <dbl>    <dbl> <int>      <dbl>    <dbl>     <dbl>   <dbl>   <dbl>
##  1  2011      0.745    0.113    12       2.20   0.0327    0.0720   0.673   0.817
##  2  2012      0.815    0.179    12       2.20   0.0517    0.114    0.701   0.929
##  3  2013      0.8      0.118    12       2.20   0.0340    0.0749   0.725   0.875
##  4  2014      0.92     0.145    12       2.20   0.0420    0.0924   0.828   1.01 
##  5  2015      1.18     0.178    12       2.20   0.0515    0.113    1.06    1.29 
##  6  2016      1.31     0.333    12       2.20   0.0961    0.212    1.10    1.52 
##  7  2017      1.18     0.226    12       2.20   0.0653    0.144    1.03    1.32 
##  8  2018      1.04     0.137    12       2.20   0.0396    0.0871   0.950   1.12 
##  9  2019      1.21     0.153    12       2.20   0.0441    0.0970   1.11    1.31 
## 10  2020      1.35     0.225    12       2.20   0.0648    0.143    1.21    1.50 
## 11  2021      1.14     0.117    12       2.20   0.0339    0.0746   1.06    1.21 
## 12  2022     NA       NA        12       2.20  NA        NA       NA      NA    
## # … with abbreviated variable names ¹​margin_of_error, ²​delta_low, ³​delta_high

In order to calculate the confidence intervals, we have calculated for every year many summary statistics including mean, SD, SE and sample size. Using these summary statistics we have calculated the t-Student distribution critical value, and multiplied the t_critical by the standard error to get the final margin of error. Finally we add and subtract the margin of error to the mean to get the limits of the interval. In order to fully understand the data we have collected, we have come up with the following visualization of the mean and confidence intervals:

formula_ci %>% 
  mutate(Year = as.factor(Year)) %>% 
  na.omit() %>% 
  ggplot(aes(color = Year)) +
  geom_pointrange(aes(x=Year, y=mean_delta, ymax=delta_high, ymin=delta_low)) +
  theme_bw() + 
  theme(legend.position = "none") +
  labs(
    title = "Weather Delta range by Year",
    subtitle = "NASA Weater Data",
    x = "Year",
    y = "Delta range", 
    color = NULL
  )

The visualization shows an increase in the mean over the years, reaching a first maximum in 2016, and an even higher maximum mean in 2020. For years with a higher mean delta, the confidence intervals seem to increase as well. With the current samples and intervals, we cannot say for sure that the mean delta in year 2016 was higher than the mean delta in 2017 for instance, but we can be sure with 95% confidence that 2011 was lower than 2016.

CI using formula (for all years)

Since there are not enough data points every year, the confidence intervals are very large and vary a lot by year. The total confidence interval for the entire period is as follows:

formula_ci_interval <- comparison %>% 

  # choose the interval 2011-present
  filter(Year >= 2011) %>% 
  na.omit() %>% 
  group_by(interval) %>% 
  summarise(
    mean_delta = mean(delta),
    sd_delta = sd(delta),
    count = n(),
    # We're choosing a 95% confidence interval:
    t_critical = qt(0.975, count-1),
    se_delta = sd(delta/sqrt(count)),
    margin_of_error = t_critical*se_delta,
    delta_low  = mean_delta - margin_of_error,
    delta_high = mean_delta + margin_of_error
  )

#print out formula_CI
formula_ci_interval %>% select(delta_low, delta_high)
## # A tibble: 1 × 2
##   delta_low delta_high
##       <dbl>      <dbl>
## 1      1.02       1.11

CI using bootstrapping

boot_dist <- comparison %>%
  # choose the interval 2011-present
  filter(Year >= 2011) %>%
  mutate(Year = as.factor(Year)) %>%
  specify(response=delta) %>%
  generate(reps=1000, type="bootstrap") %>% 
  calculate(stat = "mean")

boot_dist %>%
  # Calculate the confidence interval around the point estimate
  get_confidence_interval(
    # At the 95% confidence level; percentile method
    level = 0.95
  )
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1     1.03     1.11