Data 605 - HW12

Load Packages:

library(tidyverse)

Load Data:

my_url <- "https://raw.githubusercontent.com/geedoubledee/data605_hw12/main/who.csv"
who_df <- read.csv (my_url)

Intro:

The who.csv dataset contains real-world data from 2008. The variables included follow:

Country: name of the country
LifeExp: average life expectancy for the country in years
InfantSurvival: proportion of those surviving to one year or more
Under5Survival: proportion of those surviving to five years or more
TBFree: proportion of the population without TB.
PropMD: proportion of the population who are MDs
PropRN: proportion of the population who are RNs
PersExp: mean personal expenditures on healthcare in US dollars at average exchange rate
GovtExp: mean government expenditures per capita on healthcare, US dollars at average exchange rate
TotExp: sum of personal and government expenditures.

Question 1:

Provide a scatterplot of \(LifeExp\) ~ \(TotExp\) and run simple linear regression. Do not transform the variables. Provide and interpret the F statistics, \(R^2\), standard error,and p-values only. Discuss whether the assumptions of simple linear regression are met.

ggplot(who_df, aes(x = TotExp, y = LifeExp)) +
    geom_point() +
    labs(title = "LifeExp ~ TotExp",
         x = "TotExp", y = "LifeExp")

who_df_lm <- lm(LifeExp ~ TotExp,
                     data = who_df)
who_df_lm

## 
## Call:
## lm(formula = LifeExp ~ TotExp, data = who_df)
## 
## Coefficients:
## (Intercept)       TotExp  
##   6.475e+01    6.297e-05

\(\hat{y} = 64.75 + 0.00006297x\)

ggplot(who_df, aes(x = TotExp, y = LifeExp)) +
    geom_point() +
    geom_smooth(method = lm, se = FALSE) + 
    labs(title = "LifeExp ~ TotExp",
         x = "TotExp", y = "LifeExp")

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

summary(who_df_lm)

## 
## Call:
## lm(formula = LifeExp ~ TotExp, data = who_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.764  -4.778   3.154   7.116  13.292 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.475e+01  7.535e-01  85.933  < 2e-16 ***
## TotExp      6.297e-05  7.795e-06   8.079 7.71e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.371 on 188 degrees of freedom
## Multiple R-squared:  0.2577, Adjusted R-squared:  0.2537 
## F-statistic: 65.26 on 1 and 188 DF,  p-value: 7.714e-14

The F-statistic is more appropriate when you’re looking at multiple variables, but here it is 65.26. That is large, and it is associated with a p-value of 0.00000000000007714, which is very small. So there is a statistically significant relationship between \(TotExp\) and \(LifeExp\). The Adjusted \(R^2\) value is 0.2537, which indicates 25.37% of the variation in \(LifeExp\) can be explained by \(TotExp\). The standard error, 0.000007795, is roughly 8 times lower than the coefficient, 0.00006297, which is in line with what we want to see from a good model.

However, the assumptions of simple linear regression are not met because:

the relationship appears curvilinear based on the scatterplot

ggplot(who_df_lm, aes(x = .fitted, y = .resid)) +
    geom_point() +
    geom_hline(yintercept = 0) +
    labs(title = "Residual vs. Fitted Values",
         x = "Fitted Values", y= "Residuals")

plotting the residuals vs. the fitted values shows us that the residuals are not uniformly scattered above and below zero; there is an obvious pattern, so errors are not independent, and variance is not equal

qqnorm(resid(who_df_lm))
qqline(resid(who_df_lm))

the qq-plot shows us the residuals are not normally distributed; very few of them fall along the normal line

Question 2:

Raise life expectancy to the \(4.6\) power (i.e., \(LifeExp^{4.6}\)). Raise total expenditures to the \(0.06\) power (nearly a log transform, \(TotExp^{.06}\)). Plot \(LifeExp^{4.6}\) as a function of \(TotExp^{.06}\), and re-run the simple regression model using the transformed variables. Provide and interpret the F statistics, \(R^2\), standard error, and p-values. Which model is “better?”

who_df <- who_df %>%
    mutate(TotExp_to_0.06 = TotExp^0.06,
           LifeExp_to_4.6 = LifeExp^4.6)

ggplot(who_df, aes(x = TotExp_to_0.06, y = LifeExp_to_4.6)) +
    geom_point() +
    labs(title = "LifeExp^4.6 ~ TotExp^0.06",
         x = "TotExp^0.06", y = "LifeExp^4.6")

who_df_lm2 <- lm(LifeExp_to_4.6 ~ TotExp_to_0.06,
                     data = who_df)
who_df_lm2

## 
## Call:
## lm(formula = LifeExp_to_4.6 ~ TotExp_to_0.06, data = who_df)
## 
## Coefficients:
##    (Intercept)  TotExp_to_0.06  
##     -736527909       620060216

\(\hat{y} = -736527909 + 620060216x\)

ggplot(who_df, aes(x = TotExp_to_0.06, y = LifeExp_to_4.6)) +
    geom_point() +
    geom_smooth(method = lm, se = FALSE) + 
    labs(title = "LifeExp^4.6 ~ TotExp^0.06",
         x = "TotExp", y = "LifeExp")

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

summary(who_df_lm2)

## 
## Call:
## lm(formula = LifeExp_to_4.6 ~ TotExp_to_0.06, data = who_df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -308616089  -53978977   13697187   59139231  211951764 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -736527910   46817945  -15.73   <2e-16 ***
## TotExp_to_0.06  620060216   27518940   22.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 90490000 on 188 degrees of freedom
## Multiple R-squared:  0.7298, Adjusted R-squared:  0.7283 
## F-statistic: 507.7 on 1 and 188 DF,  p-value: < 2.2e-16

The F-statistic is more appropriate when you’re looking at multiple variables, but here it is 507.7. That is very large, and it is associated with a p-value of 0.00000000000000022, which is very small. So there is a statistically significant relationship between \(TotExp\) and \(LifeExp\) after the transformations. The Adjusted \(R^2\) value is 0.7283, which indicates 72.83% of the variation in \(LifeExp\) can be explained by \(TotExp\) after the transformations. The standard error, 27518940, is roughly 8 times smaller than the coefficient, 620060216, which is more than sufficient for what we want to see from a good model.

The model using the transformations is better because:

the relationship appears clearly linear from the scatterplot after the transformations

ggplot(who_df_lm2, aes(x = .fitted, y = .resid)) +
    geom_point() +
    geom_hline(yintercept = 0) +
    labs(title = "Residual vs. Fitted Values",
         x = "Fitted Values", y= "Residuals")

plotting the residuals vs. the fitted values shows us that the residuals are uniformly scattered above and below zero; there is no obvious pattern, so errors are independent, and variance is equal

qqnorm(resid(who_df_lm2))
qqline(resid(who_df_lm2))

the qq-plot shows us the residuals are much more normally distributed; we do still see some deviations, especially along the left tail, but also on the right tail.

Question 3:

Using the results from the previous question, forecast life expectancy when \(TotExp^{.06} = 1.5\).

a <- 1.5
forecast1 <- round((-736527909 + 620060216 * a)^(1/4.6), 1)

\(LifeExp = 63.3\)

Then forecast life expectancy when \(TotExp^{.06} = 2.5\).

b <- 2.5
forecast2 <- round((-736527909 + 620060216 * b)^(1/4.6), 1)

\(LifeExp = 86.5\)

Question 4:

Build the following multiple regression model and interpret the F Statistics, \(R^2\), standard error, and p-values. How good is the model?

\(LifeExp = b0 + b1 \times PropMd + b2 \times TotExp +b3 \times PropMD \times TotExp\)

who_df_lm3 <- lm(LifeExp ~ PropMD * TotExp, data = who_df)
who_df_lm3

## 
## Call:
## lm(formula = LifeExp ~ PropMD * TotExp, data = who_df)
## 
## Coefficients:
##   (Intercept)         PropMD         TotExp  PropMD:TotExp  
##     6.277e+01      1.497e+03      7.233e-05     -6.026e-03

\(\hat{y} = 62.77 + 1497 \times PropMD + 0.00007233 \times TotExp - 0.006026 \times PropMd \times TotExp\)

summary(who_df_lm3)

## 
## Call:
## lm(formula = LifeExp ~ PropMD * TotExp, data = who_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.320  -4.132   2.098   6.540  13.074 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    6.277e+01  7.956e-01  78.899  < 2e-16 ***
## PropMD         1.497e+03  2.788e+02   5.371 2.32e-07 ***
## TotExp         7.233e-05  8.982e-06   8.053 9.39e-14 ***
## PropMD:TotExp -6.026e-03  1.472e-03  -4.093 6.35e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.765 on 186 degrees of freedom
## Multiple R-squared:  0.3574, Adjusted R-squared:  0.3471 
## F-statistic: 34.49 on 3 and 186 DF,  p-value: < 2.2e-16

The F-statistic is 34.49, which is large, and it is associated with a p-value of 0.00000000000000022, which is very small. So there is a statistically significant relationship between \(PropMD\), \(TotExp\), their interaction term, and \(LifeExp\). The Adjusted \(R^2\) value is 0.3471, which indicates 34.71% of the variation in \(LifeExp\) can be explained by \(PropMD\), \(TotExp\), and their interaction term. The standard error for \(PropMD\) is roughly 5 times smaller than its coefficient, the standard error for \(TotExp\) is roughly 8 times smaller than its coefficient, which is in line with what we want to see from a good model. The standard error for their interaction term is only about 4 times smaller than its coefficient though.

par(mfrow=c(2,2))
plot(who_df_lm3)

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

We can see by plotting the residuals vs. the fitted values that the residuals are not uniformly scattered about zero, and there is a clear pattern. We can see from the qq-plot that many of the residuals don’t follow the normal distribution line. So we don’t have linearity, independent errors, normality of errors, or equal variance, and this model is not particularly good.

Question 5:

Forecast \(LifeExp\) when \(PropMD = .03\) and \(TotExp = 14\). Does this forecast seem realistic? Why or why not?

x <- 0.03
z <- 14
forecast3 <- 62.77 + (1497 * x) + (0.00007233 * z) - (0.006026 * x * z)

\(LifeExp = 107.6784817\)

This doesn’t seem realistic for two reasons:

it’s difficult to predict values beyond your data with any model, and the max LifeExp in our data is 83.
there are natural floors and ceilings for measurements like life expectancy; a value this high would imply there are many people in the population living even beyond this age, and the world as we know it is not there yet.