Top 10 happiest countries :

Model 1 (lm.top10):

Linear regression assumes a linear relationship between variables. It is an approach for predicting the quantitative response Y based on single or multiple predictor variables X. Our objective is to examine the impact of GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity on Happiness_score for the top 10 happiest countries. Having included five independent variables into our model, we used the multiple linear regression model to address the pertaining question.

# split dataset
set.seed(123)
split <- sample(1:nrow(lm_top10), round(nrow(lm_top10) * 0.8)) # 80-20 split
train_top10 <- lm_top10[split, ]
test_top10 <- lm_top10[-split, ]

Estimation

The estimated model (lm.top10) was performed using the “lm” function in R using the Ordinary Least Square criterion to identify the “best fitting” line that minimizes the sum of the residuals (RSS) in estimating regression coefficients.

lm.top10 <- lm(Happiness_score ~., data = train_top10)

Interpretation

summary(lm.top10)

## 
## Call:
## lm(formula = Happiness_score ~ ., data = train_top10)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.16418 -0.06959 -0.00267  0.06118  0.18424 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         6.35035    0.76747   8.274 1.17e-09 ***
## GDP_per_capita      0.56045    0.25507   2.197  0.03491 *  
## Social_support     -0.05999    0.11674  -0.514  0.61066    
## Life_expectancy    -0.07528    0.31952  -0.236  0.81516    
## Freedom_to_choices  1.07110    0.55405   1.933  0.06157 .  
## Generosity         -0.63600    0.23109  -2.752  0.00943 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09638 on 34 degrees of freedom
## Multiple R-squared:  0.3748, Adjusted R-squared:  0.2829 
## F-statistic: 4.077 on 5 and 34 DF,  p-value: 0.005247

Estimated Happiness_score = 6.3504 + 0.5605 GDP_per_capita – 0.06 Social_support – 0.0753 Life_expectancy + 1.0711 Freedom_to_choices – 0.6360 Generosity

ß0 = 6.3504 suggests that if GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity are equal to 0, we would expect the Happiness_score to increase by 6.3504.

ß1 = 0.5605 suggests that if GDP_per_capita increases by 1, we would expect the Happiness_score to increase by 0.5605, holding all other variables constant.

ß2 = – 0.06 suggests that if Social_support increases by 1, we would expect the Happiness_score to decrease by 0.06, holding all other variables constant.

ß3 = – 0.0753 suggests that if Life_expectancy increases by 1, we would expect the Happiness_score to decrease by 0.0753, holding all other variables constant.

ß4 = 1.0711 suggests that if Freedom_to_choices increases by 1, we would expect the Happiness_score to increase by 1.0711, holding all other variables constant.

ß5 = – 0.6360 suggests that if Generosity increases by 1, we would expect the Happiness_score to decrease by 0.6360, holding all other variables constant.

Diagnosis

In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (lm.top10) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

We use residual plots to check the assumptions of an OLS linear regression model to determine whether the residuals are consistent with random error. Residual plot 1 displayed the residual values on the y-axis and fitted values for estimated model (lm.top10) on the x-axis.

Having included five independent variables into our model, we used the multiple linear regression model to address the pertaining question.

Residual plot and Q-Q plot for Model(lm.top10):

autoplot(lm.top10, which = 1:2, nrow = 2, ncol = 1 )

From the Residuals vs Fitted plot, it is observed that the residuals spread around a horizontal line without any distinct pattern. This is a good indication that we do not have non-linear relationships. As for the Normal Q-Q plot, it is observed that the residuals follow a straight line well and this shows that the residuals are normally distributed.

ii) Multicollinearity

The presence of multicollinearity can have a negative impact and can severely limit the conclusions of the research study. To identify if multicollinearity exists in our model (lm.top10), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated. The higher the value of VIF, the higher the correlation between that variable with the other variables.

vif(lm.top10)

##     GDP_per_capita     Social_support    Life_expectancy Freedom_to_choices 
##           1.582743           1.897888           2.589537           2.124132 
##         Generosity 
##           1.547655

From the results, it is observed that all five independent variables (GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity) are not highly correlated because all of them met the criterion cut-off point which is below 10. However, both scatter plot matrix and correlation matrix show that GDP_per_capita and Life_expectancy has a moderately strong positive correlation.

Prediction

confint(lm.top10, level = 0.95)

##                          2.5 %     97.5 %
## (Intercept)         4.79067125  7.9100319
## GDP_per_capita      0.04208937  1.0788179
## Social_support     -0.29722682  0.1772503
## Life_expectancy    -0.72461638  0.5740592
## Freedom_to_choices -0.05486567  2.1970703
## Generosity         -1.10563317 -0.1663627

The 95% confidence interval for GDP_per_capita is 0.0421 to 1.0788. We are 95% confidence that if GDP_per_capita increases by 1, the Happiness_score will increase between 0.0421 and 1.0788. This interval does not include zero and thus slope ß1 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Social_support, Life_expectancy and Freedom_to_choices are – 0.2972 to 0.1773, – 0.7246 to 0.5741 and – 0.0549 to 2.1971 respectively. These three intervals include zero and thus the slope ß2, ß3 and ß4 are not significantly different from zero at the 5% level of significance and the effects will not be significant at the 0.05 level as well.

The 95% confidence interval for Generosity is – 1.1056 to – 0.1664. We are 95% confident that if Generosity increases by 1, the Happiness_score will decrease between 0.1664 and 1.1056. This interval does not include zero and thus slope ß5 is significantly different from zero at the 5% level of significance.

Evaluation

cv_lm.top10 <- train(
  form = Happiness_score~., 
  data = train_top10, 
  method = "lm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_lm.top10

## Linear Regression 
## 
## 40 samples
##  5 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 36, 37, 36, 35, 36, 36, ... 
## Resampling results:
## 
##   RMSE        Rsquared   MAE       
##   0.09414616  0.6331909  0.08258742
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

By comparing RMSE and MAE, we can determine whether the forecast contains large but infrequent errors. The larger the difference between RMSE and MAE, the more inconsistent the error size. By using 10-fold cross validation, it was observed that model (lm.top10) has a RMSE of 0.0941 and a MAE of 0.0826. This suggests that there is some variation in the magnitude of the errors but very large errors are unlikely to occur. MAE of 0.0826 indicates that the average difference between the predicted and the actual Happiness_score was 0.0826. Furthermore, it has a R-squared of 0.9579. This suggests that 95.79% of the variation in Happiness_score is explained by the variation in GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity. It also has an AIC value of – 66.1384 and a BIC value of – 54.3162.

Improving the Model

The best subset selection was performed in order to identify a subset of all predictor(s) that was believed to be related to our response variable (Happiness_score).

regfit.2.1 <- regsubsets(Happiness_score~ ., data = train_top10) 

par(mfrow=c(2,2))
plot(regfit.2.1, scale = "r2")
plot(regfit.2.1, scale = "adjr2")
plot(regfit.2.1, scale = "Cp")
plot(regfit.2.1, scale = "bic")

For R-squared, it was observed that all the five independent variables (GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity) are the best subset of predictors. As for adjusted R-squared, Cp and Bayesian Information Criterion (BIC), all three of them have indicated that GDP_per_capita, Freedom_to_choices and Generosity were the best subset of predictors.

# create prediction function (global)
predict.regsubsets <- function(object, newdata, id, ...) {
  form <- as.formula(object$call[[2]])
  mat <- model.matrix(form, newdata)
  coefi <- coef(object, id = id)
  xvars <- names(coefi)
  mat[ , xvars] %*% coefi
}

Besides that, 10-fold cross validation was used to identify the best number of predictors that we should include in our model.

# k-fold CV method (2nd method for subset selection)
k <- 10  
set.seed(1)
folds <- sample(1:k, nrow(train_top10), replace = TRUE) 
cv_errors <- matrix(NA, k, 5, dimnames = list(NULL, paste(1:5)))

# Perform cross validation using a for loop where the folds that are equal to j are in the testing set
# While the remaining are in the training set
# We then perform our predictions for each model size using the predict function we have created earlier
# To compute the test erros on the appropriate subset and allocate them into their appropriate slots in the matrix cv_errors

for(j in 1:k) {
  best_fit <- regsubsets(Happiness_score ~ ., data = train_top10[folds != j, ], nvmax = 10) 
  for(i in 1:5) {
    pred <- predict(best_fit, train_top10[folds == j, ], id = i) 
    cv_errors[j,i] <- mean((train_top10$Happiness_score[folds == j] - pred)^2) 
  }
}

mean_cv_errors <- apply(cv_errors, 2, mean)
plot(mean_cv_errors, type = "b")

min_point <- which.min(mean_cv_errors)
points(min_point, mean_cv_errors[min_point], col = "red", cex = 2, pch = 20)

# k-fold CV method is preferred as it takes into account of any outliers

Based on the results above, it was observed that a total of 3 predictors was the best as it has the lowest mean cross validation error.

Based on the information that we have obtained from the above methods, we have concluded that all three independent variables GDP_per_capita, Freedom_to_choices and Generosity would be selected as our predictors for our next model where Social_support and Life_expectancy will be excluded from our original model.

Model 2 (lm.top10.2):

Our objective is to examine the impact of GDP_per_capita, Freedom_to_choices and Generosity on Happiness_score for the top 10 happiest countries. Having included three independent variables from the subset selection method and 10-fold cross validation method into our model, we used the multiple linear regression model to address the pertaining question.

Estimation

The estimated model (lm.top10.2) was performed using the “lm” function in R using the Ordinary Least Square criterion to identify the “best fitting” line that minimizes the sum of the residuals (RSS) in estimating regression coefficients.

lm.top10.2 <- lm(Happiness_score ~ GDP_per_capita + Freedom_to_choices + Generosity, data = train_top10)

Interpretation

summary(lm.top10.2)

## 
## Call:
## lm(formula = Happiness_score ~ GDP_per_capita + Freedom_to_choices + 
##     Generosity, data = train_top10)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.168591 -0.075685  0.006639  0.058902  0.189988 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          6.1506     0.4368  14.081 3.26e-16 ***
## GDP_per_capita       0.6016     0.2158   2.788  0.00842 ** 
## Freedom_to_choices   1.0213     0.3999   2.554  0.01503 *  
## Generosity          -0.5790     0.1962  -2.951  0.00555 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09455 on 36 degrees of freedom
## Multiple R-squared:  0.363,  Adjusted R-squared:  0.3099 
## F-statistic: 6.838 on 3 and 36 DF,  p-value: 0.0009173

Estimated Happiness_score = 6.1506 + 0.6016 GDP_per_capita + 1.0213 Freedom_to_choices – 0.5790 Generosity

ß0 = 6.1506 suggests that if GDP_per_capita, Freedom_to_choices and Generosity are equal to 0, we would expect the Happiness_score to increase by 6.1506.

ß1 = 0.6016 suggests that if GDP_per_capita increases by 1, we would expect the Happiness_score to increase by 0.6016, holding all other variables constant.

ß2 = 1.0213 suggests that if Freedom_to_choices increases by 1, we would expect the Happiness_score to increase by 1.0213, holding all other variables constant.

ß3 = – 0.5790 suggests that if Generosity increases by 1, we would expect the Happiness_score to decrease by 0.5790, holding all other variables constant.

Diagnosis

In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (lm.top10.2) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(lm.top10.2):

autoplot(lm.top10.2, which = 1:2, nrow = 2, ncol = 1 )

From the Residuals vs Fitted plot, it is observed that the residuals spread around a horizontal line without any distinct pattern. This is a good indication that we do not have non-linear relationships. As for the Normal Q-Q plot, it is observed that the residuals follow a straight line well and this shows that the residuals are normally distributed.

ii) Multicollinearity

To identify if multicollinearity exists in our model (lm.top10.2), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated.

vif(lm.top10.2)

##     GDP_per_capita Freedom_to_choices         Generosity 
##           1.177075           1.149950           1.159661

From the results, it is observed that all three independent variables (GDP_per_capita, Freedom_to_choices and Generosity) are not highly correlated because all of them met the criterion cut-off point which is below 10.

Prediction

confint(lm.top10.2, level = 0.95)

##                         2.5 %     97.5 %
## (Intercept)         5.2647310  7.0364827
## GDP_per_capita      0.1639943  1.0392421
## Freedom_to_choices  0.2102981  1.8323842
## Generosity         -0.9769778 -0.1810234

The 95% confidence interval for GDP_per_capita is 0.1640 to 1.0392. We are 95% confident that if GDP_per_capita increases by 1, the Happiness_score will increase between 0.1640 and 1.0392. This interval does not include zero and thus slope ß1 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Freedom_to_choices is 0.2103 to 1.8324. We are 95% confident that if Freedom_to_choices increases by 1, the Happiness_score will increase between 0.2103 and 1.8324. This interval does not include zero and thus slope ß2 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Generosity is – 0.9770 to – 0.1810. We are 95% confident that if Generosity increases by 1, the Happiness_score will decrease between 0.1810 and 0.9770. This interval does not include zero and thus slope ß3 is significantly different from zero at the 5% level of significance.

Evaluation

cv_lm.top10.2 <- train(
  form = Happiness_score ~ GDP_per_capita + Freedom_to_choices + Generosity, 
  data = train_top10, 
  method = "lm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_lm.top10.2

## Linear Regression 
## 
## 40 samples
##  3 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 36, 36, 36, 37, 36, 36, ... 
## Resampling results:
## 
##   RMSE        Rsquared   MAE       
##   0.09625309  0.5317547  0.08041897
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

By using 10-fold cross validation, it was observed that model (lm.top10.2) has a RMSE of 0.0963 and a MAE of 0.0804. This suggests that there is some variation in the magnitude of the errors but very large errors are unlikely to occur. MAE of 0.0804 indicates that the average difference between the predicted and the actual Happiness_score was 0.0804. Furthermore, it has a R-squared of 0.8581. This suggests that 85.81% of the variation in Happiness_score is explained by the variation in GDP_per_capita, Freedom_to_choices and Generosity. It also has an AIC value of – 69.3889 and a BIC value of – 60.9445.

Improving the Model

regfit.2.2 <- regsubsets(Happiness_score ~ GDP_per_capita + Freedom_to_choices + Generosity, data = train_top10) 

par(mfrow=c(2,2))
plot(regfit.2.2, scale = "r2")
plot(regfit.2.2, scale = "adjr2")
plot(regfit.2.2, scale = "Cp")
plot(regfit.2.2, scale = "bic")

The best subset selection was performed and it was observed that R-squared, adjusted R-squared, Cp and Bayesian Information Criterion (BIC) have indicated that GDP_per_capita, Freedom_to_choices and Generosity were the best subset of predictors.

# k-fold CV method (2nd method for subset selection)
k <- 10  
set.seed(1)
folds <- sample(1:k, nrow(train_top10), replace = TRUE) 
cv_errors <- matrix(NA, k, 3, dimnames = list(NULL, paste(1:3)))

# Perform cross validation using a for loop where the folds that are equal to j are in the testing set
# While the remaining are in the training set
# We then perform our predictions for each model size using the predict function we have created earlier
# To compute the test erros on the appropriate subset and allocate them into their appropriate slots in the matrix cv_errors

for(j in 1:k) {
  best_fit <- regsubsets(Happiness_score ~ GDP_per_capita + Freedom_to_choices + Generosity, data = train_top10[folds != j, ], nvmax = 10) 
  for(i in 1:3) {
    pred <- predict(best_fit, train_top10[folds == j, ], id = i) 
    cv_errors[j,i] <- mean((train_top10$Happiness_score[folds == j] - pred)^2) 
  }
}

mean_cv_errors <- apply(cv_errors, 2, mean)
plot(mean_cv_errors, type = "b")

min_point <- which.min(mean_cv_errors)
points(min_point, mean_cv_errors[min_point], col = "red", cex = 2, pch = 20)

# k-fold CV method is preferred as it takes into account of any outliers

By using 10-fold cross validation, it was observed that a total of 3 predictors was the best as it has the lowest mean cross validation error.

Based on the information that we have obtained from the above methods, we have concluded that all three independent variables GDP_per_capita, Freedom_to_choices and Generosity would be selected as predictors for our model.

Model 3 (lm.top10.3):

Our objective is to examine the impact of GDP_per_capita and Generosity on Happiness_score for the top 10 happiest countries. We have made some assumptions that removing predictors with no relationship with the response variable (Happiness_score) will lead to a more effective model. As predictors that do not have a significant impact on the response variable tend to cause a deterioration in the test error rate. Having included two independent variables with the smallest p-value into our model as our second method, we used the multiple linear regression model to address the pertaining question.

Estimation

The estimated model (lm.top10.3) was performed using the “lm” function in R using the Ordinary Least Square criterion to identify the “best fitting” line that minimizes the sum of the residuals (RSS) in estimating regression coefficients.

lm.top10.3 <- lm(Happiness_score ~ GDP_per_capita + Generosity, data = train_top10)

Interpretation

summary(lm.top10.3)

## 
## Call:
## lm(formula = Happiness_score ~ GDP_per_capita + Generosity, data = train_top10)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.19199 -0.08023 -0.00747  0.09584  0.15389 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      6.9996     0.3038  23.042  < 2e-16 ***
## GDP_per_capita   0.4713     0.2248   2.097  0.04290 *  
## Generosity      -0.6814     0.2059  -3.309  0.00209 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1014 on 37 degrees of freedom
## Multiple R-squared:  0.2476, Adjusted R-squared:  0.2069 
## F-statistic: 6.087 on 2 and 37 DF,  p-value: 0.005185

Estimated Happiness_score = 6.9996 + 0.4713 GDP_per_capita – 0.6814 Generosity

ß0 = 6.9996 suggests that if GDP_per_capita and Generosity are equal to 0, we would expect the Happiness_score to increase by 6.9996.

ß1 = 0.4713 suggests that if GDP_per_capita increases by 1, we would expect the Happiness_score to increase 0.4713, holding all other variables constant.

ß2 = – 0.6814 suggests that if Generosity increases by 1, we would expect theHappiness_score to decrease by 0.6814, holding all other variables constant.

Diagnosis

In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (lm.top10.3) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(lm.top10.3):

autoplot(lm.top10.3, which = 1:2, nrow = 2, ncol = 1)

From the Residuals vs Fitted plot, it is observed that the residuals spread around a horizontal line without any distinct pattern. This is a good indication that we do not have non-linear relationships. As for the Normal Q-Q plot, it is observed that the residuals deviate slightly at the beginning and towards the end. This shows that the residuals are not normally distributed.

ii) Multicollinearity

To identify if multicollinearity exists in our model (lm.top10.3), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated.

vif(lm.top10.3)

## GDP_per_capita     Generosity 
##       1.111253       1.111253

From the results, it is observed that all two independent variables (GDP_per_capita and Generosity) are not highly correlated because all of them met the criterion cut-off point which is below 10.

Prediction

confint(lm.top10.3, level = 0.95)

##                      2.5 %     97.5 %
## (Intercept)     6.38408762  7.6150981
## GDP_per_capita  0.01588258  0.9267137
## Generosity     -1.09865002 -0.2641399

The 95% confidence interval for GDP_per_capita is 0.0159 to 0.9267. We are 95% confident that if GDP_per_capita increases by 1, the Happiness_score will increase between 0.0159 and 0.9267. This interval does not include zero and thus slope ß1 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Generosity is – 1.0987 to – 0.2641. We are 95% confident that if Generosity increases by 1, the Happiness_score will decrease between 0.2641 and 1.0987. This interval does not include zero and thus slope ß2 is significantly different from zero at the 5% level of significance.

Evaluation

cv_lm.top10.3 <- train(
  form = Happiness_score ~ GDP_per_capita + Generosity, 
  data = train_top10, 
  method = "lm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_lm.top10.3

## Linear Regression 
## 
## 40 samples
##  2 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 36, 36, 36, 37, 36, 36, ... 
## Resampling results:
## 
##   RMSE        Rsquared   MAE       
##   0.09996555  0.5070088  0.08767941
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

By using 10-fold cross validation, it was observed that model (lm.top10.3) has a RMSE of 0.09997 and a MAE of 0.0877. This suggests that there is some variation in the magnitude of the errors but very large errors are unlikely to occur. MAE of 0.0877 indicates that the average difference between the predicted and the actual Happiness_score was 0.0877. Furthermore, it has a R-squared of 0.9175. This suggests that 91.75% of the variation in Happiness_score is explained by the variation in GDP_per_capita and Generosity. It also has an AIC value of – 64.7281 and a BIC value of – 57.9726.

Improving the Model

regfit.2.3 <- regsubsets(Happiness_score ~ GDP_per_capita + Generosity, data = train_top10) 

par(mfrow=c(2,2))
plot(regfit.2.3, scale = "r2")
plot(regfit.2.3, scale = "adjr2")
plot(regfit.2.3, scale = "Cp")
plot(regfit.2.3, scale = "bic")

The best subset selection was performed and it was observed that R-squared, adjusted R-squared, Cp and Bayesian Information Criterion (BIC) have indicated that GDP_per_capita and Generosity were the best subset of predictors.

# k-fold CV method (2nd method for subset selection)
k <- 10  
set.seed(1)
folds <- sample(1:k, nrow(train_top10), replace = TRUE) 
cv_errors <- matrix(NA, k, 2, dimnames = list(NULL, paste(1:2)))

# Perform cross validation using a for loop where the folds that are equal to j are in the testing set
# While the remaining are in the training set
# We then perform our predictions for each model size using the predict function we have created earlier
# To compute the test erros on the appropriate subset and allocate them into their appropriate slots in the matrix cv_errors

for(j in 1:k) {
  best_fit <- regsubsets(Happiness_score ~ GDP_per_capita + Generosity, data = train_top10[folds != j, ], nvmax = 10) 
  for(i in 1:2) {
    pred <- predict(best_fit, train_top10[folds == j, ], id = i) 
    cv_errors[j,i] <- mean((train_top10$Happiness_score[folds == j] - pred)^2) 
  }
}

mean_cv_errors <- apply(cv_errors, 2, mean)
plot(mean_cv_errors, type = "b")

min_point <- which.min(mean_cv_errors)
points(min_point, mean_cv_errors[min_point], col = "red", cex = 2, pch = 20)

# k-fold CV method is preferred as it takes into account of any outliers

By using 10-fold cross validation, it was observed that a total of 2 predictors was the best as it has the lowest mean cross validation error.

Based on the information that we have obtained from the above methods, we have concluded that all two independent variables GDP_per_capita and Generosity would be considered as predictors for our model.

Model 4 (lm.top10.4):

Our objective is to examine the impact of Generosity on Happiness_score for the top 10 happiest countries. Having included only one independent variable which is adjusted from the p-value method into our model, we used the simple linear regression model to address the pertaining question.

Estimation

The estimated model (lm.top10.4) was performed using the “lm” function in R using the Ordinary Least Square criterion to identify the “best fitting” line that minimizes the sum of the residuals (RSS) in estimating regression coefficients.

lm.top10.4 <- lm(Happiness_score ~ Generosity, data = train_top10)

Interpretation

summary(lm.top10.4)

## 
## Call:
## lm(formula = Happiness_score ~ Generosity, data = train_top10)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.21021 -0.08858  0.01068  0.09240  0.14424 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.6191     0.0738 103.241   <2e-16 ***
## Generosity   -0.5448     0.2039  -2.672    0.011 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1058 on 38 degrees of freedom
## Multiple R-squared:  0.1581, Adjusted R-squared:  0.136 
## F-statistic: 7.139 on 1 and 38 DF,  p-value: 0.01105

Estimated Happiness_score = 7.6191 – 0.5448 Generosity

ß0 = 7.6191 suggests that if Generosity is equal to 0, we would expect the Happiness_score to increase by 7.6191.

ß1 = – 0.5448 suggests that if Generosity increases by 1, we would expect the Happiness_score to decrease by 0.5448, holding all other variables constant.

Diagnosis

In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (lm.top10.4) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model (lm.top10.4):

autoplot(lm.top10.4, which = 1:2, nrow = 2, ncol = 1 )

From the Residuals vs Fitted plot, it is observed that the residuals spread around a horizontal line without any distinct pattern. This is a good indication that we do not have non-linear relationships. As for the Normal Q-Q plot, it is observed that the residuals deviate slightly at the beginning and towards the end. This shows that the residuals are not normally distributed.

Prediction

confint(lm.top10.4, level = 0.95)

##                  2.5 %     97.5 %
## (Intercept)  7.4696712  7.7684668
## Generosity  -0.9575303 -0.1320038

The 95% confidence interval for Generosity is – 0.9575 to – 0.1320. We are 95% confident that if Generosity increases by 1, the Happiness_score will decrease between 0.1320 and 0.9575. This interval does not include zero and thus slope ß1 is significantly different from zero at the 5% level of significance.

Evaluation

cv_lm.top10.4 <- train(
  form = Happiness_score ~ Generosity, 
  data = train_top10, 
  method = "lm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_lm.top10.4

## Linear Regression 
## 
## 40 samples
##  1 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 36, 36, 36, 37, 36, 36, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE       
##   0.1049372  0.4903464  0.09337823
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

By using 10-fold cross validation, it was observed that model (lm.top10.4) has a RMSE of 0.1049 and a MAE of 0.0934. This suggests that there is some variation in the magnitude of the errors but very large errors are unlikely to occur. MAE of 0.0934 indicates that the average difference between the predicted and the actual Happiness_score was 0.0934. Furthermore, it has a R-squared of 0.9225. This suggests that 92.25% of the variation in Happiness_score is explained by the variation in Generosity. It also has an AIC value of – 62.2366 and a BIC value of – 57.1700.

Top 10 least happiest countries :

# split for bottom 10
set.seed(123)
split <- sample(1:nrow(lm_bottom10), round(nrow(lm_bottom10) * 0.8)) # 80-20 split
train_bottom10 <- lm_bottom10[split, ]
test_bottom10 <- lm_bottom10[-split, ]

Model 1 (lm.bottom10):

Our objective is to examine the impact of GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity on Happiness_score for the top 10 least happiest countries. Having included all five independent variables into our model, we used the multiple linear regression model to address the pertaining question.

Estimation

The estimated model (lm.bottom10) was performed using the “lm” function in R using the Ordinary Least Square criterion to identify the “best fitting” line that minimizes the sum of the residuals (RSS) in estimating regression coefficients.

lm.bottom10 <- lm(Happiness_score ~., data = train_bottom10)

Interpretation

summary(lm.bottom10)

## 
## Call:
## lm(formula = Happiness_score ~ ., data = train_bottom10)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50935 -0.28452  0.00991  0.25540  0.91697 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          3.3204     0.2288  14.514 6.96e-16 ***
## GDP_per_capita       0.4872     0.4530   1.076    0.290    
## Social_support       0.3632     0.2269   1.601    0.119    
## Life_expectancy     -0.2551     0.5487  -0.465    0.645    
## Freedom_to_choices   0.1810     0.3810   0.475    0.638    
## Generosity          -0.7439     0.7753  -0.960    0.344    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3573 on 33 degrees of freedom
## Multiple R-squared:  0.1865, Adjusted R-squared:  0.06329 
## F-statistic: 1.514 on 5 and 33 DF,  p-value: 0.2124

Estimated Happiness_score = 3.3204 + 0.4872 GDP_per_capita + 0.3632 Social_support – 0.2551 Life_expectancy + 0.1810 Freedom_to_choices – 0.7439 Generosity

ß0 = 3.3204 suggests that if GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity are equal to 0, we would expect the Happiness_score to increase by 3.3204.

ß1 = 0.4872 suggests that if GDP_per_capita increases by 1, we would expect the Happiness_score to increase by 0.4872, holding all other variables constant.

ß2 = 0.3632 suggests that if the Social_support increases by 1, we would expect the Happiness_score to increase by 0.3632, holding all other variables constant.

ß3 = – 0.2551 suggests that if the Life_expectancy increases by 1, we would expect the Happiness_score to decrease by 0.2551, holding all other variables constant.

ß4 = 0.1810 suggests that if the Freedom_to_choices increases by 1, we would expect the Happiness_score to increase by 0.1810, holding all other variables constant.

ß5 = – 0.7439 suggests that if the Generosity increases by 1, we would expect the Happiness_score to decrease by 0.7439, holding all other variables constant.

Diagnosis

In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (lm.bottom10) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(lm.bottom10):

autoplot(lm.bottom10, which = 1:2, nrow = 2, ncol = 1 )

From the Residuals vs Fitted plot, it is observed that the residuals spread around a horizontal line without any distinct pattern. This is a good indication that we do not have non-linear relationships. As for the Normal Q-Q plot, it is observed that the residuals deviate slightly at the beginning and towards the end. This shows that the residuals are not normally distributed.

ii) Multicollinearity

To identify if multicollinearity exists in our model (lm.bottom10), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated.

vif(lm.bottom10)

##     GDP_per_capita     Social_support    Life_expectancy Freedom_to_choices 
##           2.459384           1.475577           2.113764           1.225632 
##         Generosity 
##           1.658770

From the results, it is observed that all five independent variables (GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity) are not highly correlated because all of them met the criterion cut-off point which is below 10. However, both scatter plot matrix and correlation matrix show that GDP_per_capita and Life_expectancy has a moderately strong positive correlation.

Prediction

confint(lm.bottom10, level = 0.95)

##                          2.5 %    97.5 %
## (Intercept)         2.85495497 3.7858197
## GDP_per_capita     -0.43440708 1.4087866
## Social_support     -0.09839511 0.8247848
## Life_expectancy    -1.37146300 0.8611739
## Freedom_to_choices -0.59414031 0.9560817
## Generosity         -2.32113431 0.8333956

The 95% confidence interval for GDP_per_capita is – 0.4344 to 1.4088. We are 95% confident that if GDP_per_capita increases by 1, the Happiness_score will either decrease by 0.4344 or increase by 1.4088.

The 95% confidence interval for Social_support is – 0.0984 to 0.8248. We are 95% confident that if Social_support increases by 1, the Happiness_score will either decrease by 0.0984 or increase by 0.8248.

The 95% confidence interval for Life_expectancy is – 1.3715 to 0.8612. We are 95% confident that if Life_expectancy increases by 1, the Happiness_score will either decrease by 1.3715 or increase by 0.8612.

The 95% confidence interval for Freedom_to_choices is – 0.5941 to 0.9561. We are 95% confident that if Freedom_to_choices increases by 1, the Happiness_score will either decrease by 0.5941 or increase by 0.9561.

The 95% confidence interval for Generosity is – 2.3211 to 0.8334. We are 95% confident that if Generosity increases by 1, the Happiness_score will either decrease by 2.3211 or increase by 0.8334.

All these five intervals include zero and thus the slope ß1, ß2, ß3, ß4 and ß5 are not significantly different from zero at the 5% level of significance and the effects will not be significant at the 0.05 level as well.

Evaluation

cv_lm.bottom10 <- train(
  form = Happiness_score~., 
  data = train_bottom10, 
  method = "lm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_lm.bottom10

## Linear Regression 
## 
## 39 samples
##  5 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 35, 35, 35, 35, 35, 35, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.3653142  0.2794597  0.3121033
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

By using 10-fold cross validation, it was observed that model (lm.bottom10) has a RMSE of 0.3653 and a MAE of 0.3121. This suggests that there is some variation in the magnitude of the errors but very large errors are unlikely to occur. MAE of 0.3121 indicates that the average difference between the predicted and the actual Happiness_score was 0.3121. Furthermore, it has a R-squared of 0.9876. This suggests that 98.76% of the variation in Happiness_score is explained by the variation in GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity. It also has an AIC value of 37.8817 and a BIC value of 49.5267.

Improving the Model

regfit.2.3 <- regsubsets(Happiness_score~ ., data = train_bottom10) 

par(mfrow=c(2,2))
plot(regfit.2.3, scale = "r2")
plot(regfit.2.3, scale = "adjr2")
plot(regfit.2.3, scale = "Cp")
plot(regfit.2.3, scale = "bic")

The best subset selection was performed. For R-squared, it was observed that all the five independent variables (GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity) are the best subset of predictors. As for adjusted R-squared, Cp and Bayesian Information Criterion (BIC), all three of them have indicated that Social_support was the best subset of predictors.

# k-fold CV method (2nd method for subset selection)
k <- 5 
set.seed(1)
folds <- sample(1:k, nrow(train_bottom10), replace = TRUE) 
cv_errors <- matrix(NA, k, 5, dimnames = list(NULL, paste(1:5)))

# Perform cross validation using a for loop where the folds that are equal to j are in the testing set
# While the remaining are in the training set
# We then perform our predictions for each model size using the predict function we have created earlier
# To compute the test erros on the appropriate subset and allocate them into their appropriate slots in the matrix cv_errors

for(j in 1:k) {
  best_fit <- regsubsets(Happiness_score ~ ., data = train_bottom10[folds != j, ], nvmax = 10) 
  for(i in 1:5) {
    pred <- predict(best_fit, train_bottom10[folds == j, ], id = i) 
    cv_errors[j,i] <- mean((train_bottom10$Happiness_score[folds == j] - pred)^2) 
  }
}

mean_cv_errors <- apply(cv_errors, 2, mean)
plot(mean_cv_errors, type = "b")

min_point <- which.min(mean_cv_errors)
points(min_point, mean_cv_errors[min_point], col = "red", cex = 2, pch = 20)

# k-fold CV method is preferred as it takes into account of any outliers

By using 10-fold cross validation, it was observed that to include only one predictor was the best as it has the lowest mean cross validation error.

Based on the information that we have obtained from the above methods, we have concluded that one independent variable Social_support would be selected as our predictor for our next model where GDP_per_capita, Life_expectancy, Freedom_to_choices and Generosity will be excluded from our original model.

Model 2 (lm.bottom10.2):

Our objective is to examine the impact of Social_support on Happiness_score for the top 10 least happiest countries. Having included only one independent variable into our model, we used the simple linear regression model to address the pertaining question.

Estimation

The estimated model (lm.bottom10.2) was performed using the “lm” function in R using the Ordinary Least Square criterion to identify the “best fitting” line that minimizes the sum of the residuals (RSS) in estimating regression coefficients.

lm.bottom10.2 <- lm(Happiness_score ~ Social_support, data = train_bottom10)

Interpretation

summary(lm.bottom10.2)

## 
## Call:
## lm(formula = Happiness_score ~ Social_support, data = train_bottom10)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.59052 -0.28765  0.07949  0.20367  0.94779 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.2084     0.1113  28.818   <2e-16 ***
## Social_support   0.4557     0.1807   2.523   0.0161 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3456 on 37 degrees of freedom
## Multiple R-squared:  0.1467, Adjusted R-squared:  0.1237 
## F-statistic: 6.363 on 1 and 37 DF,  p-value: 0.01608

Estimated Happiness_score = 3.2084 + 0.4557 Social_support

ß0 = 3.2084 suggests that if Social_support is equal to 0, we would expect the Happiness_score to increase by 3.2084.

ß1 = 0.4557 suggests that if the Social_support increases by 1, we would expect the Happiness_score to increase by 0.4557, holding all other variables constant.

Diagnosis

In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (lm.bottom10.2) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(lm.bottom10.2):

autoplot(lm.bottom10.2, which = 1:2, nrow = 2, ncol = 1 )

From the Residuals vs Fitted plot, it is observed that the residuals spread around a horizontal line without any distinct pattern. This is a good indication that we do not have non-linear relationships. As for the Normal Q-Q plot, it is observed that the residuals follow a straight line well and this shows that the residuals are normally distributed.

Prediction

confint(lm.bottom10.2, level = 0.95)

##                     2.5 %    97.5 %
## (Intercept)    2.98285193 3.4340155
## Social_support 0.08966412 0.8217365

The 95% confidence interval for Social_support is 0.0897 to 0.8217. We are 95% confident that if Social_support increases by 1, the Happiness_score will increase between 0.0897 and 0.8217. This interval does not include zero and thus slope ß1 is significantly different from zero at the 5% level of significance.

Evaluation

cv_lm.bottom10.2 <- train(
  form = Happiness_score ~ Social_support, 
  data = train_bottom10, 
  method = "lm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_lm.bottom10.2

## Linear Regression 
## 
## 39 samples
##  1 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 35, 35, 35, 35, 35, 35, ... 
## Resampling results:
## 
##   RMSE       Rsquared  MAE      
##   0.3360635  0.565561  0.2894052
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

By using 10-fold cross validation, it was observed that model (lm.bottom10.2) has a RMSE of 0.3361 and a MAE of 0.2894. This suggests that there is some variation in the magnitude of the errors but very large errors are unlikely to occur. MAE of 0.2894 indicates that the average difference between the predicted and the actual Happiness_score was 0.2894. Furthermore, it has a R-squared of 0.9959. This suggests that 99.59% of the variation in Happiness_score is explained by the variation in Social_support. It also has an AIC value of 31.7447 and a BIC value of 36.7354.

Best model selection :

Comparing goodness of fit between models :

Top 10 happiest countries

Based on the table above, it is observed that Model (lm.top10) has the highest R-squared value (95.79%) and the lowest RMSE value (0.0941) as compared to other models.

Besides that, the RMSE, MAE and AIC values in Model (lm.top10.3) and Model (lm.top10.4) are higher than in Model (lm.top10) and Model (lm.top10.2), thus both Model (lm.top10.3) and Model (lm.top10.4) do not fit the data well.

Although the MAE, AIC and BIC values in Model (lm.top10.2) are the smallest which are 0.0804, – 69.39 and – 60.94 respectively, the variables (GDP_per_capita, Freedom_to_choices and Generosity) in Model (lm.top10.2) are not highly correlated. Furthermore, Model (lm.top10.2) has the lowest R-squared value (85.81%) which is approximately 10% lower than the R-squared in Model (lm.top10) which indicates that Model (lm.top10.2) do not perform better than model (lm.top10). Hence, we can conclude that Model (lm.top10) is our best model from the 4 models that we have generated in predicting the Happiness_score for the top 10 happiest countries as it has the largest R-squared value and the smallest RMSE value. These metrics have indicated that the data could fit Model (lm.top10) well where it is best to use all 5 variables which are GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity.

Top 10 least happiest countries

By comparing between Model (lm.bottom10) and Model (lm.bottom10.2), we can see that Model (lm.bottom10.2) has a higher R-squared value of 99.59% than Model (lm.bottom10). This indicates that Model (lm.bottom10.2) is better in terms of explaining the variation in Happiness_score.

Besides that, the RMSE, MAE, AIC and BIC values in Model (lm.bottom10.2) are 0.3361, 0.2894, 37.74 and 36.74 respectively which are smaller than (lm.bottom10.2) as shown in the table above. Hence, we can clearly see that Model (lm.bottom10.2) is the best model in predicting the Happiness_score for the top 10 least happiest countries where it is best to only use Social_support as the predictor for this model.

Model 1 (lm.western):

Our objective is to examine the impact of GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity on Happiness_score for the countries in Western Europe region. Having included five independent variables into our model, we used the multiple linear regression model to address the pertaining question.

# split data
set.seed(123)
split <- sample(1:nrow(lm_western), round(nrow(lm_western) * 0.8)) # 80-20 split
train_western <- lm_western[split, ]
test_western <- lm_western[-split, ]

Estimation

The estimated model (lm.western) was performed using the “lm” function in R using the Ordinary Least Square criterion to identify the “best fitting” line that minimizes the sum of the residuals (RSS) in estimating regression coefficients.

lm.western <- lm(Happiness_score ~., data = train_western)

Interpretation

summary(lm.western)

## 
## Call:
## lm(formula = Happiness_score ~ ., data = train_western)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90620 -0.28174  0.05924  0.28728  0.69152 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.6037     0.8155  -0.740 0.461448    
## GDP_per_capita       2.3344     0.4126   5.658 2.58e-07 ***
## Social_support       0.4663     0.2457   1.898 0.061501 .  
## Life_expectancy      2.2455     0.6535   3.436 0.000959 ***
## Freedom_to_choices   2.1746     0.4199   5.179 1.77e-06 ***
## Generosity           1.3392     0.4240   3.159 0.002275 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3762 on 76 degrees of freedom
## Multiple R-squared:  0.7827, Adjusted R-squared:  0.7684 
## F-statistic: 54.75 on 5 and 76 DF,  p-value: < 2.2e-16

Estimated Happiness_score = – 0.6037 + 2.3344 GDP_per_capita + 0.4663 Social_support + 2.2455 Life_expectancy + 2.1746 Freedom_to_choices + 1.3392 Generosity

ß0 = -0.6037 suggests that if GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity are equal to 0, we would expect the Happiness_score to decrease by 0.6037.

ß1 = 2.3344 suggests that if GDP_per_capita increases by 1, we would expect the Happiness_score to increase by 2.3344, holding all other variables constant.

ß2 = 0.4663 suggests that if the Social_support increases by 1, we would expect the Happiness_score to increase by 0.4663, holding all other variables constant.

ß3 = 2.2455 suggests that if the Life_expectancy increases by 1, we would expect the Happiness_score to increase by 2.2455, holding all other variables constant.

ß4 = 2.1746 suggests that if the Freedom_to_choices increases by 1, we would expect the Happiness_score to increase by 2.1746, holding all other variables constant.

ß5 = 1.3392 suggests that if the Generosity increases by 1, we would expect the Happiness_score to increase by 1.3392, holding all other variables constant.

Diagnosis

In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (lm.western) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(lm.western):

autoplot(lm.western, which = 1:2, nrow = 2, ncol = 1 )

From the Residuals vs Fitted plot, it is observed that the residuals spread around a horizontal line without any distinct pattern. This is a good indication that we do not have non-linear relationships. As for the Normal Q-Q plot, it is observed that the residuals deviate slightly towards the end. This shows that the residuals are not normally distributed.

ii) Multicollinearity

To identify if multicollinearity exists in our model (lm.western), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated.

vif(lm.western)

##     GDP_per_capita     Social_support    Life_expectancy Freedom_to_choices 
##           1.393302           1.668647           1.337811           2.337401 
##         Generosity 
##           1.797192

From the results, it is observed that all five independent variables (GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity) are not highly correlated because all of them met the criterion cut-off point which is below 10. However, both scatter plot matrix and correlation matrix show that GDP_per_capita and Life_expectancy has a moderately strong positive correlation.

Prediction

confint(lm.western, level = 0.95)

##                          2.5 %    97.5 %
## (Intercept)        -2.22789749 1.0205799
## GDP_per_capita      1.51261570 3.1562197
## Social_support     -0.02302934 0.9557007
## Life_expectancy     0.94397226 3.5470490
## Freedom_to_choices  1.33837023 3.0108981
## Generosity          0.49477505 2.1836590

The 95% confidence interval for GDP_per_capita is 1.5126 to 3.1562. We are 95% confident that if GDP_per_capita increases by 1, the Happiness_score will increase between 1.5126 and 3.1562. This interval does not include zero and thus slope ß1 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Social_support is – 0.0230 to 0.9557 respectively. This interval includes zero and thus the slope ß2 is not significantly different from zero at the 5% level of significance and the effects will not be significant at the 0.05 level as well.

The 95% confidence interval for Life_expectancy is 0.9440 to 3.5470. We are 95% confident that if Life_expectancy increases by 1, the Happiness_score will increase between 0.9440 and 3.5470. This interval does not include zero and thus slope ß3 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Freedom_to_choices is 1.3384 to 3.0109. We are 95% confident that if Freedom_to_choices increases by 1, the Happiness_score will increase between 1.3384 and 3.0109. This interval does not include zero and thus slope ß4 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Generosity is 0.4948 to 2.1837. We are 95% confident that if Generosity increases by 1, the Happiness_score will increase between 0.4948 and 2.1837. This interval does not include zero and thus slope ß5 is significantly different from zero at the 5% level of significance.

Evaluation

cv_lm.western <- train(
  form = Happiness_score~., 
  data = train_western, 
  method = "lm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_lm.western

## Linear Regression 
## 
## 82 samples
##  5 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 74, 74, 74, 73, 74, 74, ... 
## Resampling results:
## 
##   RMSE       Rsquared  MAE      
##   0.3856888  0.791963  0.3176116
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

By using 10-fold cross validation, it was observed that the model (lm.western) has a RMSE of 0.3857 and a MAE of 0.3176. This suggests that there is some variation in the magnitude of the errors but very large errors are unlikely to occur. MAE of 0.3176 indicates that the average difference between the predicted and the actual Happiness_score was 0.3176. Furthermore, it has a R-squared of 0.7920. This suggests that 79.20% of the variation in Happiness_score is explained by the variation in GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity. It also has an AIC value of 80.1321 and a BIC value of 96.9791.

Improving the Model

regfit.3.1 <- regsubsets(Happiness_score~ ., data = train_western) 

par(mfrow=c(2,2))
plot(regfit.3.1, scale = "r2")
plot(regfit.3.1, scale = "adjr2")
plot(regfit.3.1, scale = "Cp")
plot(regfit.3.1, scale = "bic")

The best subset selection was performed. For R-squared, adjusted R-squared and Cp, all three of them have indicated that GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity were the best subset of predictors. As for Bayesian Information Criterion (BIC), it was observed that four independent variables (GDP_per_capita, Life_expectancy, Freedom_to_choices and Generosity) are the best subset of predictors.

# k-fold CV method (2nd method for subset selection)
k <- 10  
set.seed(1)
folds <- sample(1:k, nrow(train_western), replace = TRUE) 
cv_errors <- matrix(NA, k, 5, dimnames = list(NULL, paste(1:5)))

# Perform cross validation using a for loop where the folds that are equal to j are in the testing set
# While the remaining are in the training set
# We then perform our predictions for each model size using the predict function we have created earlier
# To compute the test erros on the appropriate subset and allocate them into their appropriate slots in the matrix cv_errors

for(j in 1:k) {
  best_fit <- regsubsets(Happiness_score ~ ., data = train_western[folds != j, ], nvmax = 10) 
  for(i in 1:5) {
    pred <- predict(best_fit, train_western[folds == j, ], id = i) 
    cv_errors[j,i] <- mean((train_western$Happiness_score[folds == j] - pred)^2) 
  }
}

mean_cv_errors <- apply(cv_errors, 2, mean)
plot(mean_cv_errors, type = "b")

min_point <- which.min(mean_cv_errors)
points(min_point, mean_cv_errors[min_point], col = "red", cex = 2, pch = 20)

# k-fold CV method is preferred as it takes into account of any outliers

By using 10-fold cross validation, it was observed that a total of 5 predictors were the best as it has the lowest mean cross validation error.

Based on the information that we have obtained from the above methods, we have concluded that four independent variables GDP_per_capita, Life_expectancy, Freedom_to_choices and Generosity would be selected as our predictors for our next model where Social_support will be excluded from our original model.

Model 2 (lm.western.2):

Our objective is to examine the impact of GDP_per_capita, Life_expectancy, Freedom_to_choices and Generosity on Happiness_score for the countries in Western Europe region. Having included four independent variables from the subset selection method and 10-fold cross validation method into our model, we used the multiple linear regression model to address the pertaining question.

Estimation

The estimated model (lm.western.2) was performed using the “lm” function in R using the Ordinary Least Square criterion to identify the “best fitting” line that minimizes the sum of the residuals (RSS) in estimating regression coefficients.

lm.western.2 <- lm(Happiness_score ~ GDP_per_capita + Life_expectancy + Freedom_to_choices + Generosity, data = train_western)

Interpretation

summary(lm.western.2)

## 
## Call:
## lm(formula = Happiness_score ~ GDP_per_capita + Life_expectancy + 
##     Freedom_to_choices + Generosity, data = train_western)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.86317 -0.26073  0.09358  0.26050  0.64164 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.8330     0.8200  -1.016  0.31288    
## GDP_per_capita       2.5115     0.4087   6.145 3.27e-08 ***
## Life_expectancy      2.7650     0.6034   4.583 1.74e-05 ***
## Freedom_to_choices   2.5049     0.3885   6.447 9.07e-09 ***
## Generosity           1.2452     0.4281   2.908  0.00474 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3825 on 77 degrees of freedom
## Multiple R-squared:  0.7724, Adjusted R-squared:  0.7606 
## F-statistic: 65.33 on 4 and 77 DF,  p-value: < 2.2e-16

Estimated Happiness_score = – 0.8330 + 2.5115 GDP_per_capita + 2.7650 Life_expectancy + 2.5049 Freedom_to_choices + 1.2452 Generosity

ß0 = -0.8330 suggests that if GDP_per_capita, Life_expectancy, Freedom_to_choices and Generosity are equal to 0, we would expect the Happiness_score to decrease by 0.8330.

ß1 = 2.5115 suggests that if GDP_per_capita increases by 1, we would expect the Happiness_score to increase 2.5115, holding all other variables constant.

ß2 = 2.7650 suggests that if the Life_expectancy increases by 1, we would expect the Happiness_score to increase by 2.7650, holding all other variables constant.

ß3 = 2.5049 suggests that if the Freedom_to_choices increases by 1, we would expect the Happiness_score to increase by 2.5049, holding all other variables constant.

ß4 = 1.2452 suggests that if the Generosity increases by 1, we would expect the Happiness_score to increase by 1.2452, holding all other variables constant.

Diagnosis

In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (lm.western.2) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(lm.western.2):

autoplot(lm.western.2, which = 1:2, nrow = 2, ncol = 1 )

From the Residuals vs Fitted plot, it is observed that the residuals spread around a horizontal line without any distinct pattern. This is a good indication that we do not have non-linear relationships. As for the Normal Q-Q plot, it is observed that the residuals deviate slightly towards the end. This shows that the residuals are not normally distributed.

ii) Multicollinearity

To identify if multicollinearity exists in our model (lm.western.2), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated.

vif(lm.western.2)

##     GDP_per_capita    Life_expectancy Freedom_to_choices         Generosity 
##           1.322098           1.103131           1.935925           1.772675

From the results, it is observed that all four independent variables (GDP_per_capita, Life_expectancy, Freedom_to_choices and Generosity) are not highly correlated because all of them met the criterion cut-off point which is below 10. However, both scatter plot matrix and correlation matrix show that GDP_per_capita and Life_expectancy has a moderately strong positive correlation.

Prediction

confint(lm.western.2, level = 0.95)

##                         2.5 %    97.5 %
## (Intercept)        -2.4659007 0.7998601
## GDP_per_capita      1.6976830 3.3252255
## Life_expectancy     1.5635529 3.9664160
## Freedom_to_choices  1.7312519 3.2785608
## Generosity          0.3926921 2.0977654

The 95% confidence interval for GDP_per_capita is 1.6977 to 3.3252. We are 95% confident that if GDP_per_capita increases by 1, the Happiness_score will increase between 1.6977 and 3.3252. This interval does not include zero and thus slope ß1 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Life_expectancy is 1.5636 to 3.9664. We are 95% confident that if Life_expectancy increases by 1, the Happiness_scoree will increase between 1.5636 and 3.9664. This interval does not include zero and thus slope ß2 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Freedom_to_choices is 1.7313 to 3.2786. We are 95% confident that if Freedom_to_choices increases by 1, the Happiness_score will increase between 1.7313 and 3.2786. This interval does not include zero and thus slope ß3 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Generosity is 0.3927 to 2.0978. We are 95% confident that if Generosity increases by 1, the Happiness_score will increase between 0.3927 and 2.0978. This interval does not include zero and thus slope ß4 is significantly different from zero at the 5% level of significance.

Evaluation

cv_lm.western.2 <- train(
  form = Happiness_score ~ GDP_per_capita + Life_expectancy + Freedom_to_choices + Generosity, 
  data = train_western, 
  method = "lm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_lm.western.2

## Linear Regression 
## 
## 82 samples
##  4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 74, 74, 74, 73, 74, 74, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.3813104  0.7993535  0.3191809
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

By using 10-fold cross validation, it was observed that model (lm.western.2) has a RMSE of 0.3813 and a MAE of 0.3192. This suggests that there is some variation in the magnitude of the errors but very large errors are unlikely to occur. MAE of 0.3192 indicates that the average difference between the predicted and the actual Happiness_score was 0.3192. Furthermore, it has a R-squared of 0.7994. This suggests that 79.94% of the variation in Happiness_score is explained by the variation in GDP_per_capita, Life_expectancy, Freedom_to_choices and Generosity. It also has an AIC value of 81.9294 and a BIC value of 96.3697.

Improving the Model

regfit.2.2 <- regsubsets(Happiness_score ~ GDP_per_capita + Life_expectancy + Freedom_to_choices + Generosity, data = train_western)

par(mfrow=c(2,2))
plot(regfit.2.2, scale = "r2")
plot(regfit.2.2, scale = "adjr2")
plot(regfit.2.2, scale = "Cp")
plot(regfit.2.2, scale = "bic")

The best subset selection was performed. As for R-squared, adjusted R-squared, Cp and Bayesian Information Criterion (BIC), all four of them have indicated that GDP_per_capita, Life_expectancy, Freedom_to_choices and Generosity were the best subset of predictors.

# k-fold CV method (2nd method for subset selection)
k <- 10  
set.seed(1)
folds <- sample(1:k, nrow(train_western), replace = TRUE) 
cv_errors <- matrix(NA, k, 4, dimnames = list(NULL, paste(1:4)))

# Perform cross validation using a for loop where the folds that are equal to j are in the testing set
# While the remaining are in the training set
# We then perform our predictions for each model size using the predict function we have created earlier
# To compute the test erros on the appropriate subset and allocate them into their appropriate slots in the matrix cv_errors

for(j in 1:k) {
  best_fit <- regsubsets(Happiness_score ~ GDP_per_capita + Life_expectancy + Freedom_to_choices + Generosity, data = train_western[folds != j, ], nvmax = 10) 
  for(i in 1:4) {
    pred <- predict(best_fit, train_western[folds == j, ], id = i) 
    cv_errors[j,i] <- mean((train_western$Happiness_score[folds == j] - pred)^2) 
  }
}

mean_cv_errors <- apply(cv_errors, 2, mean)
plot(mean_cv_errors, type = "b")

min_point <- which.min(mean_cv_errors)
points(min_point, mean_cv_errors[min_point], col = "red", cex = 2, pch = 20)

# k-fold CV method is preferred as it takes into account of any outliers

By using 10-fold cross validation, it was observed that a total of 4 predictors were the best as it has the lowest mean cross validation error.

Based on the information that we have obtained from the above methods, we have concluded that four independent variables GDP_per_capita, Life_expectancy, Freedom_to_choices and Generosity would be considered as predictors for our model.

Model 3 (lm.western3):

Our objective is to examine the impact of GDP_per_capita, Life_expectancy and Freedom_to_choices on Happiness_score for the countries in Western Europe region. Having included three independent variables with the smallest p-value into our model, we used the multiple linear regression model to address the pertaining question.

Estimation

The estimated model (lm.western.3) was performed using the “lm” function in R using the Ordinary Least Square criterion to identify the “best fitting” line that minimizes the sum of the residuals (RSS) in estimating regression coefficients.

lm.western.3 <- lm(Happiness_score ~ GDP_per_capita + Life_expectancy + Freedom_to_choices, data = train_western)

Interpretation

summary(lm.western.3)

## 
## Call:
## lm(formula = Happiness_score ~ GDP_per_capita + Life_expectancy + 
##     Freedom_to_choices, data = train_western)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0944 -0.2057  0.1007  0.2481  0.6583 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.5611     0.8527  -0.658 0.512472    
## GDP_per_capita       2.5475     0.4276   5.958 6.96e-08 ***
## Life_expectancy      2.4115     0.6186   3.898 0.000204 ***
## Freedom_to_choices   3.1748     0.3275   9.693 4.91e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4003 on 78 degrees of freedom
## Multiple R-squared:  0.7474, Adjusted R-squared:  0.7377 
## F-statistic: 76.93 on 3 and 78 DF,  p-value: < 2.2e-16

Estimated Happiness_score = – 0.5611 + 2.5475 GDP_per_capita + 2.4115 Life_expectancy + 3.1748 Freedom_to_choices

ß0 = -0.5611 suggests that if GDP_per_capita, Life_expectancy and Freedom_to_choices and Generosity are equal to 0, we would expect the Happiness_score to decrease by 0.5611.

ß1 = 2.5475 suggests that if GDP_per_capita increases by 1, we would expect the Happiness_score to increase 2.5475, holding all other variables constant.

ß2 = 2.4115 suggests that if the Life_expectancy increases by 1, we would expect the Happiness_score to increase by 2.4115, holding all other variables constant.

ß3 = 3.1748 suggests that if the Freedom_to_choices increases by 1, we would expect the Happiness_score to increase by 3.1748, holding all other variables constant.

Diagnosis

In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (lm.western.3) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(lm.western.3):

autoplot(lm.western.3, which = 1:2, nrow = 2, ncol = 1)

From the Residuals vs Fitted plot, it is observed that the residuals spread around a horizontal line without any distinct pattern. This is a good indication that we do not have non-linear relationships. As for the Normal Q-Q plot, it is observed that the residuals deviate slightly at the beginning and towards the end. This shows that the residuals are not normally distributed.

ii) Multicollinearity

To identify if multicollinearity exists in our model (lm.western.3), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated.

vif(lm.western.3)

##     GDP_per_capita    Life_expectancy Freedom_to_choices 
##           1.320881           1.058364           1.255651

From the results, it is observed that all three independent variables (GDP_per_capita, Life_expectancy and Freedom_to_choices) are not highly correlated because all of them met the criterion cut-off point which is below 10. However, both scatter plot matrix and correlation matrix show that GDP_per_capita and Life_expectancy has a moderately strong positive correlation.

Prediction

confint(lm.western.3, level = 0.95)

##                        2.5 %   97.5 %
## (Intercept)        -2.258793 1.136573
## GDP_per_capita      1.696295 3.398749
## Life_expectancy     1.179942 3.643008
## Freedom_to_choices  2.522712 3.826810

The 95% confidence interval for GDP_per_capita is 1.6963 to 3.3987. We are 95% confident that if GDP_per_capita increases by 1, the Happiness_score will increase between 1.6963 and 3.3987. This interval does not include zero and thus slope ß1 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Life_expectancy is 1.1799 to 3.6430. We are 95% confident that if Life_expectancy increases by 1, the Happiness_score will increase between 1.1799 and 3.6430. This interval does not include zero and thus slope ß2 is significantly different from zero at the 5% level of significance.

The 95% confidence interval for Freedom_to_choices is 2.5227 to 3.8268. We are 95% confident that if Freedom_to_choices increases by 1, the Happiness_score will increase between 2.5227 and 3.8268. This interval does not include zero and thus slope ß3 is significantly different from zero at the 5% level of significance.

Evaluation

cv_lm.western.3 <- train(
  form = Happiness_score ~ GDP_per_capita + Life_expectancy + Freedom_to_choices, 
  data = train_western, 
  method = "lm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_lm.western.3

## Linear Regression 
## 
## 82 samples
##  3 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 74, 74, 74, 73, 74, 74, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.3983567  0.7799532  0.3205176
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

By using 10-fold cross validation, it was observed that model (lm.western.3) has a RMSE of 0.3984 and a MAE of 0.3205. This suggests that there is some variation in the magnitude of the errors but very large errors are unlikely to occur. MAE of 0.3205 indicates that the average difference between the predicted and the actual Happiness_score was 0.3205. Furthermore, it has a R-squared of 0.78. This suggests that 78% of the variation in Happiness_score is explained by the variation in GDP_per_capita, Life_expectancy and Freedom_to_choices. It also has an AIC value of 88.4765 and a BIC value of 100.5101.

Improving the Model

regfit.2.3 <- regsubsets(Happiness_score ~ GDP_per_capita + Life_expectancy + Freedom_to_choices, data = train_western) 

par(mfrow=c(2,2))
plot(regfit.2.3, scale = "r2")
plot(regfit.2.3, scale = "adjr2")
plot(regfit.2.3, scale = "Cp")
plot(regfit.2.3, scale = "bic")

The best subset selection was performed. As for R-squared, adjusted R-squared, Cp and Bayesian Information Criterion (BIC), all four of them have indicated that GDP_per_capita, Life_expectancy and Freedom_to_choices were the best subset of predictors.

# k-fold CV method (2nd method for subset selection)
k <- 10  
set.seed(1)
folds <- sample(1:k, nrow(train_western), replace = TRUE) 
cv_errors <- matrix(NA, k, 3, dimnames = list(NULL, paste(1:3)))

# Perform cross validation using a for loop where the folds that are equal to j are in the testing set
# While the remaining are in the training set
# We then perform our predictions for each model size using the predict function we have created earlier
# To compute the test erros on the appropriate subset and allocate them into their appropriate slots in the matrix cv_errors

for(j in 1:k) {
  best_fit <- regsubsets(Happiness_score ~ GDP_per_capita + Life_expectancy + Freedom_to_choices, data = train_western[folds != j, ], nvmax = 10) 
  for(i in 1:3) {
    pred <- predict(best_fit, train_western[folds == j, ], id = i) 
    cv_errors[j,i] <- mean((train_western$Happiness_score[folds == j] - pred)^2) 
  }
}

mean_cv_errors <- apply(cv_errors, 2, mean)
plot(mean_cv_errors, type = "b")

min_point <- which.min(mean_cv_errors)
points(min_point, mean_cv_errors[min_point], col = "red", cex = 2, pch = 20)

# k-fold CV method is preferred as it takes into account of any outliers

By using 10-fold cross validation, it was observed that a total of 3 predictors were the best as it has the lowest mean cross validation error.

Based on the information that we have obtained from the above methods, we have concluded that three independent variables GDP_per_capita, Life_expectancy and Freedom_to_choices would be considered as predictors for our model.

Best model selection :

Comparing goodness of fit between models :

Referring to the table above, we can see that although Model (lm.western.3) has the highest R-squared value (93.45%), the RMSE, MAE, AIC and BIC values in Model (lm.western.3) are the highest as compared to the other models. Hence, Model (lm.western.3) does not fit the data well.

Besides that, by comparing between Model (lm.western) and Model (lm.western.2), it is observed that Model (lm.western) has a higher R-squared value of 92.32% than Model (lm.western.2). This indicates that Model (lm.western) is better in terms of explaining the variation in Happiness_score.

Furthermore, the MAE and AIC values in Model (lm.western) are 0.3176 and 80.13 respectively which are smaller than the MAE and AIC values in Model (lm.western.2). Although Model (lm.western.2) has a lower RMSE and BIC values which are 0.3813 and 96.37 respectively as compared to Model (lm.western), Model (lm.western) is still selected as the best model because there is not much difference in RMSE and BIC values between both models. Therefore, Model (lm.western) is the best in predicting the Happiness_score for the countries in Western Europe region where it is best to use all 5 variables which are GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity.

Model 1 (glm.world):

Logistic regression is used to investigate the relationship between the binary dependent variables with other independent variables. We classify the instances into two groups (0 and 1) by calculating the probability of them in one group. In this research question, we would examine the impact of GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity on Happiness_score, where 0 represents low happiness score (with the happiness score less than the mean) and 1 represents high happiness score.

# split data into training and testing set 
set.seed(123)
split <- sample(1:nrow(lm_world), round(nrow(lm_world) * 0.8)) # 80-20 split
train_world <- lm_world[split, ]
test_world <- lm_world[-split, ]

Estimation

The estimated model (glm.world) was performed using the “glm” function in R using the logistic regression to identify the association of the binary independent variable with other dependent variables and fit the line values to the sigmoid curves.

glm.world <- glm(Happy~., data = train_world, family = binomial)

Interpretation

summary(glm.world)

## 
## Call:
## glm(formula = Happy ~ ., family = binomial, data = train_world)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.38803  -0.26160   0.05166   0.37126   2.73839  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -13.0128     1.2115 -10.741  < 2e-16 ***
## GDP_per_capita       5.3536     0.7608   7.037 1.97e-12 ***
## Social_support       1.8573     0.5791   3.207  0.00134 ** 
## Life_expectancy      4.0072     0.9887   4.053 5.06e-05 ***
## Freedom_to_choices   7.6789     1.2339   6.223 4.87e-10 ***
## Generosity           0.1316     1.2400   0.106  0.91549    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 851.08  on 613  degrees of freedom
## Residual deviance: 347.46  on 608  degrees of freedom
## AIC: 359.46
## 
## Number of Fisher Scoring iterations: 7

Estimated happiness score = -13.0128 + 5.3536 GDP per capita + 1.8573 Social support + 4.0072 Life expectancy + 7.6789 Freedom to choices + 0.1316 Generosity

ß0 = -13.0128 suggests that the probability of happiness score not being 1 is 13.01%

ß1 = 5.3536 suggests that if the GDP per capita increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(4.7110^{-4}\)%, holding all other variables constant. \(\frac{e^{(-13.0128+5.3536)}}{1+e^{(-13.0128+5.3536)}}=4.7110^{-4}\)%

ß2 = 1.8573 suggests that if the social support increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(1.4310^{-5}\)%, holding all other variables constant. \(\frac{e^{(-13.0128+1.8573)}}{1+e^{(-13.0128+1.8573)}}=1.4310 ^{-5}\)%

ß3 = 4.0072 suggests that if the life expectancy increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(1.2310^{-4}\)%, holding all other variables constant. \(\frac{e^{(-13.0128+4.0072)}}{1+e^{(-13.0128+4.0072)}}=1.2310 ^{-4}\)%

ß4 = 7.6789 suggests that if the freedom to choices increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(4.8010^{-3}\)%, holding all other variables constant. \(\frac{e^{(-13.0128+7.6789)}}{1+e^{(-13.0128+7.6789)}}=4.8010^{-3}\)%

ß5 = 0.1316 suggests that if the generosity increases by 1, the probability of having a happiness score of 1 increases multiplicatively by 4.8010-3%, holding all other variables constant. \(\frac{e^{(-13.0128+0.1316)}}{1+e^{(-13.0128+0.1316)}}=2.5510^{-6}\)%

Diagnosis

Regression diagnostics are a powerful ally in the model and data validation stages of regression model building. It allows us to assess the overall adequacy of our models. In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (glm.model) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(glm.world):

autoplot(glm.world, which = 1:2, nrow = 2, ncol = 1 )

Based on the Residuals vs Fitted plot, the residuals are equally spread without any distinct patterns, hence we could conclude that there are no non-linear relationships in this model. However, we could observe that the residuals deviate slightly at the beginning and towards the end in the Normal Q-Q plot, this shows that the residuals are not normally distributed.

ii) Multicollinearity

To identify if multicollinearity exists in our model (glm.world), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated.

vif(glm.world) 

##     GDP_per_capita     Social_support    Life_expectancy Freedom_to_choices 
##           1.382507           1.059949           1.219299           1.427642 
##         Generosity 
##           1.166290

From the results, it is observed that all five independent variables (GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity) are not highly correlated as all variables have met the criterion cut-off point which is below 10. However, both scatter plot matrix and correlation matrix show that GDP_per_capita and Life_expectancy has a moderately strong positive correlation.

Prediction

confint(glm.world, level = 0.95)

## Waiting for profiling to be done...

##                          2.5 %     97.5 %
## (Intercept)        -15.5649930 -10.801557
## GDP_per_capita       3.9371722   6.927085
## Social_support       0.7423078   3.017803
## Life_expectancy      2.1284406   6.018697
## Freedom_to_choices   5.3405752  10.192480
## Generosity          -2.2980744   2.575804

The 95% confidence interval for GDP_per_capita and Social_support are 3.94 to 6.93 and 0.74 to 3.02 respectively. Since the interval doesn’t include 0, hence we could conclude that ß1 and ß2 are significant at the 5% level of significance. Similarly, Life_expectancy and Freedom_to_choices have an interval of 2.13 to 6.02 and 5.34 to 10.19 respectively at the 5% significance level, hence ß3 and ß4 are significant since 0 is not included.

However, the 95% confidence interval for Generosity is -2.30 to 2.58. The interval includes 0 and thus ß5 is not significant at the 5% level of significance.

Evaluation

# k fold cv
cv_glm.world <- train(
  form = Happy~., 
  data = train_world, 
  method = "glm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_glm.world

## Generalized Linear Model 
## 
## 614 samples
##   5 predictor
##   2 classes: '0', '1' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 553, 553, 553, 552, 552, 552, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.8633527  0.7263887

Through 10-fold cross validation, we could observe that this model has an accuracy of 0.8634, this indicates that the accuracy for the correctly predicted observations to the total observations is as high as 86.34%. Furthermore, the kappa value for this model is 0.7264. Kappa is measured to examine the reliability of the collected data in studying the variables of the model, whether the data could represent the variables correctly. In this model, we could conclude that the collected data are reliable and could represent the variables correctly.

Improving the Model

regfit.1.1 <- regsubsets(Happy ~ ., data = train_world)

par(mfrow=c(2,2))
plot(regfit.1.1, scale = "r2")
plot(regfit.1.1, scale = "adjr2")
plot(regfit.1.1, scale = "Cp")
plot(regfit.1.1, scale = "bic")

The best subset selection is performed to find out the best fitted variables for the model. R-squared has indicated that all 5 variables would be the best subset of predictors for this model. However, adjusted R-squared, Cp and Bayesian Information Criterion (BIC) have suggested that only GDP_per_capita, Social_support, Life_expectancy and Freedom_to_choices were the best subset of predictors.

Therefore, we have concluded that three independent variables GDP_per_capita, Social_support, Life_expectancy and Freedom_to_choices would be the best predictors for our next model.

Model 2 (glm.world.2):

Our objective would be to examine the impact of the 4 significant variables, which are GDP_per_capita, Social_support, Life_expectancy and Freedom_to_choices on Happiness_score, where 0 represents low happiness score (with the happiness score less than the mean) and 1 represents high happiness score. We have excluded Generosity from the previous model as it is insignificant at 5% significance level in the previous model and its p-value is 0.92.

Estimation

The estimated model (glm.world.2) was performed using the “glm” function in R using the logistic regression to identify the association of the binary independent variable with other dependent variables and fit the line values to the sigmoid curves.

glm.world.2 <- glm(Happy ~ GDP_per_capita + Social_support + Life_expectancy + Freedom_to_choices, data = train_world, family = binomial) 

Interpretation

summary(glm.world.2)

## 
## Call:
## glm(formula = Happy ~ GDP_per_capita + Social_support + Life_expectancy + 
##     Freedom_to_choices, family = binomial, data = train_world)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.39159  -0.25881   0.05221   0.37102   2.73990  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -12.9939     1.1976 -10.850  < 2e-16 ***
## GDP_per_capita       5.3478     0.7585   7.050 1.79e-12 ***
## Social_support       1.8516     0.5764   3.212  0.00132 ** 
## Life_expectancy      4.0083     0.9880   4.057 4.98e-05 ***
## Freedom_to_choices   7.7187     1.1764   6.561 5.34e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 851.08  on 613  degrees of freedom
## Residual deviance: 347.47  on 609  degrees of freedom
## AIC: 357.47
## 
## Number of Fisher Scoring iterations: 7

Estimated happiness score = -12.9939 + 5.3478 GDP per capita + 1.8516 Social support + 4.0083 Life expectancy + 7.7187 Freedom to choices

ß0 = -13.0128 suggests that the probability of happiness score not being 1 is 12.99%

ß1 = 5.3478 suggests that if GDP per capita increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(4.7810^{-4}\)%, holding all other variables constant. \(\frac{e^{(-12.9939+5.3478)}}{1+e^{(-12.9939+5.3478)}}=4.7810^{-4}\)%

ß2 = 1.8516 suggests that if the social support increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(1.4510^{-5}\)%, holding all other variables constant. \(\frac{e^{(-12.9939+1.8516)}}{1+e^{(-12.9939+1.8516)}}=1.4510^{-5}\)%

ß3 =4.0083 suggests that if the life expectancy increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(1.2510^{-4}\)%, holding all other variables constant. \(\frac{e^{(-12.9939+4.0083)}}{1+e^{(-12.9939+4.0083)}}=1.2510^{-4}\)%

ß4 =7.7187 suggests that if the freedom to choose increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(5.0910^{-3}\)%, holding all other variables constant. \(\frac{e^{(-12.9939+7.7187)}}{1+e^{(-12.9939+7.7187)}}=5.0910^{-3}\)%

Diagnosis

Regression diagnostics are a powerful ally in the model and data validation stages of regression model building. It allows us to assess the overall adequacy of our models. In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (glm.model.2) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(glm.world.2):

autoplot(glm.world.2, which = 1:2, nrow = 2, ncol = 1 )

Based on the Residuals vs Fitted plot, there are no distinct patterns for the residuals in the plot, hence we could come to a conclusion that there are no non-linear relationships in this model. In the Normal Q-Q plot, residuals deviate at the beginning and towards the end, it only falls along the line in the middle of the plot, hence this indicates that the residuals are not normally distributed and there is a possibility that there are extreme values.

ii) Multicollinearity

To identify if multicollinearity exists in our model (glm.world.2), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated.

vif(glm.world.2)

##     GDP_per_capita     Social_support    Life_expectancy Freedom_to_choices 
##           1.374020           1.050486           1.219158           1.297832

From the results, it is observed that all four independent variables (GDP_per_capita, Social_support, Life_expectancy and Freedom_to_choices) are not highly correlated as all variables have met the criterion cut-off point which is below 10. However, both scatter plot matrix and correlation matrix show that GDP_per_capita and Life_expectancy has a moderately strong positive correlation.

Prediction

confint(glm.world.2, level = 0.95)

## Waiting for profiling to be done...

##                          2.5 %     97.5 %
## (Intercept)        -15.5182230 -10.809450
## GDP_per_capita       3.9354312   6.916437
## Social_support       0.7413638   3.006032
## Life_expectancy      2.1309436   6.018335
## Freedom_to_choices   5.4933686  10.120224

The 95% confidence interval for GDP_per_capita and Social_support are -15.52 to -10.81 and 3.94 to 6.92 respectively. As it doesn’t include 0, hence we could conclude that both ß1 and ß2 are significant at the 5% level of significance. Also, Life_expectancy and Freedom_to_choices have an interval of 2.13 to 6.02 respectively, this has indicated that both ß3 and ß4 are significant at the 5% level of significance since 0 is not included.

Evaluation

cv_glm.world.2 <- train(
  form = Happy ~ GDP_per_capita + Social_support + Life_expectancy + Freedom_to_choices,
  data = train_world, 
  method = "glm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_glm.world.2

## Generalized Linear Model 
## 
## 614 samples
##   4 predictor
##   2 classes: '0', '1' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 553, 552, 553, 552, 553, 553, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.8615812  0.7229857

Through 10-fold cross validation, we could observe that this model has an accuracy of 0.8616, this indicates that the accuracy for this the correctly predicted observations to the total observations is 86.16%, which is considered high. Furthermore, the kappa value for the model is 0.7230. Kappa is measured to examine whether the collected data could represent the variables correctly. In this model, we could conclude that the collected data are reliable and could represent the variables correctly since its kappa value is relatively high.

Improving the Model

regfit.1.2 <- regsubsets(Happy ~ GDP_per_capita + Social_support + Life_expectancy + Freedom_to_choices, data = train_world)

par(mfrow=c(2,2))
plot(regfit.1.2, scale = "r2")
plot(regfit.1.2, scale = "adjr2")
plot(regfit.1.2, scale = "Cp")
plot(regfit.1.2, scale = "bic")

The best subset selection is performed to improve and find out the best fitted variables for the model. R-squared, adjusted R-squared, Cp and Bayesian Information Criterion have indicated that all 4 variables would be the best subset of predictors for this model, which includes GDP_per_capita, Social_support, Life_expectancy and Freedom_to_choices.

Therefore, we have concluded that three independent variables GDP_per_capita, Social_support, Life_expectancy and Freedom_to_choices would be the best predictors for our model.

Model 3 (glm.world3):

Our objective would be to examine the impact of the 4 significant variables, which are GDP_per_capita, Life_expectancy and Freedom_to_choices on Happiness_score, where 0 represents low happiness score (with the happiness score less than the mean) and 1 represents high happiness score. We only include the 3 variables that have the smallest p-value into this model.

Estimation

The estimated model (glm.world.3) was performed using the “glm” function in R using the logistic regression to identify the association of the binary independent variable with other dependent variables and fit the line values to the sigmoid curves.

glm.world.3 <- glm(Happy ~ GDP_per_capita + Life_expectancy + Freedom_to_choices, data = train_world, family = binomial)

Interpretation

summary(glm.world.3) 

## 
## Call:
## glm(formula = Happy ~ GDP_per_capita + Life_expectancy + Freedom_to_choices, 
##     family = binomial, data = train_world)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.49246  -0.29281   0.04748   0.37817   2.66679  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -12.0961     1.1318 -10.687  < 2e-16 ***
## GDP_per_capita       5.9186     0.7590   7.798 6.31e-15 ***
## Life_expectancy      4.3923     0.9702   4.527 5.98e-06 ***
## Freedom_to_choices   8.6361     1.1499   7.510 5.90e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 851.08  on 613  degrees of freedom
## Residual deviance: 358.36  on 610  degrees of freedom
## AIC: 366.36
## 
## Number of Fisher Scoring iterations: 7

Estimated happiness score = -12.0961 + 5.9186 GDP per capita + 4.3923 Life expectancy + 8.6361 Freedom to choices

ß0 = -12.0961 suggests that the probability of happiness score not being 1 is 12.10%

ß1 = 5.9186 suggests that if GDP per capita increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(2.0710^{-3}\)%, holding all other variables constant. \(\frac{e^{(-12.0961+5.9186)}}{1+e^{(-12.0961+5.9186)}}=5.0910^{-3}\)%

ß2 = 4.3923 suggests that if life expectancy increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(4.5110^{-4}\)%, holding all other variables constant. \(\frac{e^{(-12.0961+4.3923)}}{1+e^{(-12.0961+4.3923)}}=4.5110^{-4}\)%

ß3 =8.6361 suggests that if the freedom to choose increases by 1, the probability of having a happiness score of 1 increases multiplicatively by \(0.03\)%, holding all other variables constant. \(\frac{e^{(-12.0961+8.6361)}}{1+e^{(-12.0961+8.6361)}}=0.03\)%

Diagnosis

Regression diagnostics are a powerful ally in the model and data validation stages of regression model building. It allows us to assess the overall adequacy of our models. In this part, we will perform relevant diagnostics to detect underlying problems and assess if the estimated model (glm.model.3) works well for the data at hand, thus coming up with a much more improved model for prediction.

i) Residual Diagnostics

Residual plot and Q-Q plot for Model(glm.world.3):

autoplot(glm.world.3, which = 1:2, nrow = 2, ncol = 1 )

Based on the Residuals vs Fitted plot, we could see that the residuals are equally distributed along the line, hence there are no non-linear relationships in this model. In contrast, residuals only fall along the line in the middle, this could indicate that there might be extreme values in the data.

ii) Multicollinearity

To identify if multicollinearity exists in our model (glm.world.3), we apply the Variance Inflation Factor (VIF) to determine if the predictors in our multiple linear regression model are correlated.

vif(glm.world.3)

##     GDP_per_capita    Life_expectancy Freedom_to_choices 
##           1.393040           1.209211           1.298513

From the results, it is observed that all three independent variables (GDP_per_capita, Life_expectancy and Freedom_to_choices) are not highly correlated as all variables have met the criterion cut-off point which is below 10. However, both scatter plot matrix and correlation matrix show that GDP_per_capita and Life_expectancy has a moderately strong positive correlation.

Prediction

confint(glm.world.3, level = 0.95)

## Waiting for profiling to be done...

##                         2.5 %     97.5 %
## (Intercept)        -14.486219 -10.036201
## GDP_per_capita       4.503541   7.487018
## Life_expectancy      2.554101   6.370489
## Freedom_to_choices   6.471742  10.993504

The 95% confidence interval for GDP_per_capita and Life_expectancy are 4.50 to 7.49 and 2.55 to 6.37 respectively. As it doesn’t include 0, hence we could conclude that both ß1 and ß2 are significant at the 5% level of significance. Also, Freedom_to_choices has an interval of 6.47 to 10.99, this has indicated that both ß3 are significant as well at the 5% level of significance since 0 is not included.

Evaluation

cv_glm.world.3 <- train(
  form = Happy ~ GDP_per_capita + Life_expectancy + Freedom_to_choices,
  data = train_world, 
  method = "glm",
  trControl = trainControl(method = "cv", number = 10)
)
cv_glm.world.3

## Generalized Linear Model 
## 
## 614 samples
##   3 predictor
##   2 classes: '0', '1' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 553, 552, 553, 551, 553, 553, ... 
## Resampling results:
## 
##   Accuracy  Kappa    
##   0.856891  0.7137284

Through 10-fold cross validation, we could observe that this model has an accuracy of 0.8569, this indicates that the accuracy for this the correctly predicted observations to the total observations is 85.69%, which is considered high. Furthermore, the kappa value for the model is 0.7137. Kappa is measured to examine whether the collected data could represent the variables correctly. In this model, we conclude that the collected data could represent the variables correctly and it is highly reliable since its kappa value is relatively high.

Improving the Model

regfit.1.3 <- regsubsets(Happy ~ GDP_per_capita + Life_expectancy + Freedom_to_choices, data = train_world)

par(mfrow=c(2,2))
plot(regfit.1.3, scale = "r2")
plot(regfit.1.3, scale = "adjr2")
plot(regfit.1.3, scale = "Cp")
plot(regfit.1.3, scale = "bic")

The best subset selection is performed to improve and find out the best fitted variables for the model. R-squared, adjusted R-squared, Cp and Bayesian Information Criterion have indicated that GDP_per_capita, Life_expectancy and Freedom_to_choices would be the best subset of predictors for this model.

Therefore, we could conclude that three independent variables, GDP_per_capita, Social_support, Life_expectancy and Freedom_to_choices would be the best predictors for our model.

Best model selection :

Comparing goodness of fit between models :

Residual deviance is used to examine if the model is a good fit for the collected data, the smaller the residual deviance, the better the model fits the data. From the table above, it is observed that both Model (glm.world) and Model (glm.world.2) have the two lowest residual deviance, which indicates that the data in these two models fit well. Besides that, based on the AIC value, we can see that Model (glm.world.2) has the lowest AIC value as compared to other models, which indicates that the model fits the data well and the penalisation is the lowest for Model (glm.world.2).

In addition, both Model (glm.world.2) and Model (glm.world.3) have a higher confusion matrix accuracy than Model (glm.world.3) which are 86.64% and 86.81% respectively. This indicates that the tested accuracy from predicting is performing well as the accuracy is quite high. Lastly, both Model (glm.world) and Model (glm.world.2) have the highest AUC value, this suggests that both models are performing well in distinguishing between the high happiness score and low happiness score since the higher the AUC value, the higher the degree of separability.

Therefore, we can conclude that Model (glm.world.2) is our best model as it not only has the smallest residual deviance and AIC value but also the highest confusion matrix accuracy and AUC value. These metrics have indicated the data could fit Model (glm.world.2) well and could highly distinguish the high happiness score and low happiness score.

Ridge Regression

Reasons for using ridge regression:

Although there are no presence of multicollinearity from all of our models based on the VIF value we have calculated earlier. It is observed that the VIF value for Model (lm.top10) is the highest as compared to other models from other research question.

Besides that, based on the VIF values obtained from Model (lm.top10), we have found that both VIF values for ‘Life_expectancy’ and ‘Freedom_to_choices’ are slightly higher than the rest. This shows that there is a small chance of multicollinearity in this model where both variables might have a correlation with the other variables. Not to mention, we have found that the R-squared value for Model (lm.top10) is substantially higher as compared to other models during our best model selection in Research Question 1, this may indicate that Model (lm.top10) is overfitted. Therefore, we would like to identify if our model will be improved by applying ridge regression to our model.

Firstly, we will split the data set into 80% training set and 20% test set.

set.seed(123)

# Split data into training and testing set
split <- sample(1:nrow(lm_top10), round(nrow(lm_top10) * 0.8)) # 80-20 split
train_top10_2 <- lm_top10[split, ]
test_top10_2 <- lm_top10[-split, ]

Next, we have to standardize the training set and test set because ridge regression assumes that the predictors are standardized and the response is centered. Therefore standardization will be required to avoid HSK-ridge regression.

train_x <- model.matrix(Happiness_score ~., train_top10_2)[,-1]  # discard intercept
train_y <- train_top10_2$Happiness_score

test_x <- model.matrix(Happiness_score ~., test_top10_2)[,-1]
test_y <- test_top10_2$Happiness_score

Furthermore, we will create a grid of values ranging from lambda 10^10 to 10^−2 to cover a full range of scenarios from the null model (only intercept) to a full model (OLS fit). This is because the glmnet package only performs a ridge regression for an automatically selected range of λ values by default.

grid <- 10^seq(10, -2, length = 100)

# Apply ridge regression to the data
ridge_reg <- glmnet(x = train_x, y = train_y, alpha = 0, lambda = grid)

# The higher the lambda, the more the coefficients are shrinking towards zero.
plot(ridge_reg, xvar= "lambda", label = TRUE, main = "Ridge penalty")
legend("topright", lwd = 1, col = 1:6, legend = colnames(train_x), cex = .7)

λ is a tuning parameter that allows us to control our model from overfitting to the training data. We will perform the 10-fold cross-validation to identify the optimal λ value.

# Perform ridge regression with 10-fold CV
ridge_reg_cv <- cv.glmnet(
  x = train_x,
  y = train_y,
  alpha = 0        # alpha = 0 for ridge penalty
)

# Plot the results of penalty versus CV MSE
plot(ridge_reg_cv,  main = "MSE of a range of 10-fold cv-lambda")

Based on the figure above, the first vertical dashed line represents the λ value with the minimum MSE, while the second vertical dashed lines represents the largest λ value within one standard error of the minimum MSE. We can see that we can minimize our MSE by applying approximately -3.8 ≤ log(λ) ≤ −1.

# select optimal cv-lambda with the minimum MSE
ridge.pred <- predict(ridge_reg, s = ridge_reg_cv$lambda.min, newx = test_x) # unstandardized best cv-lambda value

Mean squared error (MSE) is a measure of how close a regression line is to a set of points. The lower the mean squared error, the better the model fits the data.

# calculate test MSE
mean((ridge.pred - test_y)^2) # test MSE of best cv-lambda

## [1] 0.01657888

# fit ridge regression model on full dataset
final.mod <- glmnet(train_x, train_y, alpha = 0)

# display coefficient using lambda chosen by cv
ridge.coef <- predict(final.mod, type = "coefficients", s = ridge_reg_cv$lambda.min)[1:6, ] 
ridge.coef

##        (Intercept)     GDP_per_capita     Social_support    Life_expectancy 
##         6.71894800         0.38034972        -0.01564159        -0.11875855 
## Freedom_to_choices         Generosity 
##         0.76863135        -0.49792503

We can see that all variables remain in the model with some penalty being applied to them as ridge regression does not perform variable selection where none of the estimated coefficients are zero.

# CV-RMSE using 1-SE rule lambda (1-standard error of minimum MSE)
ridge_reg_cv$cvm[ridge_reg_cv$lambda == ridge_reg_cv$lambda.1se] 

## [1] 0.0119372

# The lambda for the MSE using 1-SE rule lambda
ridge_reg_cv$lambda.1se

## [1] 0.3887182

ridge_reg_min <- glmnet(
  x = train_x,
  y = train_y,
  alpha = 0
)

plot(ridge_reg_min, xvar = "lambda")
legend("topright", lwd = 1, col = 1:6, legend = colnames(train_x), cex = .7)
abline(v = log(ridge_reg_cv$lambda.1se), col = "red", lty = "dashed")

We have plot the coefficients across the λ values and the vertical dashed red line represents the largest λ value within one standard error of the minimum MSE. This allows us to identify how much we can apply constraints to our coefficients while maximizing the predictive accuracy of our model.

Based on the plot above, we can see that the top influential variables in order are ‘GDP_per_capita’, ‘Freedom_to_choices’, ‘Generosity’, followed by ‘Life_expectancy’ and ‘Social_support’. The higher the coefficient value, the higher the influence.

Lasso Regression

Reasons for using lasso regression:

Given that Model (lm.top10) is the best in predicting the Happiness_score for the top 10 happiest countries, where it has used all 5 variables which are GDP_per_capita, Social_support, Life_expectancy, Freedom_to_choices and Generosity as shown earlier. We would like to identify if we can further automate certain parts of model selection such as the variable selection or parameter elimination where it will arbitrality select features which highly correlates to Happiness_score and reduced the coefficients of the rest to zero.

# Apply lasso regression to the data
lasso_reg <- glmnet(x = train_x, y = train_y, alpha = 1)

# The higher the lambda, the more the coefficients are shrinking towards zero.
plot(lasso_reg, xvar= "lambda", label = TRUE)
legend("topright", lwd = 1, col = 1:6, legend = colnames(train_x), cex = .7)

We will perform the 10-fold cross-validation to identify the optimal λ value for our lasso regression.

# Perform lasso regression with 10-fold CV
lasso_reg_cv <- cv.glmnet(
  x = train_x,
  y = train_y,
  alpha = 1        # the alpha = 1 for lasso penalty
)

# Plot the results of penalty versus CV MSE
plot(lasso_reg_cv)

Based on the plot above, we can see that we can minimize our MSE by applying approximately -4.8 ≤ log(λ) ≤ −3.8.

# select optimal cv-lambda with the minimum MSE
lasso.pred <- predict(lasso_reg, s = lasso_reg_cv$lambda.min, newx = test_x) # manually set s 
mean((lasso.pred - test_y)^2) # test MSE of best cv-lambda

## [1] 0.01568489

# display coefficient using lambda chosen by cv
final.mod2 <- glmnet(train_x, train_y, alpha = 1)
lasso.coef <- predict(final.mod2, type = "coefficients", s = lasso_reg_cv$lambda.min)[1:6, ] 
lasso.coef

##        (Intercept)     GDP_per_capita     Social_support    Life_expectancy 
##         6.60547825         0.39246383         0.00000000        -0.03367441 
## Freedom_to_choices         Generosity 
##         0.74920307        -0.46879125

We can see that lasso regression has reduced the estimated coefficients of some of the less important variables to zero as lasso regression performs variable selection unlike ridge regression. Hence, some variables will only remain in the model.

# CV-RMSE using 1-SE rule lambda (1-standard error of minimum MSE)
lasso_reg_cv$cvm[lasso_reg_cv$lambda == lasso_reg_cv$lambda.1se] 

## [1] 0.01229725

# The lambda for the MSE using 1-SE rule lambda
lasso_reg_cv$lambda.1se

## [1] 0.02330304

lasso_reg_min <- glmnet(
  x = train_x,
  y = train_y,
  alpha = 1
)

plot(lasso_reg_min, xvar = "lambda")
legend("topright", lwd = 1, col = 1:6, legend = colnames(train_x), cex = .7)
abline(v = log(lasso_reg_cv$lambda.1se), col = "red", lty = "dashed")

lasso_reg_min <- glmnet(
  x = train_x,
  y = train_y,
  alpha = 1
)

plot(lasso_reg_min, xvar = "lambda")
legend("topright", lwd = 1, col = 1:6, legend = colnames(train_x), cex = .7)
abline(v = log(lasso_reg_cv$lambda.1se), col = "red", lty = "dashed")

Based on the plot above, we can see that only the estimated coefficients for ‘Freedom_to_choices’, ‘GDP_per_capita’ and ‘Generosity’ are not penalize to zero, while the estimated coefficients for ‘Social_support’ and ‘Life_expectancy’ has been penalized to zero in order to maximize the predictive accuracy of our model.

Based on the plot above, we can see that the top most influential variables are ‘Freedom_to_choices’ followed by ‘GDP_per_capita’ and ‘Generosity’.

Insights drawn from both methods :

By comparing between ridge regression and lasso regression, we can see that lasso regression performs better as the MSE value for lasso regression is 0.0157 which is lower as compared to ridge regression with a MSE value of 0.0166. Therefore, we do not have to retain all variables in our model and also reduce the noise caused by less influential variables and minimize the chances of multicollinearity through ridge regression in this case.

Ultimately, we can see that the variables which strongly correlates with the ‘Happiness_score’ for the top 10 happiest countries are only ‘Freedom_to_choices’, ‘GDP_per_capita’ and ‘Generosity’ from our lasso regression as compared to using all 5 variables as a predictor for ‘Happiness_score’ in our Model (lm.top10) using a linear regression method.