BANA 7052 Assignment 4 - RFischer

The R code is hidden by default. Please click on Show to view all codes. Thank you.

Assignment 4

Question 1. Alumni Donation Data

Import Data and load libraries from URL: https://bgreenwell.github.io/uc-bana7052/data/alumni.csv

url <- "https://bgreenwell.github.io/uc-bana7052/data/alumni.csv"
alumni <- read.csv(url)
str(alumni)
attach(alumni)
library(tidyverse)
library(investr)
library(car)
library(lmtest)
library(gridExtra)

1) Does the assumption of error normality appear to be violated?

To see if the error normality is violated or not, I am going to use a QQ plot to diagnose it. First I am going to fit the data to a linear regression model.

# fit the model
fit <- lm(alumni_giving_rate ~ percent_of_classes_under_20 + student_faculty_ratio)

# residual diagnostics
qqnorm(residuals(fit), main = "Right-skewed Distribution")
qqline(residuals(fit), col = "blue")

From the QQ plot we can see that the majority of the points are on the straight line while some of the of points are above the line. This tells us that it’s right-skewed.

We can also use the Shapiro-Wilk test to exam the p value:

shapiro.test(residuals(fit))

    Shapiro-Wilk normality test

data:  residuals(fit)
W = 0.94577, p-value = 0.02721

The P value from the Shapiro-Wilk normality test is 0.027, which is less than the threshold value of 0.05, so it indicates the non-normality.

2) Does the assumption of constant error variance appear to be violated?

To diagnose that, I am going to compare the Residuals VS Fitted Values plot. I will be getting the studentized residual and fitted values.

# studentized residual
r.stu <- rstandard(fit)
# fitted value
f.val <- fitted(fit)
plot(f.val, y = r.stu, xlab = "Fitted Value", ylab = "Studentized Residual")
abline(h = c(-2, 0,2), lty = 2, col = 2)

From the Residual vs Fitted value plot, we can see that most points are scattered around the 0 reference line and within the -2 to +2 range. It suggests the assumtpion of constant variance may hold.

We can also use Breush-Pagan test to see if there is non-constant variance.

bptest(fit)

    studentized Breusch-Pagan test

data:  fit
BP = 0.2723, df = 2, p-value = 0.8727

We can see that the P value is 0.8727 which is larger than the 0.05 threshold, this suggests that there is no strong evidence of heteroscedasticity, so we can say that the constant variance is not violated.

3) Does there appear to be any Y or X outliers or influential points?

To detect if there are outliers or influential points, I am going to use the Cook’s distance and leverage plots.

h <- hatvalues(fit)
out <- which(h > 2 * mean(h))
plot(h, type = "h", ylim = extendrange(h, f = 0.15), main = "High Leverage Points")
abline(h = 2 * mean(h), col = 2)

text(out, y = h[out], labels = out, pos = 3, col = 3)

influencePlot(fit)

We can see from the High Leverage Points plot that there are three values that are greater than \(2\bar{h}\), which can be considered as outliers with respect to the predictor values. Also from the influence measure, We can see there are two more values that can be considered influential and maybe outliers.

4) Is there any concern about multicollinearity?

To access the multicollinearity, I will use the vaiance inflation factors (VIFs) to detect. VIF quantifies how much the variance of a regression coefficient is inflated due to multicollinearity with other predictors.

car::vif(fit)

percent_of_classes_under_20       student_faculty_ratio 
                   2.611671                    2.611671 

From the result we can see that the VIF values for both predictors are around 2.61, which is not large enough (normal 10) to be a major concern. Therefore, we can say that there is not much of a concern about multicollinearity.

e) What is the predicted alumni giving rate for an observation with (𝑋1 = 40, 𝑋2 = 5)? Is there any concern about this prediction? Please explain.

To find the predicted alumni giving rate, we are going to first estimate the model.

MLR model with normal errors is given – \(Y_i = \beta_0 + \sum_{j=1}^{p-1} \beta_j X_{ij} + \epsilon_i, i = 1,2, ..., n\).

coef(fit)

                (Intercept) percent_of_classes_under_20       student_faculty_ratio 
                 39.6555835                   0.1661686                  -1.7021103 

From the coefficients result above, we can estimate that the final model is \(\hat{Y} = 39.66 + 0.166 X_1 - 1.702 X_2\).

To find the predicted alumni giving rate:

X1 <- 40
X2 <- 5
Y <- 39.6555835 + 0.1661686 * X1 - 1.7021103 * X2
print(paste("The predicted alumni giving rate is",Y))

[1] "The predicted alumni giving rate is 37.791776"

boxplot(alumni$percent_of_classes_under_20, main = "Boxplot of percent of classes under 20", xlab = "percent of calsses under 20", ylab = "percentage")

boxplot(alumni$student_faculty_ratio, main = "Boxplot of student faculty ratio", xlab = "student faculty ratio", ylab = "ratio")

To determine if there is a concern about the prediction, I am going to see it from two different aspects.

1. From the boxplot we can see the mean and the 1st and 3rd quartiles of each predictors, which we can see that X1 = 40 and X2 equals to 5 is not within that range. If those are outliers of the observed dataset, then the prediction may not be most reliable.

2. We have diagnosed that the residuals do not seem to have normality. The model assumptions are linearity, normality of residuals and constant variance. Since we have a violation against normality, the prediction may not be reliable.

Question 2, Simulation Study. Assume mean function \(E(Y|X = 10 + 5X - 2X^2)\).

a) Generate data with \(X_1 \sim N(\mu = 3, \sigma = 0.5)\), sample size \(n = 100\), and error term \(\epsilon \sim N(\mu = 0, \sigma = 0.5)\)

set.seed(7052)
x1 <- rnorm(100, mean = 3, sd = 0.5)
error <- rnorm(100, mean = 0, sd = 0.5)
y <- 10 + 5 * x1 - 2 * x1^2 + error
data <- data.frame(x1, y)
head(mdata, n=3)

b) Fit a SLR using just X. What is the estimated regression equation? Please conduct model estimation, inference, and residual diagnostics. What do you conclude?

To fit the simple linear regression model \(Y = \beta_0 + \beta_1 X +\epsilon\), I will first use the lm() function to fit the model and run a summary statistic on it to get the coefficient.

model.fit <- lm(y ~ x1, data = data)
summary(model.fit)

Call:
lm(formula = y ~ x1, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2342 -0.2929  0.1358  0.5694  1.9273 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  27.5188     0.4870   56.51   <2e-16 ***
x1           -6.9848     0.1588  -43.99   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8617 on 98 degrees of freedom
Multiple R-squared:  0.9518,    Adjusted R-squared:  0.9513 
F-statistic:  1935 on 1 and 98 DF,  p-value: < 2.2e-16

The intercept \(\beta_0 = 27.52\) and slope value is \(\beta_1 = -6.98\), so we can say the estimated model is \(\hat{y} = 27.52 -6.98x_1\).

From the summary statistics, we can also see that the p-value for both intercept and x1 are very small, both are less than \(2e^-16\), which is less than the 0.05 threshold, so it is statistically significant, there is a negative linear relationship between the response and the predictor.

par(mfrow = c(2,2))
plot(model.fit)

From the residual diagnostic plots, we can see that:

1. There is a violation to the normality. From the QQ plot, we can see points below the trend line, which indicats a left-skewed distribution.

2. There is a violation of constant variance. While the fitted value between 5-10 has consistent small residuals, there are also quite a lot points that below the reference line and drop out of the standardized residual threshold \(r_i^{(stan)} < 3\).

3. From the Residuals vs Leverage plot we can see that there are presence of outliers and influential points. Such as 14, 85, 68 labed in the plot.

c) Update the model from part b by adding a quadratic term. Conduct model estimation, inference, and residual diagnostics. What do you conclude? Does this model seem to fit the data better?

To fit the quadratic model, we can say \(\hat{y} = \beta_0 + \beta_1 x_1 + \beta_2 x^2 + \epsilon\). I am going to assign x_squared to be \(x_1^2\).

data2 <- mutate(data, x_squared = x1^2)
head(data2, n=3)

fit2 <- lm(y ~ x1 + x_squared, data = data2)
summary(fit2)

Call:
lm(formula = y ~ x1 + x_squared, data = data2)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.19046 -0.30582  0.03546  0.31477  1.38463 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.7438     1.0707  10.034  < 2e-16 ***
x1            4.4770     0.7154   6.258 1.06e-08 ***
x_squared    -1.8949     0.1175 -16.131  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4513 on 97 degrees of freedom
Multiple R-squared:  0.9869,    Adjusted R-squared:  0.9866 
F-statistic:  3656 on 2 and 97 DF,  p-value: < 2.2e-16

After refitting the linear model and run the summary statistics, we can estimate based on the \(\beta_0, \beta_1, and \beta_2\) value. \(\beta_0 = 10.74, \beta_1 = 4.48, \beta_2 = -1.89\), therefore, the estimated model is \(\hat{y} = 10.74 + 4.48x_1 - 1.89x_1^2\).

The p value of three coefficients are very small, indicates statistical significance. We can also look at the \(R^2\) value, the \(R^2\) value from the single linear regression model is 0.95 and this model’s \(R^2\) value is increased to 0.99, which increased quite a bit. When \(R^2\) value increase, it indicates a better fit of the model.

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

From the residual diagnostic plots, we can also see that:

1. The QQ plot is behaving better, most of the points are right on the line, indicates a normal distribution, there is one value slightly outside of the line, but there is not big influence for that one.

2. The fitted value vs standardized residual also seems to be better, most of the values are within the -1.5 - +1.5 range, which indicates a constant variance.

3. For cook’s distance and leverages, we can see there are still a couple of outliers but it’s also behaving better than the previous model. Most of the points are on the lefts ide which has a lower score. We could do something to the outliers to make the data more accurate, but these outliers may not be as significant.

From all the reasons above, the 2nd model fit is a better fit of the data.

d) What is the variance inflation factor (VIF) for the quadratic model? Any concern of multicollinearity?

car::vif(fit2)

       x1 x_squared 
 73.97863  73.97863 

By using the vif() function, we can see that the VIF value for both \(x_1\) and \(x_1^2\) are very high, the values are almost 74, which is way pass the normal value of 10. It strongly suggests that there are multicollinearity, which could affect codfficient stability.

e) Now enter the X variable and compare the VIF from d. What did you find? Which VIF is smaller? Please brifly explain the reason.

#method 1:
fit3 <- lm(y ~ x1 + I(x1^2), data = data2)
fit_centered1 <- lm(y ~ I(x1 - mean(x1)) + I((x1 - mean(x1))^2), data = data2)
# Method 2: center x in the training data
data2$x1_centered <- data2$x1 - mean(data2$x1)
fit_centered2 <- lm(y ~ x1_centered + I(x1_centered^2), data = data2)

car::vif(fit_centered1)

    I(x1 - mean(x1)) I((x1 - mean(x1))^2) 
            1.000295             1.000295 

car::vif(fit_centered2)

     x1_centered I(x1_centered^2) 
        1.000295         1.000295 

After center the X variable, the VIF scores dropped tremendously. I used both methods to calculate the VIF and the results are almost equal to 1, which is small enough to say that there is now no multicollinearity, becuase the quadratic term is less correlated with the centered X.

Question 3.

a) What are the assumptions of the normal linear regression model?

The typical assumptions are:

1. Independence: The observation are independent of each other.

2. Constant variance: the variance of the error term is constant across all levels of the predictors. Therefore the residuals should be roughly the same at all fitted values.

3. Normally distributed error. The error term are normally distributed, particularly for inference purposes. When we do t test and obtain confidence intervals, we normally assume that the error is normally distributed.

4. We assume that we have the correct model or at least a useful one.

5. Other assumptions like linearity, meaning the response variable is linearly related to the predictors; and no multicollinearity, meaning predictor variables are not highly correlated with each other are also some common assumptions.

b) Are the residuals from linear regression uncorrelated in general? Please explain. What is the distribution of residuals?

Yes, in a correctly specified linear regression model, the residuals are assumed to be uncorrelated with each other and the assumption of independence holds. When we assume the observations are independent, we imply that the errors should be uncorrelated across observations. The residuals don’t necessarily be normally distributed, but if the error terms are normally distributed, then the residuals will approximate a normal distribution. Residuals should ideally have a mean of zero and unit variance.

c) Why do we plot the residuals against the fitted value rather than against the observed response value?

1. We can spot the pattern in a residuals vs fitted value plot, it can indicate non-linearity or non-constant variance in the data. If residuals are plotted against observed values, it’s harder to detect htese issues related to model fit.

2. Fitted values are the best estimates provided by the model for y. Checking how the residuals behave relative to these values helps assess whether the model’s linear fit is appropriate.

3. The plot can also help spot outliers. The points are deviate significantly from the test of the residuals at high or low fitted values may have a large impact on the model.

d) What is multicollinearity? Computationally, what kind of problems can multicollinearity cause when using ordinary lease squares?

Multicollinearity happens when predictor variables are correlated among themselves. When the correlation is high, there is a significant multicollinearity. Although, multicollinearity does not always prevent us from obtaining a good fit, it could cause some issues:

1. It can cause some of the estimated coefficients to become unstable, like high standard errors. The estimates may become highly sensitive to small changes in the data. It could also increase the variance of the estimated regression coefficients which can make them unstable and sensitive to minor changes in the model.

2. It can complicate the interpretation of the estimated coefficients, because a change in one predictor will almost always correspond with a change in another. It becomes hard to detect the true effect of predictors because their p value become less reliable.