library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bexar_data <- read.csv("bexar_schools.csv")
clean_data <- bexar_data %>%
  select(DISTNAME, DPETECOP, DA0CC21R) %>%
  filter(!is.na(DA0CC21R))
summary_stats <- clean_data %>%
  summarize(
    mean_disadvantaged = mean(DPETECOP),
    sd_disadvantaged = sd(DPETECOP),
    median_disadvantaged = median(DPETECOP),
    mean_college_ready = mean(DA0CC21R, na.rm = TRUE),
    sd_college_ready = sd(DA0CC21R, na.rm = TRUE),
    median_college_ready = median(DA0CC21R, na.rm = TRUE),
    n = n()
  )
print(summary_stats)
##   mean_disadvantaged sd_disadvantaged median_disadvantaged mean_college_ready
## 1           62.82069         27.75069                 68.7           26.67241
##   sd_college_ready median_college_ready  n
## 1         22.49381                 19.4 29
ggplot(clean_data, aes(x = DPETECOP, y = DA0CC21R)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, color = "blue") +
  labs(
    title = "College Readiness vs Socioeconomically Disadvantaged Students",
    x = "Percentage of Economically Disadvantaged Students",
    y = "College Readiness Rate (%)",
    caption = "Source: Bexar County School Districts"
  ) +theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

model <- lm(DA0CC21R ~ DPETECOP, data = clean_data)
summary(model)
## 
## Call:
## lm(formula = DA0CC21R ~ DPETECOP, data = clean_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.887 -10.962  -1.040   5.053  37.508 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  65.1533     6.9964   9.312 6.40e-10 ***
## DPETECOP     -0.6126     0.1022  -5.996 2.14e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15 on 27 degrees of freedom
## Multiple R-squared:  0.5711, Adjusted R-squared:  0.5552 
## F-statistic: 35.95 on 1 and 27 DF,  p-value: 2.139e-06
summary(model)
## 
## Call:
## lm(formula = DA0CC21R ~ DPETECOP, data = clean_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.887 -10.962  -1.040   5.053  37.508 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  65.1533     6.9964   9.312 6.40e-10 ***
## DPETECOP     -0.6126     0.1022  -5.996 2.14e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15 on 27 degrees of freedom
## Multiple R-squared:  0.5711, Adjusted R-squared:  0.5552 
## F-statistic: 35.95 on 1 and 27 DF,  p-value: 2.139e-06
confint(model)
##                  2.5 %     97.5 %
## (Intercept) 50.7978617 79.5087292
## DPETECOP    -0.8221701 -0.4029321
plot(model, which = 2)

raintest(model)
## 
##  Rainbow test
## 
## data:  model
## Rain = 0.56078, df1 = 15, df2 = 12, p-value = 0.8557
plot(model, which = 1)

durbinWatsonTest(model)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.1892619      1.596642   0.224
##  Alternative hypothesis: rho != 0
bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 13.262, df = 1, p-value = 0.0002708
plot(model, which = 3)

shapiro.test(residuals(model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.97464, p-value = 0.6904
plot(model, which = 2)

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

par(mfrow = c(1, 1))
std_resid <- rstandard(model)
cooks_dist <- cooks.distance(model)
outliers <- clean_data[abs(std_resid) > 2, ]
print("Potential outliers:")
## [1] "Potential outliers:"
print(outliers)
##              DISTNAME DPETECOP DA0CC21R
## 13        BASIS TEXAS      8.1     97.7
## 18 RANDOLPH FIELD ISD      3.7     35.0