1) Temperature

Temperature is the response variable in this exercise.

# adding Year-Month variable as date
climate1 <- climate %>%
  mutate(year_month = ymd(paste(climate$Year, climate$Month, truncated = 1))) 

# plotting Tempearture-Year 
ggplot(climate1, aes(year_month, Temp)) + 
  geom_line() + 
  geom_smooth(se=FALSE, linetype = "dotted") + 
  labs(title = "Temperature (1983-2010)",
       x = "Year", 
       y = "Temperature") +
  theme_minimal()

We can see from the plot that there has been a steady rise in temperature over the years but since 2005 the curve has plateaued off. If its a permanent shift in the trend or seasonality, we don’t know and discussing it is beyond the scope of this post. However, lets see if there’s any seasonality to temperature across the year.

# adding right and left hand sided labels
climate1 <- climate1 %>%
  mutate(label_rt = if_else(Month == 12 & Year%%2 == 0, 
                            as.character(Year), NA_character_)) %>%
  mutate(label_lt = if_else(Month == 1 & Year%%2 != 0, 
                            as.character(Year), NA_character_))

# creating label zones for left and right sided
x_lt <- c(NA,1)
x_rt <- c(12,NA)

# Temperature Month-wise plot for each year in data
ggplot(climate1, aes(as.factor(Month), Temp)) + 
  geom_point(aes(color = as.factor(Year))) + 
  geom_line(aes(group = as.factor(Year), 
                color = as.factor(Year)), 
            alpha = 0.7) + 
  labs(title = 'Temperature by month') +
  xlab("Months") +
  ylab("Temperature") + 
  geom_text_repel(aes(label = label_rt, color = as.factor(Year)), 
                  nudge_x = 1, xlim = x_rt) +
  geom_text_repel(aes(label = label_lt, color = as.factor(Year)), 
                  nudge_x = 1, xlim = x_lt) +
  scale_x_discrete(expand=c(0.1,0),
                   breaks=c("1", "2", "3", "4", "5", "6", 
                            "7", "8", "9", "10", "11", "12"),
                   labels=c("1", "2", "3", "4", "5", "6", 
                            "7", "8", "9", "10", "11", "12"),
                   limits=c("-1","1", "2", "3", "4", "5", "6", 
                            "7", "8", "9", "10", "11", "12", "13"), drop=FALSE) +
  theme(legend.position = "none")

We can see that generally the temperature has been steadily rising across years (same as last plot) but I find this plot is little bit cluttered, so lets try another approach and plot a ‘Temperature-density’ distribution.

ggplot(climate1, aes(x = Temp, y = as.factor(Year))) + 
  geom_density_ridges_gradient(aes(fill = ..x..), 
                               scale = 3, size = 0.3, alpha = 0.) +
  geom_vline(xintercept = 0.5, alpha = 0.5, color = "red", linetype = "dotted") +
   scale_fill_gradientn(colours = c("#87CEFA", "#FFFFE0", "#FF0000"),
                       name = "Temp") +
  labs(title = 'Temperature density',
       subtitle = "Temperature is difference in degrees Celsius between the \n average global temperature in that period and a reference value") + 
  theme(legend.position = c(0.9,0.2)) +
  xlab("Temperature") + 
  ylab("Year") + 
  theme_classic()

The plot reveals that for last few decades we have a permanent shift towards higher temperature. Although there is a dip from 2005 through 2008 as we did see in other plots, we notice that that extreme temperatures (>0.5°C) occurred in all years since 2000.

2) MEI

MEI i.e. “Multivariate El Nino Southern Oscillation index by Year” is a prime predictor for global climate disruptions. MEI is a time series type of data and its a combination of multiple variables which contains both atmospheric and oceanic variables. Real time monitoring of MEI enables goverments to tackle regional issue affected by climate and plan for food and water supply, health and safety etc.

climate1 %>%
  arrange(year_month) %>%
  ggplot() + 
  geom_line(aes(year_month, MEI, color = MEI), size = 1) +
  scale_color_gradient2(low = 'blue', high = 'red') + 
  labs(title = "Multivariate El Nino Southern Oscillation index by Year", 
       caption = "Source: ESRL/NOAA Physical Sciences Division") + 
  xlab("Year") +
  ylab("MEI") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme_minimal()

There is definitely a seasonality to MEI but not a continuous positive or negative trend over the years.

3) Carbon dioxide

Carbon dioxide is a Greenhouse house which traps the solar enegery and helps warm up the planet. Since Industrial revolution humans have significantly contributed to natrually occurring \(CO_2\) levels which might have significant impact on increasing Global temperature. The \(CO_2\) levels have been constantly increasing in atmosphere since this data is collected. In model building it’d be interesting to know how does temperature correlate with \(CO_2\) levels.

climate1 %>%
  arrange(year_month) %>%
  ggplot(aes(year_month, CO2)) + 
  geom_line() +
  geom_smooth(se = FALSE, linetype = "dotted", color = "red") +
  labs(title = "Carbon-dioxide by Year", 
       caption = "Source: ESRL/NOAA Global Monitoring Division") + 
  xlab("Year") +
  ylab("CO2 (ppmv)") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme_classic()

4) Methane

Methane is also a Greenhouse gas and traps more energy than Carbon dioxide. We can see that its levels are on constantly rise too, however we notice that the curve has somewhat flattened around 2000.

climate1 %>%
  arrange(year_month) %>%
  ggplot(aes(year_month, CH4)) + 
  geom_line() +
  geom_smooth(se=FALSE, linetype = "dotted", color = "red") +
  labs(title = "Methane by Year", 
       caption = "Source: ESRL/NOAA Global Monitoring Division") + 
  xlab("Year") +
  ylab("CH4 (ppmv)") + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme_classic()

5) Nitrous Oxide

Nitrous Oxide, another Green House gas, has lot more Global Warming Potential (GWP) than \(CO_2\) or \(N_2O\). As we can see all 3 biggest man-made contributor to Greenhouse gases are on constant rise.

climate1 %>%
  arrange(year_month) %>%
  ggplot() + 
  geom_line(aes(year_month, N2O)) +
  labs(title = "Nitrous Oxide by Year", 
       caption = "Source: ESRL/NOAA Global Monitoring Division") + 
  xlab("Year") +
  ylab("N2O (ppmv)") + 
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))

6) CFC-11 and CFC-12

We can see from the plots below that after rising steadily through 1995-2000 there is a decline in levels of CFC-11 and CFC-12. Read more about Kyoto protocol and Montreal protocol to understand the drivers behind this decline.

cfc11 <- climate1 %>%
  arrange(year_month) %>%
  ggplot() + 
  geom_line(aes(year_month, `CFC-11`)) +
  labs(title = "CFC-11",
       caption = "Source: ESRL/NOAA Global Monitoring Division") +
  xlab("Year") +
  ylab("CFC-11 (ppbv)") +
  theme_classic()

cfc12 <- climate1 %>%
  arrange(year_month) %>%
  ggplot() + 
  geom_line(aes(year_month, `CFC-12`)) +
  labs(title = "CFC-12",
       caption = "Source: ESRL/NOAA Global Monitoring Division") + 
  xlab("Year") +
  ylab("CFC-12 (ppbv)") + 
  theme_classic()

gridExtra::grid.arrange(cfc11, cfc12, nrow=1)

6) TSI

TSI pattern exhibits a 10 year cycle or at least that’s what it looks like from the plot. I was interested to know more about it so I googled “total solar irradiance cyclic?” and it took me to Solar Cycle wiki page which explains that “The solar cycle or solar magnetic activity cycle is the nearly periodic 11-year change in the Sun’s activity (including changes in the levels of solar radiation and ejection of solar material) and appearance (changes in the number and size of sunspots, flares, and other manifestations).”

climate1 %>%
  arrange(year_month) %>%
  ggplot(aes(year_month, TSI)) + 
  geom_line() + 
  labs(title = "Total Solar Irradiance (TSI)",
       caption = "Source: SOLARIS-HEPPA project") + 
  xlab("Year") +
  ylab("TSI (W/m^2)") +
  theme_classic()

7) Aerosols

Well at first glance it looks like something is off. The peak around 1992 looks like an anomaly or some data issue. However, if you go to Godard Institure for Space Studies at NASA website, you’ll notice that there is no problem with this data. It is accurate and a bit of research will tell you that in 1991 Mount Pinatubo erupted causing Aerosols to form a global layer of sulfuric acid haze which caused global temperatures to drop by 0.5°C.

climate1 %>%
  arrange(year_month) %>%
  ggplot() + 
  geom_line(aes(year_month, Aerosols)) + 
  labs(title = "Mean Stratospheric Aerosol Optical Depth at 550 nm",
       caption = "Source: Godard Institute for Space Studies at NASA") + 
  xlab("Year") +
  ylab("AOD")

Model Building

Before we build a Simple Linear model, lets understand what a “Baseline model” is and build one. A Baseline model predicts a unique value for a quantitative response variable or the most common value for qualitative/classification response variables. For example, if you’re predicting rainfall for next month, a baseline model will use average of all monthly rainfall data as a prediction, however, if you’re trying to predict who wins the next soccer world cup, a baseline model will look at team with highest wins and it’ll simply predict that team to win. In this particular example, baseline model will always predict the average temperature, which is:

mean(climate_train$Temp)

## [1] 0.2477993

Plotting Baseline model against \(CO_2\) will look like:

climate_train_baseline <- climate_train %>%
  mutate(baseline_temp = mean(Temp)) 

ggplot(climate_train_baseline, aes(CO2, baseline_temp)) + 
  geom_point(alpha = 0.2) + 
  stat_smooth(geom = "line", se = FALSE, color = "red", alpha = 0.8) +
  labs(title = "Baseline Temperature-CO2 model",
       x = "CO2",
       y = "Baseline Temperature") + 
  theme(plot.title = element_text(hjust = 0.5))

Basically, as stated above, this model will always predict same value for Temperature regardless of \(CO_2\) value. The error of this model, defined as difference between actual and predicted values, can be calculated as following:

climate_train_baseline <- climate_train_baseline %>%
  mutate(error = Temp - baseline_temp)

ggplot(climate_train_baseline, aes(CO2, error)) + 
  geom_point(aes(color = ifelse(error > 0, "red", "blue"))) +
  scale_color_identity() + 
  labs(title = "Error plot for Baseline model",
       x = "CO2",
       y = "Error")

We can see that we have, both, the positive and negative type of errors and they are distributed uniformly on both sides of zero. Above is a great plot to analyze errors visually but to compare the accuracy and strength of models we use few other statistics such as ‘Sum of Squared Errors’, ‘Root Mean Squared Errors’, ‘R-square’ etc.

The ‘Sum of Squared Errors’ (SSE) for baseline model is calculated as:

\[SSE = \sum_{i=1}^N(X_i - \bar{X})^{2} \]

(SSE_baseline <- sum(climate_train_baseline$error^2))

## [1] 9.285256

Note that SSE of Baseline model is also known as SST i.e Sum of Squared Total Errors.

SSE is a good indicator of model accuracy but it’s dependent on number of observations (N). So if we double the observation the model SSE will double as well which is NOT good for assessing accuracy of model, hence we use ‘Root Mean Squared Error’ (RMSE) which is basically square root of SSE divided by N.

\[RMSE = \sqrt[]{\frac{SSE}{N}}\]

(RMSE_baseline <- sqrt(SSE_baseline/nrow(climate_train_baseline)))

## [1] 0.1808164

Now we understand our baseline model so we’re ready to build our Simple Linear model.

\[Temp = \beta_0 + \beta_1*CO_2 + \epsilon \]

model1 <- lm(Temp ~ CO2, data = climate_train)

We can look at summary of model by using summary() function

summary(model1)

## 
## Call:
## lm(formula = Temp ~ CO2, data = climate_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31912 -0.07819 -0.00129  0.05805  0.43420 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.2646515  0.2096815  -20.34   <2e-16 ***
## CO2          0.0124855  0.0005799   21.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1116 on 282 degrees of freedom
## Multiple R-squared:  0.6218, Adjusted R-squared:  0.6204 
## F-statistic: 463.6 on 1 and 282 DF,  p-value: < 2.2e-16

We observe that coefficient of \(CO_2\) is 0.012 which means that for every increase of 10 ppmv of \(CO_2\), the temperature difference will increase by 0.12°C. The p value establishes that coefficient is significant.

Note: \(R^2\) ( = 1 - \(\frac{SSE}{SST}\)) depends on SSE and SST (SST is SSE of baseline). \(R^2\) value of model built using training dataset is anywhere between 0 and 1. A value of 0 means no improvement over Baseline and 1 means perfect Predictive model. In other words, \(R^2\) tells us how better our model is than Baseline model. The \(R^2\) for Simple model is 0.6218 which is significant improvement over Baseline model.

The summary() function also outputs ‘Adjusted \(R^2\)’, which adjusts the ‘Multiple \(R^2\)’ to account for the number of independent variables used relative to the number of data points. ‘Multiple \(R^2\)’ will always increase if you add independent variable whereas ‘Adjusted \(R^2\)’ will decrease if you add independent variable which decreases the quality of the model.

Lets calculate SSE and RMSE for Simple Linear model:

climate_train <- climate_train %>%
  add_predictions(model1) %>%
  mutate(error_simple = Temp - pred)

(SSE_simple_model = sum(climate_train$error_simple^2))

## [1] 3.511885

(RMSE_simple_model = sqrt(SSE_simple_model)/nrow(climate_train))

## [1] 0.0065986

We can see that Simple Model has significantly less SSE and RMSE compared to Baseline model.

Interpret model

We’re working with a single variable model and we can already see its much better than baseline model. Now lets plot the estimate versus actual data points for a Temperature-Carbon-dioxide plot.

aug_model1 <- augment(model1, climate_train)

ggplot(aug_model1, aes(CO2, Temp)) +
  geom_point(color = "red") + 
  geom_smooth(method = "lm", color = "blue", se = FALSE) +
  geom_segment(aes(x = CO2, xend = CO2,
                   y = Temp, yend = aug_model1$.fitted), alpha = 0.5) + 
  ggtitle("Regression Error") +
  theme_classic()

The summary() function is great to see relevant stats of a model but accessing the results from summary() is difficult as it’s not in a data frame or tibble format. Luckily, we can just do that using tidy() function from broom package. It’s a very handy tool for writing reports or futher analysis.

tidy(model1)

## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  -4.26    0.210        -20.3 3.13e-57
## 2 CO2           0.0125  0.000580      21.5 1.74e-61

Assessing accuracy

There are few ways we can assess accuracy of the model built against the raw data.

1. RSE
2. \(R^2\)
3. F-statistics

We can see that the \(R^2\) is 0.6217783 and adjusted \(R^2\) is 0.6204371 but we don’t know if its better or not since we haven’t build any other model. But we know its an improvement of ~62% over baseline model. We can look up all three values using glance() function from broom() package.

glance(model1)

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.622         0.620 0.112      464. 1.74e-61     2   221. -436. -425.
## # … with 2 more variables: deviance <dbl>, df.residual <int>

As we can see R-squared value is 0.622 (1 is perfect prediction), Residual Standard Error (RSE) is 0.112 (minimize as much as possible), and F-statistics is 463.594 (larger the better).

Making Prediction

We have our model ready, now we can use this model on to predict temperature changes on unseen data and assess its accuracy.

climate_test <- climate_test %>%
  add_predictions(model1)

climate_test %>%
  select(Temp, pred)

## # A tibble: 24 x 2
##     Temp  pred
##    <dbl> <dbl>
##  1 0.601 0.516
##  2 0.498 0.527
##  3 0.435 0.537
##  4 0.466 0.560
##  5 0.372 0.562
##  6 0.382 0.555
##  7 0.394 0.536
##  8 0.358 0.505
##  9 0.402 0.491
## 10 0.362 0.494
## # … with 14 more rows

Calculate Mean Squared Error (MSE)

# MSE (Mean Square Error)
climate_test %>%
  summarise(MSE = mean((Temp - pred)^2))

## # A tibble: 1 x 1
##      MSE
##    <dbl>
## 1 0.0438

We’re going to compare model1 against other models in next section.

Model Building

We can build a model where Temperature depends on all other variables except Year and Month.

Lets take a look at the summary() of this model.

summary(model2)

## 
## Call:
## lm(formula = Temp ~ MEI + CO2 + CH4 + N2O + `CFC-11` + `CFC-12` + 
##     TSI + Aerosols, data = climate_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25888 -0.05913 -0.00082  0.05649  0.32433 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.246e+02  1.989e+01  -6.265 1.43e-09 ***
## MEI          6.421e-02  6.470e-03   9.923  < 2e-16 ***
## CO2          6.457e-03  2.285e-03   2.826  0.00505 ** 
## CH4          1.240e-04  5.158e-04   0.240  0.81015    
## N2O         -1.653e-02  8.565e-03  -1.930  0.05467 .  
## `CFC-11`    -6.631e-03  1.626e-03  -4.078 5.96e-05 ***
## `CFC-12`     3.808e-03  1.014e-03   3.757  0.00021 ***
## TSI          9.314e-02  1.475e-02   6.313 1.10e-09 ***
## Aerosols    -1.538e+00  2.133e-01  -7.210 5.41e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09171 on 275 degrees of freedom
## Multiple R-squared:  0.7509, Adjusted R-squared:  0.7436 
## F-statistic: 103.6 on 8 and 275 DF,  p-value: < 2.2e-16

This is not a bad model. Our \(R^2\) increased from 0.622 to 0.751 and there are few significant predictors. It appears to be a better model than the simple linear model, however, there is a problem associated with this model, Collinearity.

Collinearity or Multicollinearity exists when there are two or more highly correlated variables in the predictor set. Because of this high collinearity among variables, it’s difficult to explain the variation in the dataset and which variable caused it. In short, we need to deal with this issue and the easiest way is to drop one or few highly collinear variables. To calculate collinearity we can use inbuilt cor() function.

round(cor(climate_train), digits = 2)

##               Year Month   MEI   CO2   CH4   N2O CFC-11 CFC-12   TSI Aerosols
## Year          1.00 -0.03 -0.04  0.98  0.92  0.99   0.57   0.90  0.17    -0.35
## Month        -0.03  1.00  0.00 -0.11  0.02  0.01  -0.01   0.00 -0.03     0.01
## MEI          -0.04  0.00  1.00 -0.04 -0.03 -0.05   0.07   0.01 -0.15     0.34
## CO2           0.98 -0.11 -0.04  1.00  0.88  0.98   0.51   0.85  0.18    -0.36
## CH4           0.92  0.02 -0.03  0.88  1.00  0.90   0.78   0.96  0.25    -0.27
## N2O           0.99  0.01 -0.05  0.98  0.90  1.00   0.52   0.87  0.20    -0.34
## CFC-11        0.57 -0.01  0.07  0.51  0.78  0.52   1.00   0.87  0.27    -0.04
## CFC-12        0.90  0.00  0.01  0.85  0.96  0.87   0.87   1.00  0.26    -0.23
## TSI           0.17 -0.03 -0.15  0.18  0.25  0.20   0.27   0.26  1.00     0.05
## Aerosols     -0.35  0.01  0.34 -0.36 -0.27 -0.34  -0.04  -0.23  0.05     1.00
## Temp          0.79 -0.10  0.17  0.79  0.70  0.78   0.41   0.69  0.24    -0.38
## pred          0.98 -0.11 -0.04  1.00  0.88  0.98   0.51   0.85  0.18    -0.36
## error_simple  0.02 -0.03  0.33  0.00  0.02  0.01   0.00   0.02  0.17    -0.17
##               Temp  pred error_simple
## Year          0.79  0.98         0.02
## Month        -0.10 -0.11        -0.03
## MEI           0.17 -0.04         0.33
## CO2           0.79  1.00         0.00
## CH4           0.70  0.88         0.02
## N2O           0.78  0.98         0.01
## CFC-11        0.41  0.51         0.00
## CFC-12        0.69  0.85         0.02
## TSI           0.24  0.18         0.17
## Aerosols     -0.38 -0.36        -0.17
## Temp          1.00  0.79         0.61
## pred          0.79  1.00         0.00
## error_simple  0.61  0.00         1.00

We can clearly see that there are few highly collinear variables in the dataset. However, there is a better way to visualize collinearity matrix using corrplot package.

library(corrplot)

x <- cor(climate_train)
corrplot.mixed(x, lower = "number", upper = "square", number.cex = 0.7,
         order = "hclust")

We can see that variable \(N_2O\) and \(CH_4\) have high correlation with other variables, so we can drop them from the dataset and build another model.

model3 <- lm(Temp ~ MEI + CO2 + `CFC-11` + `CFC-12` + TSI + Aerosols, 
                 data = climate_train)
summary(model3)

## 
## Call:
## lm(formula = Temp ~ MEI + CO2 + `CFC-11` + `CFC-12` + TSI + Aerosols, 
##     data = climate_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.26799 -0.06002 -0.00323  0.05651  0.33074 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.223e+02  1.991e+01  -6.139 2.86e-09 ***
## MEI          6.421e-02  6.465e-03   9.932  < 2e-16 ***
## CO2          4.061e-03  1.928e-03   2.106 0.036112 *  
## `CFC-11`    -4.314e-03  1.109e-03  -3.891 0.000125 ***
## `CFC-12`     2.430e-03  6.431e-04   3.778 0.000194 ***
## TSI          8.852e-02  1.461e-02   6.060 4.44e-09 ***
## Aerosols    -1.567e+00  2.132e-01  -7.347 2.28e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09201 on 277 degrees of freedom
## Multiple R-squared:  0.7475, Adjusted R-squared:  0.742 
## F-statistic: 136.6 on 6 and 277 DF,  p-value: < 2.2e-16

There is another method in R to work through this. R provides a step function which returns a model with optimized AIC. If we use model2 which contains all variables, step will figure out which variables to drop to maximize the accuracy of the model.

model4 <- step(model2)

## Start:  AIC=-1348.16
## Temp ~ MEI + CO2 + CH4 + N2O + `CFC-11` + `CFC-12` + TSI + Aerosols
## 
##            Df Sum of Sq    RSS     AIC
## - CH4       1   0.00049 2.3135 -1350.1
## <none>                  2.3130 -1348.2
## - N2O       1   0.03132 2.3443 -1346.3
## - CO2       1   0.06719 2.3802 -1342.0
## - `CFC-12`  1   0.11874 2.4318 -1335.9
## - `CFC-11`  1   0.13986 2.4529 -1333.5
## - TSI       1   0.33516 2.6482 -1311.7
## - Aerosols  1   0.43727 2.7503 -1301.0
## - MEI       1   0.82823 3.1412 -1263.2
## 
## Step:  AIC=-1350.1
## Temp ~ MEI + CO2 + N2O + `CFC-11` + `CFC-12` + TSI + Aerosols
## 
##            Df Sum of Sq    RSS     AIC
## <none>                  2.3135 -1350.1
## - N2O       1   0.03133 2.3448 -1348.3
## - CO2       1   0.06672 2.3802 -1344.0
## - `CFC-12`  1   0.13023 2.4437 -1336.5
## - `CFC-11`  1   0.13938 2.4529 -1335.5
## - TSI       1   0.33500 2.6485 -1313.7
## - Aerosols  1   0.43987 2.7534 -1302.7
## - MEI       1   0.83118 3.1447 -1264.9

Assessing Models Numerically

We can numerically compare quality of all four models.

list(model1 = broom::glance(model1), 
     model2 = broom::glance(model2),
     model3 = broom::glance(model3),
     model4 = broom::glance(model4))

## $model1
## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.622         0.620 0.112      464. 1.74e-61     2   221. -436. -425.
## # … with 2 more variables: deviance <dbl>, df.residual <int>
## 
## $model2
## # A tibble: 1 x 11
##   r.squared adj.r.squared  sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl>  <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.751         0.744 0.0917      104. 1.94e-78     9   280. -540. -504.
## # … with 2 more variables: deviance <dbl>, df.residual <int>
## 
## $model3
## # A tibble: 1 x 11
##   r.squared adj.r.squared  sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl>  <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.747         0.742 0.0920      137. 9.15e-80     7   278. -540. -511.
## # … with 2 more variables: deviance <dbl>, df.residual <int>
## 
## $model4
## # A tibble: 1 x 11
##   r.squared adj.r.squared  sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl>  <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.751         0.745 0.0916      119. 1.77e-79     8   280. -542. -509.
## # … with 2 more variables: deviance <dbl>, df.residual <int>

It is interesting to note that the step function does not address the collinearity of the variables, except that adding highly correlated variables will not improve the \(R^2\) significantly. The consequence of this is that the step function will not necessarily produce a very interpretable model - just a model that has balanced quality and simplicity for a particular weighting of quality and simplicity (AIC). Looking at the data above its clear that model2, model3 and model4 are roughly same in terms of quality.

Personally I’d use model3 for its interpretability and also because it addresses collinearity issue head on.

Assessing Models Visually

So far we compared the models numerically which is great but we should also assess model strength visually. We will choose model3 for this analysis. In each model a visual study of residuals can provide valuable insights.

aug_model3 <- augment(model3, climate_train)

ggplot(aug_model3, aes(.resid)) +
  geom_histogram(aes(y = ..density..), color = "black", fill = "white") +
  stat_function(fun = dnorm, color = "red",
                args = list(mean = mean(aug_model3$.resid, na.rm = TRUE),
                            sd = sd(aug_model3$.resid, na.rm = TRUE))) +
  xlab("Residuals from model3")

ggplot(aug_model3, aes(.fitted, .resid)) +
  geom_ref_line(h = 0) +
  geom_point(aes(color = as.factor(Year)), alpha = 0.6) + 
  geom_smooth(se = FALSE) + 
  ggtitle("Residuals vs Fitted") +
  xlab("Fitted values (Predictions)") + 
  ylab("Residuals (Actual - Predicted values)") +
  labs(color = "Year")

ggplot(aug_model3, aes(.fitted, .std.resid)) +
  geom_ref_line(h = 0) +
  geom_point(aes(color = as.factor(Year)), alpha = 0.6) +
  geom_smooth(se = FALSE) + 
  ggtitle("Standardized Residual vs Fitted") +
  xlab("Fitted values (Predictions)") + 
  ylab("Standardized Residuals") +
  labs(color = "Year")

ggplot(aug_model3, aes(.fitted, sqrt(abs(.std.resid)))) +
  geom_ref_line(h = 0) +
  geom_point(aes(color = as.factor(Year)), alpha = 0.6) +
  geom_smooth(se = FALSE) + 
  ggtitle("Scale-Location")  +
  xlab("Fitted values (Predictions)") + 
  ylab("sqrt(Standardized Residuals)") +
  labs(color = "Year")

ggplot(aug_model3, aes(sample = .resid)) +
  stat_qq(alpha = 0.3) +
  stat_qq_line(color = "blue") + 
  ggtitle("Normal Q-Q Plot") + 
  xlab("Theoretical Quantiles") + 
  ylab("Sample Quantiles")

plot(model3, which = 4, id.n = 5) 

# id.n identifies "n" top outlier observations

plot(model3, which = 5, id.n = 5)

Making Prediction

Having analyzed the model on training data both, numerically and visually, we’re now ready to make prediction on the testing dataset. We will use model3 for predictions.

(climate_test <- climate_test %>%
  add_predictions(model3)) %>%
  select(Temp, pred, everything())

## # A tibble: 24 x 12
##     Temp  pred  Year Month    MEI   CO2   CH4   N2O `CFC-11` `CFC-12`   TSI
##    <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl> <dbl>
##  1 0.601 0.488  2007     1  0.974  383. 1800.  321.     248.     539. 1366.
##  2 0.498 0.462  2007     2  0.51   384. 1803.  321.     248.     539. 1366.
##  3 0.435 0.442  2007     3  0.074  385. 1803.  321.     248.     539. 1366.
##  4 0.466 0.440  2007     4 -0.049  386. 1802.  321.     248.     539. 1366.
##  5 0.372 0.454  2007     5  0.183  387. 1796.  320.     247.     538. 1366.
##  6 0.382 0.423  2007     6 -0.358  386. 1782.  320.     247.     537. 1366.
##  7 0.394 0.421  2007     7 -0.290  384. 1772.  320.     246.     537. 1366.
##  8 0.358 0.403  2007     8 -0.44   382  1779.  320.     246.     537. 1366.
##  9 0.402 0.349  2007     9 -1.16   381. 1794.  321.     246.     537. 1366.
## 10 0.362 0.354  2007    10 -1.14   381. 1802.  321.     246.     537. 1366.
## # … with 14 more rows, and 1 more variable: Aerosols <dbl>

Comparing MSE of model3 versus model1:

climate_test %>%
  gather_predictions(model1, model3) %>%
  group_by(model) %>%
  summarise(MSE = mean((Temp-pred)^2))

## # A tibble: 2 x 2
##   model     MSE
##   <chr>   <dbl>
## 1 model1 0.0438
## 2 model3 0.0104

This clearly shows that when we use model1 and model3 for out of sample prediction, model3 is still a much better model with very low MSE.