Let’s find an equation to predict the median house value of houses in a neighborhood (medv) based on the percentage of households with low socio-economic status in the neighborhood (lstat)

We are going to use the data in the Boston data set (download it from Canvas). Read the Boston.csv file in R. Call it Boston_df.

Boston_df = read.csv("Boston.csv", header = T)
Boston_df
# Take a look at the Boston_df
str(Boston_df)
## 'data.frame':    506 obs. of  13 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : int  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

Let’s define the random variables that we are focusing on:

Y (dependent variable): medv (median house value in the neighborhood) (IT IS GIVEN IN UNITS OF $1000)

X (independent variable): lstat (percent of households with low socioeconomic status in the neighborhood)

Check if the linear relationship between these two variables is reasonable

  1. Obtain the sample correlation coefficient between X and Y
cor (Boston_df$lstat, Boston_df$medv)
## [1] -0.7376627
cor (Boston_df$medv, Boston_df$lstat)
## [1] -0.7376627
  1. Do a scatter plot between X and Y
plot(Boston_df$lstat, Boston_df$medv, xlab = "Socioeco Status", ylab = "Value")

# If you want to include the least square line in the plot, add the following:

abline(lm(medv ~ lstat, data = Boston_df), col = "red")

Interpretation of r and scatter plot_

The value of r suggests a fairly strong negative linear relationship between medv and lstat. We can confirm in the scatter plot that the relationship is indeed negative.

The scatter plot shows us that the linear relationship is acceptable (it fits the data points quite well), especially, in the range from 5< lstat< 30. However, towards the two ends, the linear equation does not fit the data very well.

A closer look at the graph shows that the data points exhibit some non-linearity,i.e., the data points follow non-linear trend. Then, maybe an exponential relationship or a second degree polynomial will do a better job than a linear equation.

In conclusion: It is reasonable to continue with the linear regression analysis. Even if a linear equation might not be the BEST, it might give us good results.

Apply Simple Linear Regression

We apply linear regression in R by using the lm() function, where lm stands for linear model. It is customary to save the result of applying lm() in a user-defined object

medv_simple_out = lm(medv ~ lstat, data=Boston_df)

# Then, we call the summary() function on the defined object

summary(medv_simple_out)
## 
## Call:
## lm(formula = medv ~ lstat, data = Boston_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.168  -3.990  -1.318   2.034  24.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.55384    0.56263   61.41   <2e-16 ***
## lstat       -0.95005    0.03873  -24.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.216 on 504 degrees of freedom
## Multiple R-squared:  0.5441, Adjusted R-squared:  0.5432 
## F-statistic: 601.6 on 1 and 504 DF,  p-value: < 2.2e-16

The equation we got is the following:

Predicted medv of the neighborhood = 34.55 - 0.95*(lstat of the neighborhood)

Predicted medv = 34.55 - 0.95*(lstat)

Is this equation good?

Part 1: Statistical significance

The t test for B1:

Ho: The slope of the estimated equation is NOT statistically different from zero; thus, there is no linear relationship between lstat and medv.

Ha: The slope of the estimated equation is statistically different from zero; thus, there is a linear relationship between lstat and medv.

The PV that results from the test for the b1 coefficient is very small (PV < 2e-16). This PV is less than alpha. Therefore, we can reject Ho and support Ha. The data give us evidence to conclude that there is a statistically significant relationship between lstat and medv.

t.test(Boston_df$medv, Boston_df$lstat, alternative = c("greater"), var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  Boston_df$medv and Boston_df$lstat
## t = 19.086, df = 1010, p-value < 2.2e-16
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  9.027525      Inf
## sample estimates:
## mean of x mean of y 
##  22.53281  12.65306

OPTIONAL !!!: You could also analyze the ANOVA table and the PV related to the F statistic. To get the ANOVA table, you would use the command: anova (medv_simple_out). You would get the same conclusion we got by using the t test for B1.

anova_test = aov(medv_simple_out)
summary(anova_test)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## lstat         1  23244   23244   601.6 <2e-16 ***
## Residuals   504  19472      39                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpret the slope of the equation

Note: It only makes sense to interpret the slope after you have already showed statistical significance.

Interpretation of the slope: When the percent of households with low socioeconomic status in a neighborhood increases by one percent, we expect the median household value to decrease on average by 950 dollars (0.95* 1000= 950)

Another way of saying the same thing: When we compared two neighborhoods that differ by 1% with respect to the percent of households with low socioeconomic status, we should expect a difference between them in median household value of around 950 dollars.

Is this equation good?

Part 2: Practical significance

R squared is 0.5441 (54.4%)

Interpretation: If we use the linear equation instead of the sample mean to make predictions of Y, we can reduce 54.4 % of the variability in the values of Y.

Since we reduce a little more than 50% of the variability compared to using the sample mean, using the equation is better than using the sample mean to make predictions.

Note: In Statistics, the term variability (or error) refers to the difference between the predicted values of Y and real values of Y.

The RSE is 6.216 (= $ 6,216)

On average, the values of medv deviate from the regression equation by $6,216.

Stated in another way, we can expect the predictions of medv based on the equation to be off by around $6,216.

To know whether the value of RSE is small or not, let’s compare RSE with the values of medv. Are we going to compare RSE with all the values of medv?

No, that’ unfeasible. A good strategy is to compare RSE with the average of medv. That’s precisely what the coefficient of variation does. Let’s compute the coefficient of variation:

coefficient of variation = RSE/sample mean

# As a proportion
6.216/mean(Boston_df$medv)
## [1] 0.2758644
# As a percentage
(6.216/mean(Boston_df$medv))*100
## [1] 27.58644

Converted into a %, the coefficient of variation is 27.6%. An RSE of $6,216 does not seem to be as low as we would desire IN THIS PROBLEM. It is more desirable to get a coefficient of variation close to 10% (we want the magnitude of the error to be around 10% compared to a typical value of Y). We should attempt to get a smaller RSE.

Although there is NOT a universally valid criterion to know when the coefficient of variation is low enough, in this class, we will consider a coefficient of variation <= 20% as a desirable outcome ((we want the magnitude of the error to be at most 20% compared to a typical value of Y). A value <= 10% would be a very desirable outcome.

In this problem, an RSE of $ 6,216 does NOT seem to be as low as we would desire. We should attempt to get a smaller RSE (How? Using a different regression equation).

Overall conclusion AFTER assessing the statistical and practical significance:

The equation is statistically significant. But, is it practically significant?

It is better than the sample mean to make predictions of medv (this statement is based on the R sq value that we got). However, the RSE is not as good as we would like (we need to try to get a lower RSE by using multiple regression or non-linear regression).

In any case, for the purpose of this class session, we are going to assume that we made the decision to use the equation. How are going to use the equation? Next…

Using the equation to make predictions of medv

Let’s use the equation to predict medv for three new neighborhoods with lstat of 4.5%, 5%, and 7%.

Let’s use the predict() function to do this.

predict(medv_simple_out, data.frame(lstat=c(4.5, 5, 7)))
##        1        2        3 
## 30.27862 29.80359 27.90350

These are the predicted values of medv GIVEN by the equation for these three new neighborhoods.

OPTIONAL: What if you plug in the values of lstat directly in the equation?

34.55 - 0.95*(4.5)
## [1] 30.275
34.55 - 0.95*(5)
## [1] 29.8
34.55 - 0.95*(7)
## [1] 27.9

Analysis of residuals to check the validity assumptions

Residuals are the errors of the equation.

Residuals= real values of Y - predicted values of Y

You can get the values of the residuals by using the function residuals()

# Let's use head() to only print the value of the first six residuals (just to keep the output short. You do not want to print 506 residuals )

head(residuals(medv_simple_out))
##          1          2          3          4          5          6 
## -5.8225951 -4.2703898  3.9748580  1.6393042  6.7099222 -0.9040837

You can also get the residuals by subtracting the predicted values from the actual values. Next:

head(Boston_df$medv - predict(medv_simple_out))
##          1          2          3          4          5          6 
## -5.8225951 -4.2703898  3.9748580  1.6393042  6.7099222 -0.9040837

To check assumptions 1 (linearity) and 4 (all residuals have equal variance) we do a graph to plot the residuals versus the predicted values of Y.

plot(predict(medv_simple_out), residuals(medv_simple_out), xlab = "Predict", ylab="Residuals")

# abline() allows us to draw a horizontal line at y=0 (because 0 is the mean of the residuals)

abline(h=0, col="red")

Interpretation: In this case, we can clearly see that the residuals show a nonlinear pattern. The data points seem to follow a parabolic shape; therefore, a quadratic model might be more appropriate than a linear model ( assumption 1 is not satisfied). We had already observed this nonlinear pattern when we did the original scatter plot. However, this graph shows the non-linear pattern more clearly.

When assumption 1 is not satisfied, it is really hard to evaluate assumption 4. So, NO need to assess assumption 4 in this case.

To check assumption 3 (whether the residuals follow a normal distribution), we are going to conduct a hypothesis test: the Shapiro test to check for normality

Ho: The residuals follow a Normal distribution Ha: The residuals do NOT follow a Normal distribution

shapiro.test(residuals(medv_simple_out))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(medv_simple_out)
## W = 0.87857, p-value < 2.2e-16

PV= 2.2e-16, which is a lot smaller than alpha (0.05). Thus, we reject Ho and support Ha. The data is giving us evidence that the residuals do NOT follow a Normal distribution. Assumption 3 is not satisfied either!

This usually happens. When assumption 1 (the linearity assumption) is not satisfied and a non-linear model is obviously better than a linear model, assumption 3 tends to be compromised too.