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(readxl)
library(dplyr)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(MASS)
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
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
wd<-getwd()
ECD<-read_xlsx("ECD.xlsx")
###2) Create a linear model "lm()" from the variables, with a continuous dependent variable as the outcome

ECD_multiple<-lm(ECD$`2023 Total Usage` ~ ECD$`Total people in household` + ECD$`Annual Household Income`, data=ECD)

###3) Check the following assumptions:(a-e)

###a)Linearity (plot and raintest)

plot(ECD_multiple, which = 1)

###PLOT: POINTS ARE SCATTERED SOMEWHAT AROUND THE STRAIGHT LINE. ALTHOUGH THERE ARE DEVIATIONS, THE LINE IS MOSTLY STRAIGHT, SUGGESTING THAT THE RELATIONSHIP BETWEEN VARIABLES IS MOSTLY LINEAR- ASSUMPTION MET. 
raintest(ECD_multiple)
## 
##  Rainbow test
## 
## data:  ECD_multiple
## Rain = 0.62729, df1 = 150, df2 = 146, p-value = 0.9976
###RAINBOW TEST: HERE THE P VALUE IS HIGH, MEANING THE DATA IS LINEAR- ASSUMPTION MET.

###b)Independence of errors (durbin-watson)

durbinWatsonTest(ECD_multiple)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.1444348      1.709165   0.014
##  Alternative hypothesis: rho != 0
###DURBIN-WATSON TEST: HERE THE P-VALUE IS LESS THAN 0.5 (0.01) ,THIS MEANS THE ERRORS ARE NOT INDEPENDENT. HOWEVER THE D-W STATISTIC IS   1.709165, WHICH IS VERY CLOSE TO 2 MEANING ASSUMPTION IS MOST LIKELY MET.

###c)Homoscedasticity (plot, bptest)

plot(ECD_multiple,which=3)

###FITTED PLOT: LINE IS PRETTY STRAIGHT AND THE POINTS CLUSTERED SOMEWHAT EVENLY AROUND IT, I WOULD SAY THIS MEETS THE ASSUMPTIONS. 
bptest(ECD_multiple)
## 
##  studentized Breusch-Pagan test
## 
## data:  ECD_multiple
## BP = 1.308, df = 2, p-value = 0.52
###BP TEST: HERE THE P VALUE IS GREATER THAN 0.05 (P-VALUE= 0.52). SINCE THE P VALUE IS NOT SIGNIFICANT, HETEROSCEDASCITY CAN BE REJECTED AND HOMOSCEDASTICITY CAN BE ASSUMED- ASSUMPTION MET

###d)Normality of residuals (QQ plot, shapiro test)

plot(ECD_multiple,which=2)

###QQ PLOT: HERE THERE ARE OBSERVATIONS THAT DO NOT LIE ALONG THE DOTTED LINE, ASSUMPTION IS VIOLATED
shapiro.test(ECD_multiple$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  ECD_multiple$residuals
## W = 0.85562, p-value = 4.706e-16
###SHAPIRO WILK TEST: TEST SHOWS P VALUE IS LESS THAN 0.05 MEANING THE DATA IS PROBABLY NOT NORMAL AND ASSUMPTION IS VIOLATED

###e)No multicolinarity

vif(ECD_multiple)
## ECD$`Total people in household`   ECD$`Annual Household Income` 
##                        1.004333                        1.004333
###VIF: HERE THERE IS NO VIF THAT IS GREATER THAN 10 WHICH COULD  MEAN THEY ARE NOT CORRELATED WITH OTHER VARIABLES- ASSUMPTION IS MET
###5) If your model violates an assumption, which one?
###MODEL VIOLATES: NORMALITY OF RESIDUALS

###6) What would you do to mitigate this assumption? Show your work.

###IN ORDER TO MITIGATE I WOULD ADJUST MODEL USING LOG 

###HAVE TO FILTER DUE TO 0'S IN VARIABLE "ANNUAL HOUSEHOLD INCOME"
ECD_filtered <- ECD %>% filter(`2023 Total Usage` > 0, `Total people in household` > 0, `Annual Household Income` > 0)

ECD_multiple_log <- lm(log(`2023 Total Usage`) ~ log(`Total people in household`) + log(`Annual Household Income`), data = ECD_filtered)

plot(ECD_multiple, which=2)

plot(ECD_multiple_log, which=2)