GroupProject2.Rmd

---
title: "Group Project 2"
subtitle: <center> Intro and EDA </center>
author: <center> Joshua Carpenter, Yong-Nan Chan, Andy Phillips, and Brandon Fletcher <center>
output: html_document
---

<style type="text/css">
h1.title {
font-size: 40px;
text-align: center;
}
h4.author 
font-size: 40px;
text-align: center;
}
</style>

```{r setup, include=FALSE}
# load packages here
library(corrplot)
library(tidyverse)
library(ggfortify)
library(car)
library(bestglm)
library(glmnet)
library(gridExtra)
set.seed(89954)

# Useful functions for later
resid_vs_fitted <- function(model) {
  autoplot(model, which = 1, ncol = 1) +
    theme_minimal() +
    theme(aspect.ratio = 1)
}
jcreg_qq <- function(model) {
  autoplot(model, which = 2, ncol = 1) +
    theme_bw() +
    theme(aspect.ratio = 1)
}
jcreg_hist <- function(model) {
  residuals <- data.frame(residuals = resid(model))
  ggplot(data = residuals, mapping = aes(x = residuals)) +
    geom_histogram(binwidth = sd(residuals$residuals / 4),
                   mapping = aes(y = ..density..)) +
    stat_function(fun = dnorm,
                  color = "blue",
                  args = list(mean = 0,
                              sd = sd(residuals$residuals)),
                  size = 1.2) +
    xlab("Residuals") +
    ylab("Density") +
    theme_light()
}
jcreg_av <- function(model) {
  predictors <- attr(model$terms, "term.labels")
  rows <- ifelse(length(predictors > 3),
                 floor(sqrt(length(predictors))),
                 1)
  cols <- ifelse(length(predictors > 3),
                 length(predictors) / rows,
                 3)
  par(pty = "s")
  car::avPlots(model, layout = c(rows, cols), pch = 19)
}
point_matrix <- function(data) {
  par(pty = "s", las = 1)
  pairs(data, pch = ".", lower.panel = NULL)
}
rpred_col <- function(data, residuals, predictor) {
  ggplot2::ggplot(data = data,
         mapping = ggplot2::aes(x = pull(data, predictor),
                       y = residuals)) +
    ggplot2::geom_point() +
    ggplot2::geom_smooth(se = FALSE, span = 0.95, n = 7, size = 0.5) +
    ggplot2::geom_abline(slope = 0, intercept = 0, linetype = "dashed") +
    ggplot2::theme_minimal() +
    ggplot2::theme(aspect.ratio = 1) +
    ggplot2::xlab(predictor) +
    ggplot2::ylab("Residuals")
}
resid_vs_pred <- function(model) {
  data <- model.frame(model)
  predictors <- attr(model$terms, "term.labels")
  plots <- lapply(predictors, rpred_col, data = data, residuals = resid(model))
  plots["ncol"] <- ceiling(sqrt(length(plots)))
  plots["top"] <- "Residuals vs Predictors"
  do.call(gridExtra::grid.arrange, plots)
}
jcreg_boxplot <- function(model) {
  residuals <- data.frame(residuals = resid(model))
  ggplot2::ggplot(data = residuals, mapping = ggplot2::aes(y = residuals)) +
    ggplot2::geom_boxplot() +
    ggplot2::stat_summary(mapping = ggplot2::aes(x = 0),
                 fun = mean, geom = "point",
                 shape = 4, size = 2, color = "darkred") +
    ggplot2::theme_classic() +
    ggplot2::theme(aspect.ratio = 2,
          axis.text.x = ggplot2::element_blank(),
          axis.ticks.x = ggplot2::element_blank()) +
    #  scale_y_continuous(limits = c(-20000, 30000), breaks = seq(-20000, 30000, 10000)) +
    ggplot2::ylab("Residuals") +
    ggplot2::xlab("")
}
jcreg_cooksd <- function(model, nLabels = 3) {
  cooks_d <- cooks.distance(model)
  top_cd <- as.numeric(names(sort(cooks_d, decreasing = TRUE)[1:nLabels]))

  ggplot2::ggplot() +
    ggplot2::geom_point(data = tibble::tibble(cooks_d),
               mapping = ggplot2::aes(x = as.numeric(names(cooks_d)),
                             y = cooks_d)) +
    ggplot2::geom_text(mapping = ggplot2::aes(x = top_cd,
                            y = cooks_d[top_cd] + max(cooks_d) / 40,
                            label = top_cd)) +
    ggplot2::theme_bw() +
    ggplot2::ylab("Cook's Distance") +
    ggplot2::xlab("Observation Number") +
    ggplot2::geom_hline(mapping = ggplot2::aes(yintercept = 4 / length(cooks_d)),
               color = "red", linetype = "dashed") +
    ggplot2::theme(aspect.ratio = 1)
}
dfb_col <- function(df_betas, predictor, nLabels = 3) {
  require(tibble)
  # Find which observations have the highest dfbetas
  top_vals <- df_betas[predictor] %>%
    arrange(desc(abs(eval(parse(text = predictor))))) %>%
    .[1:nLabels,] %>%
    pull(predictor)
  top_ind <- which(pull(df_betas, predictor) %in% top_vals)

  out <- ggplot2::ggplot() +
    ggplot2::geom_point(data = df_betas,
               mapping = ggplot2::aes(x = as.numeric(rownames(df_betas)),
                             y = abs(pull(df_betas, predictor)))) +
    ggplot2::geom_text(mapping = ggplot2::aes(x = top_ind,
                            y = abs(pull(df_betas, predictor)[top_ind]) + 0.07,
                            label = top_ind)) +
    ggplot2::theme_bw() +
    ggplot2::theme(aspect.ratio = 1) +
    ggplot2::ylab("Abs of DFBETAS") +
    ggplot2::xlab("Observation Number") +
    ggtitle(predictor)

  if(length(dfbetas) <= 30) {
    out <- out +
      ggplot2::geom_hline(mapping = ggplot2::aes(yintercept = 1),
                 color = "red", linetype = "dashed")
  }else {
    out <- out +
      ggplot2::geom_hline(mapping = ggplot2::aes(yintercept = 2 / sqrt(length(dfbetas))),
                 color = "red", linetype = "dashed")
  }
  return(out)
}
jcreg_dfbetas <- function(model, nLabels = 3) {
  predictors <- attr(model$terms, "term.labels")
  df_betas <-  tibble::as_tibble(dfbetas(model))[, predictors]

  plots <- lapply(predictors, dfb_col, df_betas = df_betas, nLabels = nLabels)
  plots["ncol"] <- ceiling(sqrt(length(plots)))
  do.call(gridExtra::grid.arrange, plots)
}
jcreg_dffits <- function(model, nLabels = 3) {
  df_fits <- dffits(model)
  top_dff <- as.numeric(names(sort(abs(df_fits), decreasing = TRUE)[1:nLabels]))

  df_fits_plot <- ggplot2::ggplot() +
    ggplot2::geom_point(data =tibble::tibble(df_fits),
               mapping = ggplot2::aes(x = as.numeric(names(df_fits)),
                             y = abs(df_fits))) +
    ggplot2::geom_text(mapping = ggplot2::aes(x = top_dff,
                            y = abs(df_fits[top_dff]) + max(df_fits) / 40,
                            label = top_dff)) +
    ggplot2::theme_bw() +
    ggplot2::ylab("Absolute Value of DFFITS for Y") +
    ggplot2::xlab("Observation Number") +
    ggplot2::theme(aspect.ratio = 1)
  if(length(df_fits) <= 30) {
    df_fits_plot +
      ggplot2::geom_hline(mapping = ggplot2::aes(yintercept =
                                 2 * sqrt(length(model$coefficients) /
                                            length(df_fits))),
                 color = "red", linetype = "dashed")
  }else {
    df_fits_plot +
      ggplot2::geom_hline(mapping = ggplot2::aes(yintercept = 1),
                 color = "red", linetype = "dashed")
  }
}
cor_graphic <- function(data, show_key = FALSE, title = TRUE) {
  par(mfrow = c(1, 2))
  corrplot::corrplot(cor(data), method = "number", type = "upper", diag = F,
                     tl.col = "#1f3366", cl.pos = "n")
  if(title) title("Correlation Coefficients")
  if (show_key) {
    corrplot::corrplot(cor(data), type = "upper", diag = F, tl.col = "#1f3366")
  }
  else {
    corrplot::corrplot(cor(data), type = "upper", diag = F, tl.col = "#1f3366",
                       cl.pos = "n")
  }
  if(title) title("Correlation Matrix")
}
```

# Background and Introduction

The World Health Organization (WHO) and United Nations (UN) were interested in factors affecting life expectancy in countries around the world. In 2015 they collected demographic and immunization data from 130 countries.\

The purpose of this analysis is to find which factors are correlated with life expectancy and their relationship. We expect to find that having a greater proportion of the population immunized will lead to increased life expectancy and that greater schooling will also be associated with greater life expectancy, although especially in that case we do not believe there is necessarily a causal relationship.\

To analyze these hypotheses, we will perform linear regression using `life.expectancy` as our response variable. We will check the assumptions of linear regression and remove or transform variables as necessary. If transformations are necessary, we will use the Box-Cox method and plots of different transformations to help us determine the best one. In order to avoid multicolinearity and over-fitting, we will use several variable selection and shrinkage methods to choose the most appropriate subset of variables to keep in the model. After performing regression, we will examine the model and the hypothesis tests to see which variables are significant and how they affect life expectancy. Finally, we will provide our results along with confidence intervals.

# Methods and Results

The data set we will be using contains measurements, collected in 2015, of 130 different countries. This was collected by the World Health Organization and the United Nations. The data set contains measurements for both developed and developing countries.

The following table displays the variable names in this data set, along with their descriptions.

Variable       | Description
-------------- | -------------
Status         | Country Status (Developed or Developing)
Life.expectancy| Average life expectancy in years
Adult.Mortality| Probability of dying between 15 and 60 years per 1000 population
Hepatitis.B    | Immunization coverage among 1-year-olds (%)
Measles        | Number of reported cases er 1000 population
BMI            | Average Body Mass Index of entire population
Polio          | Immunization coverage among 1-year-olds (%)
Diphtheria     | Immunization coverage among 1-year-olds (%)
GDP            | Gross Domestic Product per Capita (In USD)
Population     | Population of the country
Schooling      | Average number of years of schooling

## Exploratory Statistics

We start by applying basic summary and exploratory statistics to this data to better understand the data and identify trends.

```{r, message=FALSE}
life_expect <- read_csv("lifeexpect.csv") %>%
  # Change Status into a numeric variable where 1 means
  #   developed and 0 means developing
  mutate(Status = factor(Status),
         # Fix capitalization on Life.expectancy
         Life.Expectancy = Life.expectancy) %>%
  # Remove un-needed columns
  select(Life.Expectancy, Adult.Mortality, Hepatitis.B, BMI, Polio,
         Diphtheria, GDP, Population, Schooling, Status) %>% 
  filter(Diphtheria < 100)
life_expect
summary(life_expect)
```
```{r, fig.width=8, fig.height=8}
# create data set only with continuous variable
life_expect_cont <- life_expect %>% select(-Status)

### scatterplot matrix (only with continuous)
point_matrix(life_expect_cont)
```
```{r}
### correlation matrix (only used for continuous variables)
round(cor(life_expect_cont), 2)
corrplot(cor(life_expect_cont), type = "upper", diag = F,
         tl.col = "#1f3366", cl.pos = "n")
```

```{r, results=FALSE, message=FALSE}
#### Histogram
histLifeExpect <- function(variable, name, width) {
  ggplot(data = life_expect, mapping = aes(x = variable)) +
  geom_histogram(mapping = aes(y = ..density..), binwidth = width) +
  xlab(name) +
  ylab("Density") +
  theme_classic() +
  theme(aspect.ratio = 1)
}
# We remove GDP because the hist is useless and takes a long time to run
to_hist <- life_expect_cont %>%
  select(-c(Life.Expectancy, GDP))
hists <- mapply(histLifeExpect,
       variable = to_hist,
       name = names(to_hist),
       width = c(25, 5, 5, 10, 7, 10000000, 2),
       SIMPLIFY = FALSE)
do.call(grid.arrange, hists)

# Individual scatterplots for easier viewing
scatter <- function(data, x_ind, y_ind) {
  x <- pull(data, x_ind)
  y <- pull(data, y_ind)
  ggplot(mapping = aes(x = x, y = y)) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    theme_minimal() +
    theme(aspect.ratio = 1) +
    xlab(names(data)[x_ind]) +
    ylab(names(data)[y_ind])
}
splots <- map(2:length(life_expect_cont), scatter, data = life_expect_cont, y_ind = 1)
do.call(grid.arrange, splots)

### Box Plot (for categorical: Status)
ggplot(data = life_expect, mapping = aes(x = Status, y = Life.Expectancy)) +
  geom_boxplot() +
  theme_bw() +
  theme(aspect.ratio = 1)
```

## Variable Selection
```{r, fig.align='center', results=FALSE, message = FALSE}
#create new variable and put Life.Expectancy last
life_expect_var <- life_expect %>%
  select(Adult.Mortality : Status, Life.Expectancy) %>%
  as.data.frame()
head(life_expect_var)

life_expect_var_lm <- lm(Life.Expectancy ~ ., life_expect_var)
summary(life_expect_var_lm)

#BEST SUBSETS
life_best_subsets_bic <- bestglm(life_expect_var, IC = "BIC", method = "exhaustive")

life_best_subsets_bic$BestModels
summary(life_best_subsets_bic$BestModel)

# BACKWARD
life_backward_bic <- bestglm(life_expect_var, IC = "BIC", method = "backward", t = 100)
summary(life_backward_bic$BestModel)

# SEQUENTIAL
life_seqrep_bic <- bestglm(life_expect_var, IC = "BIC", method = "seqrep", t = 100)
summary(life_seqrep_bic$BestModel)


# Note: had trouble with the data turning into a matrix here
life_expect_var_x <- as.matrix(life_expect_var[,1:7])
life_expect_var_y <- life_expect_var[,8]

#LASSO
life_expect_var_lasso_cv <- cv.glmnet(x = life_expect_var_x,
                              y = life_expect_var_y,
                              type.measure = "mse",
                              alpha = 1)
life_expect_var_lasso_cv$lambda.min
life_expect_var_lasso_cv$lambda.1se

coef(life_expect_var_lasso_cv, s = "lambda.min")
coef(life_expect_var_lasso_cv, s = "lambda.1se")

#ELASTIC NET
life_expect_var_elastic_cv <- cv.glmnet(x = life_expect_var_x,
                              y = life_expect_var_y,
                              type.measure = "mse",
                              alpha = 0.5)
life_expect_var_elastic_cv$lambda.min
life_expect_var_elastic_cv$lambda.1se

coef(life_expect_var_elastic_cv, s = "lambda.min")
coef(life_expect_var_elastic_cv, s = "lambda.1se")

```


Variable            | Best Subset | Backward | Sequential Replacement | LASSO  | Elastic Net
--------------------| ----------- | -------- | ---------------------- | ------ | -----------
  Adult.Mortality   |      X      |     X    |          X             |   X    |      X
  Hepatitis.B       |             |          |                        |   X    |      X
  BMI               |             |          |                        |   X    |      X
  Polio             |             |          |                        |   X    |      X
  Diphtheria        |      X      |     X    |          X             |   X    |      X
  GDP               |             |          |                        |   X    |      X
  Population        |             |          |                        |        |
  Schooling         |             |          |                        |        |
  Status            |             |          |          X             |        |


Given the results from all of the variable selection procedures, shown in the table above, we choose to keep Adult.Mortality and Diphtheria.

## Initial Linear Model
Next, we will run an initial linear model and take a look at the residuals to see what transformations might need to be done.
```{r, warning=FALSE, fig.width=9.5}
init_model <- lm(Life.Expectancy ~ Adult.Mortality + Diphtheria, data = life_expect)
rvf <- resid_vs_fitted(init_model)@plots[[1]] +
  theme_bw() + theme(aspect.ratio = 1)
hist <- jcreg_hist(init_model) +
  theme_bw() + theme(aspect.ratio = 1) +
  ggtitle("Histogram of Residuals")
qq <- jcreg_qq(init_model)@plots[[1]]
grid.arrange(rvf, hist, qq, ncol = 3)
shapiro.test(init_model$residuals)
```

The residuals are fairly normal, with some outliers at the lower tail. There does appear to be a slight cone shape to the residuals, which could be a problem. Based on these graphs we will attempt a transformation on the response.

## Check if we need to include an interaction
```{r}
interaction_model <- lm(Life.Expectancy ~ Adult.Mortality + Diphtheria + Adult.Mortality:Diphtheria, data = life_expect)
anova(init_model, interaction_model)
```

The p-value for the difference between models is highly insignificant, so we will not include an interaction.

## Transformations
```{r}
# Apply boxCox to give us a starting point
boxCox(init_model, lambda = seq(-2, 4, 1/10))
```

```{r, message=FALSE, fig.width=9.5}
# Squared transformation
life_expect_sq <- life_expect %>%
  mutate(Life.Expectancy = Life.Expectancy^2)
sq_model <- lm(Life.Expectancy ~ Adult.Mortality + Diphtheria, data = life_expect_sq)
rvf <- resid_vs_fitted(sq_model)@plots[[1]] +
  theme_bw() + theme(aspect.ratio = 1)
hist <- jcreg_hist(sq_model) +
  theme_bw() + theme(aspect.ratio = 1) +
  ggtitle("Histogram of Residuals")
qq <- jcreg_qq(sq_model)@plots[[1]]
grid.arrange(rvf, hist, qq, ncol = 3)
shapiro.test(sq_model$residuals)

# Cubed transformation
life_expect_cu <- life_expect %>%
  mutate(Life.Expectancy = Life.Expectancy^3)
cu_model <- lm(Life.Expectancy ~ Adult.Mortality + Diphtheria, data = life_expect_cu)
rvf <- resid_vs_fitted(cu_model)@plots[[1]] +
  theme_bw() + theme(aspect.ratio = 1)
hist <- jcreg_hist(cu_model) +
  theme_bw() + theme(aspect.ratio = 1) +
  ggtitle("Histogram of Residuals")
qq <- jcreg_qq(cu_model)@plots[[1]]
grid.arrange(rvf, hist, qq, ncol = 3)
shapiro.test(cu_model$residuals)
```

We tried a squared, power of 1.5, cubed, log, and square root transformation. Based on the transformations attempted, the best one appears to be a squared or cubic transformation. The cubic transformation sees great improvement in homoscedasticity and some improvement in normality. The squared transformation has somewhat worse homoscedasticity the the cubic, but better normality. For simplicity, we will complete the rest of the analysis on a squared scale. We will attempt to further improve homoscedasticity with transformations of the predictors.

```{r, message=FALSE, results=FALSE}
life_expect_sq <- life_expect_sq %>%
  select(Life.Expectancy, Adult.Mortality, Diphtheria)
to_hist <- life_expect_sq %>%
  select(-Life.Expectancy)
hists <- mapply(histLifeExpect,
       variable = to_hist,
       name = names(to_hist),
       width = c(25, 5),
       SIMPLIFY = FALSE)
hists[["ncol"]] <- 2
do.call(grid.arrange, hists)

# Log transformation
life_expect_t <- life_expect_sq %>%
  mutate(Adult.Mortality = log(Adult.Mortality),
        Diphtheria = log(Diphtheria))
to_hist <- life_expect_t %>%
  select(-Life.Expectancy)
hists <- mapply(histLifeExpect,
       variable = to_hist,
       name = names(to_hist),
       width = c(0.3, 0.07),
       SIMPLIFY = FALSE)
hists[["ncol"]] <- 2
do.call(grid.arrange, hists)

# Sqrt transformation
life_expect_t <- life_expect_sq %>%
  mutate(Adult.Mortality = sqrt(Adult.Mortality),
        Diphtheria = log(Diphtheria))
to_hist <- life_expect_t %>%
  select(-Life.Expectancy)
hists <- mapply(histLifeExpect,
       variable = to_hist,
       name = names(to_hist),
       width = c(1, 0.1),
       SIMPLIFY = FALSE)
hists[["ncol"]] <- 2
do.call(grid.arrange, hists)

# Cbrt transformation
life_expect_t <- life_expect_sq %>%
  mutate(Adult.Mortality = Adult.Mortality^(1/3),
        Diphtheria = log(Diphtheria))
to_hist <- life_expect_t %>%
  select(-Life.Expectancy)
hists <- mapply(histLifeExpect,
       variable = to_hist,
       name = names(to_hist),
       width = c(0.3, 0.07),
       SIMPLIFY = FALSE)
hists[["ncol"]] <- 2
do.call(grid.arrange, hists)
```

Adult Mortality was skewed, so we attempted a transformation. Clearly the log transform looks more normal than the others. For Diphtheria, we will stick to the original data.
```{r}
life_expect_t <- life_expect_sq %>%
  mutate(Adult.Mortality = log(Adult.Mortality))
```

## Final Linear Model

```{r, fig.align='center'}
life_expect_lm <- lm(Life.Expectancy ~ ., data = life_expect_t)
summary(life_expect_lm)
```

Our final fitted model is: 

$$(\widehat{\text{Life.Expectancy}}_i)^2=-567.8-210.6\cdot\log(\text{Adult.Mortality}_i)+86.0\cdot\text{Diphtheria}_i$$

```{r, fig.height=7, fig.width=7}
point_matrix(life_expect_t)
```
```{r}
jcreg_av(life_expect_lm)
cor_graphic(life_expect_t, title = FALSE)
vif(life_expect_lm)
```

The linearity assumption is met. The partial regression plots appear to be linear. 
There is a strong negative correlation between Adult.Mortality and the Percentage of 1 year olds vaccinated against Diptheria. However, this correlation is below the multicollinearity cutoff of 0.8. The average VIF is 1.94 which is above the cutoff of 1. However we feel that the multicollinearity assumption is met.
The independence assumption is met, the data was collected randomly.

```{r, fig.width=9.5}
# Plot residuals
rvf <- resid_vs_fitted(life_expect_lm)@plots[[1]] +
  theme_bw() + theme(aspect.ratio = 1)
hist <- jcreg_hist(life_expect_lm) +
  theme_bw() + theme(aspect.ratio = 1) +
  ggtitle("Histogram of Residuals")
qq <- jcreg_qq(life_expect_lm)@plots[[1]]
boxp <- jcreg_boxplot(life_expect_lm)
grid.arrange(rvf, hist, qq, boxp, ncol = 2)
```
```{r, message=FALSE}
resid_vs_pred(life_expect_lm)
shapiro.test(life_expect_lm$residuals)
```
The Normality assumption is met. The residuals appear to be fairly normally distributed around 0. The p-value of the shapiro-wilk test is greater than 0.05 so we can say that the data is normally distributed.

There does not appear to be any distinctive cone shape to the data in the residual vs. fitted value plot so we can assume that the homoskedastic assumption is met.
```{r}
# Check for influential points
jcreg_dfbetas(life_expect_lm)
dffits <- jcreg_dffits(life_expect_lm)
cooksd <- jcreg_cooksd(life_expect_lm)
grid.arrange(dffits, cooksd, ncol = 2)
```
There are a few points (110 and 117) that show up in a few of the plots but these are not too far out from the other points so the assumption of no influential points is met

## Confidence Intervals

```{r, fig.align='center'}

confint(life_expect_lm, level = 0.95)

#confidence interval with mean adult mortality and mean diphtheria
predict(life_expect_lm,
        newdata = data.frame(Adult.Mortality = log(169),
                             Diphtheria = 72),
        interval = "confidence",
        level = 0.95)

#prediction interval with mean adult mortality and mean diphtheria
predict(life_expect_lm,
        newdata = data.frame(Adult.Mortality = log(169),
                             Diphtheria = 72),
        interval = "prediction",
        level = 0.95)

#confidence interval with mean adult mortality and everyone vaccinated against Diphtheria
predict(life_expect_lm,
        newdata = data.frame(Adult.Mortality = log(169),
                             Diphtheria = 99),
        interval = "confidence",
        level = 0.95)

#confidence interval with mean adult mortality and no one vaccinated against Diphtheria
predict(life_expect_lm,
        newdata = data.frame(Adult.Mortality = log(169),
                             Diphtheria = 46),
        interval = "confidence",
        level = 0.95)
```

Confidence intervals for the coefficients:

We are 95% confident that, holding all else constant, for every additional unit of adult mortality (log scale) the average life expectancy will decrease by between 5.457 and 19.783 years.

We are 95% confident that, holding all else constant that for every additional 1% of 1 year olds vaccinated for diphtheria, the average life expectancy will increase by between 8.797 and 9.731 years.

Confidence and Prediction intervals for specific values:

We are 95% confident that the average (over all countries) life expectancy in countries with an average adult mortality of 169 per 1000 and a diphtheria vaccination rate of 72% is between 66.912 years and 67.946 years

We are 95% confident that the life expectancy of a country with an average adult mortality of 169 per 1000 and a diphtheria vaccination rate of 72% will be between 61.650 years and 72.754 years.

We are 95% confident that the average (over all countries) life expectancy in countries with an average adult mortality of 169 per 1000 and a diphtheria vaccination rate of 99% will be between 81.301 years and 84.440 years. 

We are 95% confident that the average life expectancy in a country with an average adult mortality of 169 per 1000 and a diphtheria vaccination rate of 46% will be between 45.737 years and 50.278 years.

## Model Assesment

We are also interested in how well the model fits the data. To do this, we look at metrics such as $R^2$, the RMSE, and.... These metrics are important to check and understand because...

```{r, fig.align='center'}
summary(life_expect_lm)

# MSE
anova <- aov(life_expect_lm)  # get ANOVA components
life_expect_anova <- summary(anova)[[1]]  # save data in a usable form
(mse = life_expect_anova["Residuals", "Sum Sq"] / 
    life_expect_anova["Residuals", "Df"])

# RMSE
(sqrt(mse))

# MAE
life_expect_t$fitted_life.expectancy <- life_expect_lm$fitted.values
sum(abs(life_expect_t$Life.Expectancy - life_expect_t$fitted_life.expectancy) / 
  life_expect_lm$df.residual)
```

It turns out that the model fits the data pretty well. An $R^2$ of 87.5% is 
rather high, and indicates that the variables we chose have done a good job 
describing the model. The RMSE is 375.4. Since we transform the response variable
to a square (the scale is transforms to the scale about 3000 to 7000), the RMSE 
is relatively low comparing to the scale. In other words, the average error 
performed by the model is small. The MAE also have a relative small number
comparing to the scale of the transformed model, which also indicated that the
average absolute difference between the outcome and the prediction is small.

### Summary and Conclusions

Life expectancy is one of the most frequently used metrics to determine health status among countries. It is also a good indicator of the development within countries. We conducted an analysis to find out which variables significantly affected life expectancy in 130 different countries. After fitting a multiple linear regression model, we found that the percent of infants immunized against diphtheria has a significant positive impact on life expectancy (the higher the immunization rate the higher the life expectancy). We also found that adult mortality had a negative impact on life expectancy. Based on our results, we can make a reasonable prediction of the life expectancy of a country based on only those two factors.