HW 12
Daria Dubovskaia
1. Assignment
- 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.
- 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 r 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?”
- Using the results from 3, forecast life expectancy when TotExp^.06 =1.5. Then forecast life expectancy when TotExp^.06=2.5.
- 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
- Forecast LifeExp when PropMD=.03 and TotExp = 14. Does this forecast seem realistic? Why or why not?
2. Load data
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
The file will be downloaded from the Github repository to csv file using read.csv function.
who <- read.csv('https://raw.githubusercontent.com/ex-pr/Dataset/main/who.csv', header=TRUE, sep=",", check.names=FALSE)
summary(who)## Country LifeExp InfantSurvival Under5Survival
## Length:190 Min. :40.00 Min. :0.8350 Min. :0.7310
## Class :character 1st Qu.:61.25 1st Qu.:0.9433 1st Qu.:0.9253
## Mode :character Median :70.00 Median :0.9785 Median :0.9745
## Mean :67.38 Mean :0.9624 Mean :0.9459
## 3rd Qu.:75.00 3rd Qu.:0.9910 3rd Qu.:0.9900
## Max. :83.00 Max. :0.9980 Max. :0.9970
## TBFree PropMD PropRN PersExp
## Min. :0.9870 Min. :0.0000196 Min. :0.0000883 Min. : 3.00
## 1st Qu.:0.9969 1st Qu.:0.0002444 1st Qu.:0.0008455 1st Qu.: 36.25
## Median :0.9992 Median :0.0010474 Median :0.0027584 Median : 199.50
## Mean :0.9980 Mean :0.0017954 Mean :0.0041336 Mean : 742.00
## 3rd Qu.:0.9998 3rd Qu.:0.0024584 3rd Qu.:0.0057164 3rd Qu.: 515.25
## Max. :1.0000 Max. :0.0351290 Max. :0.0708387 Max. :6350.00
## GovtExp TotExp
## Min. : 10.0 Min. : 13
## 1st Qu.: 559.5 1st Qu.: 584
## Median : 5385.0 Median : 5541
## Mean : 40953.5 Mean : 41696
## 3rd Qu.: 25680.2 3rd Qu.: 26331
## Max. :476420.0 Max. :482750glimpse(who)## Rows: 190
## Columns: 10
## $ Country <chr> "Afghanistan", "Albania", "Algeria", "Andorra", "Angola~
## $ LifeExp <int> 42, 71, 71, 82, 41, 73, 75, 69, 82, 80, 64, 74, 75, 63,~
## $ InfantSurvival <dbl> 0.835, 0.985, 0.967, 0.997, 0.846, 0.990, 0.986, 0.979,~
## $ Under5Survival <dbl> 0.743, 0.983, 0.962, 0.996, 0.740, 0.989, 0.983, 0.976,~
## $ TBFree <dbl> 0.99769, 0.99974, 0.99944, 0.99983, 0.99656, 0.99991, 0~
## $ PropMD <dbl> 0.000228841, 0.001143127, 0.001060478, 0.003297297, 0.0~
## $ PropRN <dbl> 0.000572294, 0.004614439, 0.002091362, 0.003500000, 0.0~
## $ PersExp <int> 20, 169, 108, 2589, 36, 503, 484, 88, 3181, 3788, 62, 1~
## $ GovtExp <int> 92, 3128, 5184, 169725, 1620, 12543, 19170, 1856, 18761~
## $ TotExp <int> 112, 3297, 5292, 172314, 1656, 13046, 19654, 1944, 1907~3. Life expectancy vs. Total expenditures
3.1 Scatter plot
ggplot(who, aes(x=TotExp, y=LifeExp)) +
geom_point (color="blue") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))+
labs(x = 'Total expenditures', y = "Life expectancy", title = "Life expectancy vs. Total expenditures")
3.2 Linear regression
Using lm() function, we will build the linear regression model for Life expectancy ~ Total expenditures.
Summary() function will show all the details about the model.
F statistics, 65.26, compares the current model to a model that has only the intercept parameter, it gives the same information as that on the line for the slope estimate. It suggests we can reject a regression model with a zero coefficient.
R^2 is 0.2577, it is a measure of how well the model describes the measured data, the value is pretty small. Our model describes only 25% of the data, so model doesn’t really fit the data.
The Standard error column shows the statistical standard error for each of the coefficients, it should be at least five to ten times smaller than the corresponding coefficient for a good model. The estimated value is 6.297e-05, the std. error is 7.795e-06 which is 8 times smaller 6.297e-05/7.795e-06=8.079 which is almost 10 times smaller which is the sign for a good model.
p-value is 7.714e-14 which is a tiny value which can help us to assume that there is strong evidence of a linear relationship between Life expectancy and Total expenditures.
The low R^2 doesn’t let us to assume that linear model fits our observation, we need additional residual analysis to make a final conclusion.
model <- lm(who$LifeExp ~ who$TotExp)
summary(model)##
## Call:
## lm(formula = who$LifeExp ~ who$TotExp)
##
## 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 ***
## who$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-14The residual analysis doesn’t look right. The q-q plot and residual plot shows that residuals are not normally distributed, linear modeal doesn’t fit the data.
par(mfrow=c(2,2))
plot(model)
4. Linear regression with transformed variables
We will raise life expectancy to the 4.6 power, total expenditures to the 0.06 power and plot the results.
We see that now it looks more like linear relationship.
LifeExp_4.6 <- who$LifeExp^4.6
TotExp_0.06 <- who$TotExp^0.06
ggplot(who, aes(x=TotExp_0.06, y=LifeExp_4.6 )) +
geom_point (color="blue") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))+
labs(x = 'Total expenditures^4.6', y = "Life expectancy^0.06", title = "Life expectancy vs. Total expenditures transformed") +
geom_smooth(method = "lm")
Let’s re-run the linear model for transformed values.
F statistics, 507.7 is much better comparing to the previous model (65.26). R^2 is 0.7298, it is a measure of how well the model describes the measured data, the value is pretty small. Our model describes only 72.98% of the data, the model highly fits the data.
The estimated value is 620060216, the std. error is 27518940 which is 22 times smaller 620060216/27518940=22.53, it is too much which is not the sign for a good model. The previous model had 8 times difference which is much better.
p-value is < 2.2e-16 which is a tiny value which can help us to assume that there is strong evidence of a linear relationship between Life expectancy and Total expenditures.The same we had for the previous model.
This model has a much better value for R^2 but the standard error is too high, this model looks better than the previous one, the scatter plot shows linear relationship between transformed values (we didn’t have it for the previous model) but I would analyze residuals to make the final conclusion.
model_transformed <- lm(LifeExp_4.6 ~ TotExp_0.06)
summary(model_transformed)##
## Call:
## lm(formula = LifeExp_4.6 ~ TotExp_0.06)
##
## 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_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-16The residual analysis looks much better for this transformed model. The residuals normally distributed and reside around the zero, uniformly scattered above and below zero though there is some variance on the left of the q-q plot. This model fits the data much better.
par(mfrow=c(2,2))
plot(model_transformed)
5. Forecast life expectancy
Let’s forecast life expectancy when TotExp^.06 =1.5 and when TotExp^.06=2.5.
We will use model summary (intercept and estimate) to build an equation of our model.
The life expectancy values are 63.311 for TotExp^.06 =1.5 and 86.506 for TotExp^.06=2.5.
TotExp_values = c(1.5, 2.5)
model_summary <- summary(model_transformed)
intercept <-model_summary$coefficients[1,1]
estimate <- model_summary$coefficients[2,1]
for (i in seq_along(TotExp_values)) {
print((intercept + estimate * TotExp_values[[i]])^(1/4.6))
}## [1] 63.31153
## [1] 86.506456. Multiple regression model
The equation for the multiple regression model: LifeExp = b0+b1 x PropMd + b2 x TotExp +b3 x PropMD x TotExp.
F statistics, 34.49, which is good for a 3 regression degrees of freedom and 120 residual degrees of freedom.
R^2 is 0.3574. Our model describes only 35.74% of the data, so model doesn’t really fit the data.
The Standard error is 8 times smaller than the estimated value for TotExp and 5 times for PropMD, which is a good sign.
p-value is < 2.2e-16 which is a tiny value which can help us to assume that there is strong evidence of a linear relationship between Life expectancy and Total expenditures.
The low R^2 doesn’t let us to assume that linear model fits our observation, we need additional residual analysis to make a final conclusion.
multiple_model <- lm(LifeExp ~ PropMD + TotExp + PropMD*TotExp, data=who)
summary(multiple_model)##
## Call:
## lm(formula = LifeExp ~ PropMD + TotExp + PropMD * TotExp, data = who)
##
## 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-16The residuals are not normally distributed, there is a certain pattern on the residual plot. The model doesn’t fit the data.
par(mfrow=c(2,2))
plot(multiple_model)## Warning in sqrt(crit * p * (1 - hh)/hh): ñîçäàíû NaN
## Warning in sqrt(crit * p * (1 - hh)/hh): ñîçäàíû NaN
Let’s try to use LifeExp_4.6 and TotExp_0.06.
F statistics, 180.3, which is much better.
R^2 is 0.7441. Our model describes only 74.41 of the data, the model fits the data in a much better way.
The Standard error is not 5 to 10 times smaller that the estimate value, which is not the sign for a good model.
p-value is < 2.2e-16 which is a tiny value which can help us to assume that there is strong evidence of a linear relationship between Life expectancy and Total expenditures.
This model fits thr data much better but we need additional residual analysis to make a final conclusion.
multiple_model_transformed <- lm(LifeExp_4.6 ~ PropMD + TotExp_0.06 + PropMD*TotExp_0.06, data=who)
summary(multiple_model_transformed)##
## Call:
## lm(formula = LifeExp_4.6 ~ PropMD + TotExp_0.06 + PropMD * TotExp_0.06,
## data = who)
##
## Residuals:
## Min 1Q Median 3Q Max
## -296470018 -47729263 12183210 60285515 212311883
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.244e+08 5.083e+07 -14.253 <2e-16 ***
## PropMD 4.727e+10 2.258e+10 2.094 0.0376 *
## TotExp_0.06 6.048e+08 3.023e+07 20.005 <2e-16 ***
## PropMD:TotExp_0.06 -2.121e+10 1.131e+10 -1.876 0.0622 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 88520000 on 186 degrees of freedom
## Multiple R-squared: 0.7441, Adjusted R-squared: 0.74
## F-statistic: 180.3 on 3 and 186 DF, p-value: < 2.2e-16The residual analysis looks much better for this transformed model. The residuals normally distributed and reside around the zero, uniformly scattered above and below zero though there is some variance on the left of the q-q plot. This model fits the data much better.
par(mfrow=c(2,2))
plot(multiple_model_transformed)
7. Forecast life expectancy
The life expectancy is 107.696 years for PropMD=.03 and TotExp = 14. It doesn’t look realistic, 107 years, when there is only 3% of population are doctors and with $14 for a total expenditures.
PropMD_value = 0.03
TotExp_value = 14
multiple_summary <- summary(multiple_model)
b0 <-multiple_summary$coefficients[1,1]
b1 <- multiple_summary$coefficients[2,1]
b2 <- multiple_summary$coefficients[3,1]
b3 <- multiple_summary$coefficients[4,1]
b0+ b1*PropMD_value + b2*TotExp_value + b3*PropMD_value*TotExp_value## [1] 107.696We can look at the plot below to see that if total expenditures are below $200, the life expectancy is below 65 years. Also, the highest life expectancy in the data frame is 83 years.
max(who$LifeExp)## [1] 83ggplot(who, aes(x=TotExp, y=LifeExp)) +
geom_point (color="blue") +
theme_light() +
theme(plot.title = element_text(hjust = 0.5))+
labs(x = 'Total expenditures', y = "Life expectancy", title = "Life expectancy vs. Total expenditures") +
xlim(0, 200)## Warning: Removed 173 rows containing missing values (geom_point).