Part I: Gauss-Markov Assumptions

The assumptions provided below are listed in the same sequence from matrix_OLS_NYU_notes.pdf and use info from statisticshowto.com. These are not necessarily definitions, rather they are descriptions

Linearity

Non-Technical: Changes in input variables have a straightforward and consistent impact on the outcome variable. It’s like assuming a straight line on a graph.

Technical: Each parameter in the model must have a linear relationship with the dependent variable

Acceptable models:

\[ Y=\beta_0+\beta_1X_1+\beta_2X_2+u_i \]

\[ Y=\beta_0+\beta_1(X_1)^2+\beta_2X_2+u_i \]

\[ Y=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3(X_1X_2)+u_i \]

Unacceptable models:

\[ Y=\beta_0+(\beta_1)^2X_1+\beta_2X_2+u_i \]

\[ Y=\beta_0+\beta_1X_1\beta_2X_2+u_i \]
Full Rank

Non-Technical: All variables provide unique and valuable information that does not overlap with other variables.

Technical: The design matrix, represented as X has the maximum possible rank, meaning all its columns are linearly independent. No to columns are the exact same in X.
Zero Conditional Mean

Non-Technical: On average, our predictions are neither consistently too high nor too low compared to the actual values. We want our model to be unbiased.

Technical: The expected value of the error term conditional on the independent variables is zero, indicating that the errors are not systematically over- or under-predicting the dependent variable.

\[ E[u|X]=0 \]
Homoskedastic

Non-Technical: The model’s predictions are equally scattered across the range of possible outcomes, without any noticeable pattern or change in how spread out they are. It’s like saying the model’s guesses are consistently “even” in their uncertainty.

Technical: The error term is constant across all levels of the independent variables, indicating consistent scatter of residuals.

\[ E(\epsilon\epsilon'|X) = \sigma^2I \]
No Autocorrelation

Non-Technical: We assume that each piece of information we have is like a separate, independent piece of the puzzle. There’s no connection between them; they don’t rely on each other. It’s like each new bit of data is a fresh start with no memory of the past.

Technical: The error terms are not correlated with each other across different observations or time periods. It is backed by the same mathematical assumption for homoskedasticity

\[ E(\epsilon\epsilon'|X) = \sigma^2I \]
Exogeneity

Non-Technical: We assume that the things we’re using to make predictions aren’t affected by our predictions themselves. It’s like saying our crystal ball doesn’t change the events it’s trying to predict; they’re independent of each other.

Technical: The independent variables are not influenced by the error term, ensuring that they are determined by external factors and not by the model itself.
Central Limit Theorem

Non-Technical: When we gather a lot of data and calculate averages, those averages tend to become more reliable and predictable, making it easier to understand and work with data in many practical situations.

Technical: As the sample size approaches infinity, the distribution of the sample mean approaches a normal distribution, regardless of the underlying population distribution.

\[ \epsilon|X ∼N[0,\sigma^2I] \]

Part II: Regression

Imports (skip)

library(AER)

## Loading required package: car

## Loading required package: carData

## Loading required package: lmtest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

## Loading required package: survival

library(knitr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

data(MASchools)
MASchools <- MASchools %>%
  select(-district)
MASchools <- MASchools %>%
  mutate(special = special / 100,
         lunch = lunch / 100,
         english = english / 100)
row.names(MASchools)<-MASchools$municipality
MASchools<-MASchools[, -1]
library(sjPlot)
library(sjmisc)
library(sjlabelled)

## 
## Attaching package: 'sjlabelled'

## The following object is masked from 'package:dplyr':
## 
##     as_label

library(ggplot2)

## 
## Attaching package: 'ggplot2'

## The following object is masked from 'package:sjlabelled':
## 
##     as_label

Data Description

Dataframe: MASchools

Format: cross-sectional data frame containing 220 observations on 16 variables

Background: MASchools contains data on test performance, school characteristics and student demographic backgrounds for school districts in Massachusetts.

Additional details: The Massachusetts data are district-wide averages for public elementary school districts in 1998. The test score is taken from the Massachusetts Comprehensive Assessment System (MCAS) test, administered to all fourth graders in Massachusetts public schools in the spring of 1998. The test is sponsored by the Massachusetts Department of Education and is mandatory for all public schools. The data analyzed here are the overall total score, which is the sum of the scores on the English, Math, and Science portions of the test. Data on the student-teacher ratio, the percent of students receiving a subsidized lunch and on the percent of students still learning english are averages for each elementary school district for the 1997--1998 school year and were obtained from the Massachusetts department of education. Data on average district income are from the 1990 US Census.

Link: MASchools: Massachusetts Test Score Data

Data at a glance

kable(head(MASchools))

	expreg	expspecial	expbil	exptot	scratio	special	lunch	stratio	income	score4	score8	salary	english
Abington	4201	7375.69	0	4646	16.6	0.146	0.118	19.0	16.379	714	691	34.3600	0.0000000
Acton	4129	8573.99	0	4930	5.7	0.174	0.025	22.6	25.792	731	NA	38.0630	0.0124611
Acushnet	3627	8081.72	0	4281	7.5	0.121	0.141	19.3	14.040	704	693	32.4910	0.0000000
Agawam	4015	8181.37	0	4826	8.6	0.211	0.121	17.9	16.111	704	691	33.1060	0.0032258
Amesbury	4273	7037.22	0	4824	6.1	0.168	0.174	17.5	15.423	701	699	34.4365	0.0000000
Amherst	5183	10595.80	6235	6454	7.7	0.172	0.268	15.7	11.144	714	NA	NA	0.0392157

str(MASchools)

## 'data.frame':    220 obs. of  14 variables:
##  $ expreg    : int  4201 4129 3627 4015 4273 5183 4685 5518 5009 3823 ...
##  $ expspecial: num  7376 8574 8082 8181 7037 ...
##  $ expbil    : int  0 0 0 0 0 6235 0 0 0 12943 ...
##  $ expocc    : int  0 0 0 0 0 0 0 0 0 11519 ...
##  $ exptot    : int  4646 4930 4281 4826 4824 6454 5537 6405 5649 4814 ...
##  $ scratio   : num  16.6 5.7 7.5 8.6 6.1 ...
##  $ special   : num  0.146 0.174 0.121 0.211 0.168 ...
##  $ lunch     : num  0.118 0.025 0.141 0.121 0.174 ...
##  $ stratio   : num  19 22.6 19.3 17.9 17.5 ...
##  $ income    : num  16.4 25.8 14 16.1 15.4 ...
##  $ score4    : num  714 731 704 704 701 714 725 717 702 701 ...
##  $ score8    : num  691 NA 693 691 699 NA 728 715 705 688 ...
##  $ salary    : num  34.4 38.1 32.5 33.1 34.4 ...
##  $ english   : num  0 0.01246 0 0.00323 0 ...

summary(MASchools, na.count = TRUE)

##      expreg       expspecial        expbil           expocc          exptot    
##  Min.   :2905   Min.   : 3832   Min.   :     0   Min.   :    0   Min.   :3465  
##  1st Qu.:4065   1st Qu.: 7442   1st Qu.:     0   1st Qu.:    0   1st Qu.:4730  
##  Median :4488   Median : 8354   Median :     0   Median :    0   Median :5155  
##  Mean   :4605   Mean   : 8901   Mean   :  3037   Mean   : 1104   Mean   :5370  
##  3rd Qu.:4972   3rd Qu.: 9722   3rd Qu.:     0   3rd Qu.:    0   3rd Qu.:5789  
##  Max.   :8759   Max.   :53569   Max.   :295140   Max.   :15088   Max.   :9868  
##                                                                                
##     scratio          special           lunch           stratio     
##  Min.   : 2.300   Min.   :0.0810   Min.   :0.0040   Min.   :11.40  
##  1st Qu.: 6.100   1st Qu.:0.1338   1st Qu.:0.0530   1st Qu.:15.80  
##  Median : 7.800   Median :0.1545   Median :0.1055   Median :17.10  
##  Mean   : 8.107   Mean   :0.1597   Mean   :0.1532   Mean   :17.34  
##  3rd Qu.: 9.800   3rd Qu.:0.1792   3rd Qu.:0.2003   3rd Qu.:19.02  
##  Max.   :18.400   Max.   :0.3430   Max.   :0.7620   Max.   :27.00  
##  NA's   :9                                                         
##      income           score4          score8          salary     
##  Min.   : 9.686   Min.   :658.0   Min.   :641.0   Min.   :24.96  
##  1st Qu.:15.223   1st Qu.:701.0   1st Qu.:685.0   1st Qu.:33.80  
##  Median :17.128   Median :711.0   Median :698.0   Median :35.88  
##  Mean   :18.747   Mean   :709.8   Mean   :698.4   Mean   :35.99  
##  3rd Qu.:20.376   3rd Qu.:720.0   3rd Qu.:712.0   3rd Qu.:37.96  
##  Max.   :46.855   Max.   :740.0   Max.   :747.0   Max.   :44.49  
##                                   NA's   :40      NA's   :25     
##     english        
##  Min.   :0.000000  
##  1st Qu.:0.000000  
##  Median :0.000000  
##  Mean   :0.011177  
##  3rd Qu.:0.008859  
##  Max.   :0.244939  
##

Variables

municipality: character. Municipality name.
expreg: Expenditures per pupil, regular.
expspecial: Expenditures per pupil, special needs.
expbil: Expenditures per pupil, bilingual.
expocc: Expenditures per pupil, occupational.
exptot: Expenditures per pupil, total.
scratio: Students per computer.
special: Special education students (per cent).
lunch: Percent qualifying for reduced-price lunch.
stratio: Student-teacher ratio.
income: Per capita income.
score4: 4th grade score (math + English + science).
score8: 8th grade score (math + English + science).
salary: Average teacher salary.
english: Percent of English learners.

Regression model

\[ score8=\beta_0+\beta_1(log(salary))+\beta_2(stratio)+\beta_3(log(income))+\beta_4(special)+\beta_5(log(exptot))+\beta_6(english)+u_i \]

For the sake of simplicity, we’ll focus on grade 8 students. The primary reason being that they are older and more mentally developed and have been through more of the education system. Each of the above variables seemed like good external factors influencing test scores. The log of salary, income, and exptot were taken to examine the impact a marginal change in the variables would have on test scores.

Hypotheses

salary: Positive

A higher salary should mean teachers are held to a higher standard or do higher quality work
stratio: Negative

As the class size gets larger, teachers have less time to dedicate to individual students
income: Positive

Higher income communities may be able to pay for learning opportunities outside of the classroom
special: Negative

Students who have special needs may have a more difficult time performing on exams. Thus a high concentration may bring down overall test scores
exptot: Positive

More expenses per student means that more resources are provided for students to learn from

Model

my_reg <- lm(score8 ~ log(salary) + stratio + log(income) + special + log(exptot), data=MASchools)
tab_model(my_reg)

	score8
Predictors	Estimates	CI	p
(Intercept)	858.06	751.79 – 964.33	<0.001
salary [log]	-9.31	-33.41 – 14.79	0.447
stratio	-1.98	-2.81 – -1.16	<0.001
income [log]	72.57	64.49 – 80.64	<0.001
special	-78.57	-128.49 – -28.64	0.002
exptot [log]	-33.65	-46.23 – -21.06	<0.001
Observations	163
R² / R² adjusted	0.779 / 0.772

Part III: Plots and Analysis

Auto plots

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

Manual Plots

Residual vs Fitted

resdf <- data.frame(fitted   = my_reg$fitted.values, 
                 residual = my_reg$residuals)

# Create a residual plot using ggplot2
ggplot(data = resdf, 
       mapping = aes(x = fitted, 
                     y = residual)
       ) +
  geom_point() +
  geom_smooth(method = "loess") + 
  xlab("Fitted Values") +
  ylab("Residuals") +
  ggtitle("Residual Plot")

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

The curve pattern of the residual plot suggests that the model may be nonlinear. This may be fixed by omitting variables that may skew the model or suggest an entire change in the model is necessary

Normal Q-Q

qqdf <- data.frame(residuals = residuals(my_reg))
ggplot(data = qqdf, 
       mapping = aes(sample = residuals)
       ) +
  stat_qq() +
  stat_qq_line() +
  xlab("Theoretical Quantiles") +
  ylab("Sample Quantiles") +
  ggtitle("Normal QQ Plot of Residuals")

The data on the plot keeps a linear shape that remains mostly even with the theoretical line. This suggests that the data is indeed normally distributed. One cause for concern is a variable around theoretical quantile 2.75. It appears that it could possibly be an outlier given that it is significantly further from the line than any other variable.

scale_location_data <- data.frame(
  Fitted_Values = fitted(my_reg),
  Sqrt_Abs_Std_Residuals = sqrt(abs(rstandard(my_reg))))

scale_location_plot <- ggplot(data = scale_location_data, 
                              mapping = aes(x = Fitted_Values, 
                                            y = Sqrt_Abs_Std_Residuals)
                              ) +
  geom_point() +
  geom_smooth(method = "lm", 
              se = FALSE, 
              color = "blue") + 
  labs(x = "Fitted Values", 
       y = "Square Root of Absolute Standardized Residuals") +
  ggtitle("Scale-Location (Spread-Location) Plot")

print(scale_location_plot)

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

This plot focuses on identifying heteroskedasticity. There is an ever-so-faint but still noticeable trend of an upward slope. This suggests that there is heteroskedasticity in the model.

Residuals vs Leverage

cooksd <- cooks.distance(my_reg)

residuals_vs_leverage_data <- data.frame(
  Leverage = hatvalues(my_reg),
  Residuals = residuals(my_reg),
  Cooks_Distance = cooksd)

residuals_vs_leverage_plot <- ggplot(residuals_vs_leverage_data, 
                                     aes(x = Leverage, y = Residuals)) +
  geom_point(aes(size = Cooks_Distance), alpha = 0.7) +  
  labs(x = "Leverage", 
       y = "Residuals") +
  ggtitle("Residuals vs. Leverage Plot with Cook's Distance") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +  
  scale_size_continuous(range = c(2, 10), breaks = c(0.05, 0.1, 0.2), labels = c("0.05", "0.1", "0.2")) +  
  theme_minimal()
label_threshold <- 0.05
residuals_vs_leverage_plot +
  geom_text(data = subset(residuals_vs_leverage_data, Cooks_Distance > label_threshold), 
            aes(label = rownames(subset(residuals_vs_leverage_data, Cooks_Distance > label_threshold))),
            position = position_nudge(y = 0.15), col = "red", size = 3)

The residual-vs-leverage plot is especially relevant given previous concerns that there may be an outlier. This plot further backs that suspicion of an outlier with Holyoke and potentially Harvard. The point has significant leverage and seems to have a decent error. Omission of the variable could have a decent impact on cleaning up the model.

Current conclusion

My hypotheses for salary and expenditures are supposedly false given the model. However, the model appears to be flawed given the nonlinearity suggested in the residual-fitted plot. Additionally, the p-value for salary shows the variable may not be significant anyway.

Refining the model

\[ score8=\beta_0+\beta_1(stratio)+\beta_2(log(income))+\beta_3(special)+\beta_4(log(exptot))+u_i \]

After discovering that salary isn’t linked, it was omitted from the model. Both Holyoke and Harvard were also omitted as potential outliers to observe the outcome.

MASchools2 <- subset(MASchools, !(row.names(MASchools) %in% c("Holyoke", "Harvard")))
reg2 <- lm(score8 ~ stratio + log(income) + special + log(exptot), data=MASchools2)
tab_model(reg2)

	score8
Predictors	Estimates	CI	p
(Intercept)	851.19	758.23 – 944.15	<0.001
stratio	-1.99	-2.72 – -1.26	<0.001
income [log]	68.65	62.27 – 75.03	<0.001
special	-61.92	-108.01 – -15.83	0.009
exptot [log]	-35.75	-46.64 – -24.86	<0.001
Observations	178
R² / R² adjusted	0.777 / 0.772

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

resdf2 <- data.frame(fitted   = reg2$fitted.values, 
                 residual = reg2$residuals)

# Create a residual plot using ggplot2
ggplot(data = resdf, 
       mapping = aes(x = fitted, 
                     y = residual)
       ) +
  geom_point() +
  geom_smooth(method = "loess") + 
  xlab("Fitted Values") +
  ylab("Residuals") +
  ggtitle("Residual Plot")

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

Conclusion

This is not an adequate model to determine test scores due to its nonlinearity. Further analysis with tools not yet discussed in class may be required to truly refine the model. Otherwise, it may just need to be completely rebuilt.

Turner_Discussion4

Cam Turner

2023-09-23