Temperature Regression Analysis

Temps <- read_excel("Temp Data.xlsx", 
  col_types = c("numeric", "numeric", "numeric", 
                "numeric", "numeric", "numeric", 
                "numeric", "numeric", "numeric", 
                "numeric", "numeric", "numeric", 
                "numeric", "numeric", "numeric", 
                "numeric", "numeric", "numeric", 
                "numeric"))

## Warning: Expecting numeric in O2 / R2C15: got '***'

## Warning: Expecting numeric in P2 / R2C16: got '***'

## Warning: Expecting numeric in E148 / R148C5: got '***'

## Warning: Expecting numeric in F148 / R148C6: got '***'

## Warning: Expecting numeric in G148 / R148C7: got '***'

## Warning: Expecting numeric in H148 / R148C8: got '***'

## Warning: Expecting numeric in I148 / R148C9: got '***'

## Warning: Expecting numeric in J148 / R148C10: got '***'

## Warning: Expecting numeric in K148 / R148C11: got '***'

## Warning: Expecting numeric in L148 / R148C12: got '***'

## Warning: Expecting numeric in M148 / R148C13: got '***'

## Warning: Expecting numeric in N148 / R148C14: got '***'

## Warning: Expecting numeric in O148 / R148C15: got '***'

## Warning: Expecting numeric in Q148 / R148C17: got '***'

## Warning: Expecting numeric in R148 / R148C18: got '***'

## Warning: Expecting numeric in S148 / R148C19: got '***'

Temps <- Temps %>%
  filter(Year < 2026)
Temps <- Temps %>%
  mutate(across(where(is.numeric), ~ round(., 2)))

annual_model <- lm(`J-D` ~ Year, data = Temps)
summary(annual_model)

## 
## Call:
## lm(formula = `J-D` ~ Year, data = Temps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37706 -0.13891 -0.03332  0.12425  0.60653 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.607e+01  7.481e-01  -21.48   <2e-16 ***
## Year         8.273e-03  3.831e-04   21.60   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1951 on 144 degrees of freedom
## Multiple R-squared:  0.7641, Adjusted R-squared:  0.7625 
## F-statistic: 466.5 on 1 and 144 DF,  p-value: < 2.2e-16

ggplot(Temps, aes(x = Year, y = `J-D`)) +
  geom_point(color = "steelblue", alpha = 0.6) +
    labs(
    title = "Global Temperature Anomaly Over Time",
    x = "Year",
    y = "Temperature Anomaly (C)",
    caption = "Source: NASA GISS"
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold"),
    axis.text = element_text(face = "bold"),
    plot.title = element_text(face = "bold", size = 16)
  )

Adding Trend line to the data

ggplot(Temps, aes(x = Year, y = `J-D`)) +
  geom_point(color = "steelblue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  labs(
    title = "Global Temperature Anomaly Over Time",
    x = "Year",
    y = "Temperature Anomaly (C)",
    caption = "Source: NASA GISS"
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold"),
    axis.text = element_text(face = "bold"),
    plot.title = element_text(face = "bold", size = 16)
  )

## `geom_smooth()` using formula = 'y ~ x'

Now we can look at the trend lines per season

Temps_seasonal <- Temps %>%
  select(Year, Winter, Spring, Summer, Fall) %>%
  pivot_longer(cols = c(Winter, Spring, Summer, Fall),
               names_to = "Season",
               values_to = "Anomaly")

# Plot all four seasons with trend lines
ggplot(Temps_seasonal, aes(x = Year, y = Anomaly, color = Season)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1) +
  labs(
    title = "Seasonal Temperature Anomaly Trends Over Time",
    x = "Year",
    y = "Temperature Anomaly (C)",
    caption = "Source: NASA GISS"
  ) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold"),
    axis.text = element_text(face = "bold"),
    plot.title = element_text(face = "bold", size = 16)
  )

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

Summary of Seasonal trends

# Run regression for each season and extract key stats
seasonal_models <- Temps_seasonal %>%
  group_by(Season) %>%
  summarise(
    Slope = coef(lm(Anomaly ~ Year))[2],
    R_squared = summary(lm(Anomaly ~ Year))$r.squared,
    P_value = summary(lm(Anomaly ~ Year))$coefficients[2,4]
  ) %>%
  mutate(
    Slope = round(Slope, 5),
    R_squared = round(R_squared, 3),
    P_value = round(P_value, 6),
    Warming_per_decade = round(Slope * 10, 3)
  ) %>%
  arrange(desc(Slope))

print(seasonal_models)

## # A tibble: 4 × 5
##   Season   Slope R_squared P_value Warming_per_decade
##   <chr>    <dbl>     <dbl>   <dbl>              <dbl>
## 1 Winter 0.00868     0.725       0              0.087
## 2 Spring 0.00852     0.741       0              0.085
## 3 Fall   0.00801     0.725       0              0.08 
## 4 Summer 0.00793     0.755       0              0.079