data_605_hw12
The attached 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.
Setup
df <- read.csv("who.csv", stringsAsFactors = FALSE)head(df)## Country LifeExp InfantSurvival Under5Survival TBFree PropMD
## 1 Afghanistan 42 0.835 0.743 0.99769 0.000228841
## 2 Albania 71 0.985 0.983 0.99974 0.001143127
## 3 Algeria 71 0.967 0.962 0.99944 0.001060478
## 4 Andorra 82 0.997 0.996 0.99983 0.003297297
## 5 Angola 41 0.846 0.740 0.99656 0.000070400
## 6 Antigua and Barbuda 73 0.990 0.989 0.99991 0.000142857
## PropRN PersExp GovtExp TotExp
## 1 0.000572294 20 92 112
## 2 0.004614439 169 3128 3297
## 3 0.002091362 108 5184 5292
## 4 0.003500000 2589 169725 172314
## 5 0.001146162 36 1620 1656
## 6 0.002773810 503 12543 13046str(df)## 'data.frame': 190 obs. of 10 variables:
## $ Country : chr "Afghanistan" "Albania" "Algeria" "Andorra" ...
## $ LifeExp : int 42 71 71 82 41 73 75 69 82 80 ...
## $ InfantSurvival: num 0.835 0.985 0.967 0.997 0.846 0.99 0.986 0.979 0.995 0.996 ...
## $ Under5Survival: num 0.743 0.983 0.962 0.996 0.74 0.989 0.983 0.976 0.994 0.996 ...
## $ TBFree : num 0.998 1 0.999 1 0.997 ...
## $ PropMD : num 2.29e-04 1.14e-03 1.06e-03 3.30e-03 7.04e-05 ...
## $ PropRN : num 0.000572 0.004614 0.002091 0.0035 0.001146 ...
## $ PersExp : int 20 169 108 2589 36 503 484 88 3181 3788 ...
## $ GovtExp : int 92 3128 5184 169725 1620 12543 19170 1856 187616 189354 ...
## $ TotExp : int 112 3297 5292 172314 1656 13046 19654 1944 190797 193142 ...Q1
- 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 met.
with(df, plot(TotExp, LifeExp, xlab = "Total Expenditures", ylab = "Life Expectancy"))
mod1 <- lm(LifeExp ~ TotExp, data = df)
summary(mod1)##
## Call:
## lm(formula = LifeExp ~ TotExp, data = 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-14From above simple linear regression model summary, we see that the F-statistic is 65.26 with a p-value of 7.714e-14. The F test is actually comparing the model against the null model. The extremely small p-value allows us to reject the null hypothesis, i.e. there is no difference between the simple linear regression model and the null model. However, the adjusted R-squared (0.2537) pointed to a weak relationship between the two variables.
par(mfrow = c(2, 2))
plot(mod1)
In addition, above residual and Q-Q plots clearly indicate that the model does not meet the assumption for linear regression model at all.
Q2
- 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?”
mod2 <- lm((LifeExp^4.6)~I(TotExp^.06), data = df)
with(df, plot(TotExp^0.06, LifeExp^4.6, xlab = "Total Expenditures^0.06", ylab = "Life Expectancy^4.6"))
abline(mod2, col = "red", lty = "dashed")
summary(mod2)##
## Call:
## lm(formula = (LifeExp^4.6) ~ I(TotExp^0.06), data = 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 ***
## I(TotExp^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-16par(mfrow = c(2, 2))
plot(mod2)
Now the adjusted R-squared is 0.7283 and it is considerably better than model 1 that is without transformation. The F-statistic also significantly improved to 507.7 with same degree of freedom as model 1. The extremely small p-value allows us to reject the null hypothesis again. The residual and Q-Q plots suggest constant variance and that the assumptions for regression model are met.
Q3
- Using the results from 3, forecast life expectancy when TotExp^.06 = 1.5. Then forecast life expectancy when TotExp^.06 = 2.5.
lifeExp1 = (mod2$coefficients[[2]] * 1.5 + mod2$coefficients[[1]]) ^(1/4.6)
lifeExp2 = (mod2$coefficients[[2]] * 2.5 + mod2$coefficients[[1]]) ^(1/4.6)
print(glue("When TotExp^.06 = 1.5, the life expectancy is approximately ", round(lifeExp1, 1)))## When TotExp^.06 = 1.5, the life expectancy is approximately 63.3print(glue("When TotExp^.06 = 2.5, the life expectancy is approximately ", round(lifeExp2, 1)))## When TotExp^.06 = 2.5, the life expectancy is approximately 86.5Q4
- 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 x PropMd + b2 x TotExp +b3 x PropMD x TotExp
mod3 <- lm(LifeExp ~ PropMD+TotExp+(PropMD*TotExp), data = df)
summary(mod3)##
## Call:
## lm(formula = LifeExp ~ PropMD + TotExp + (PropMD * TotExp), data = 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-16par(mfrow = c(2, 2))
plot(mod3)
The additional predictors and interaction terms improves the prediction and performs slightly better than model 1. The adjusted R squared and F-statistic are higher. And the extremely small p-value also rejected the null hypothesis. However, the residual and Q-Q plots once again indicate that the model (data without transformation) is not suitable for regression model.
hist(mod3$residuals)
The residuals are heavily skewed to the left.
Q5
- Forecast LifeExp when PropMD=.03 and TotExp = 14. Does this forecast seem realistic? Why or why not?
x = data.frame(PropMD = .03, TotExp = 14)
predLifeExp = predict(mod3, x)[[1]]
print(glue("When PropMD = .03 and TotExp = 14, the life expectancy is approximately ", round(predLifeExp, 1), " which seems very unrealistic. The maximum life expectancy is only ", max(df$LifeExp), " in our dataset."))## When PropMD = .03 and TotExp = 14, the life expectancy is approximately 107.7 which seems very unrealistic. The maximum life expectancy is only 83 in our dataset.