library(readxl)
district <- read_excel("district.xls")
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(pastecs)
## 
## Attaching package: 'pastecs'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
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
#Linear model "lm()" from the variables.
data_multiple <- lm(DDH00A001S22R~DPETBILP+DPSTBIFP+DPFPABILP,data = district)
#Assumptions: Linearity(plot and raintest)
plot(data_multiple,which = 1)

raintest(data_multiple)
## 
##  Rainbow test
## 
## data:  data_multiple
## Rain = 0.87303, df1 = 599, df2 = 594, p-value = 0.9513
#Independence of errors (durbin-watson)
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
durbinWatsonTest(data_multiple)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.01915531      1.961555     0.5
##  Alternative hypothesis: rho != 0
#Homoscedasticity(plot,bptest)
plot(data_multiple,which = 3)

bptest(data_multiple)
## 
##  studentized Breusch-Pagan test
## 
## data:  data_multiple
## BP = 4.4104, df = 3, p-value = 0.2204
#Normality of residuals (QQ plot, shapiro test)
plot(data_multiple,which = 2)

shapiro.test(data_multiple$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  data_multiple$residuals
## W = 0.89531, p-value < 2.2e-16
#No multicollinearity (VIF, cor)
vif(data_multiple)
##  DPETBILP  DPSTBIFP DPFPABILP 
##  1.989514  1.459168  1.598077
data_multiple_vars <- district %>% dplyr::select(DPETBILP,DPSTBIFP,DPFPABILP)
cor(data_multiple_vars)
##           DPETBILP DPSTBIFP DPFPABILP
## DPETBILP         1       NA        NA
## DPSTBIFP        NA        1        NA
## DPFPABILP       NA       NA         1
#Log transformation - Normality of Residuals
district_drop <- district %>% drop_na(DDH00A001S22R,DPETBILP,DPSTBIFP,DPFPABILP) %>% filter(DDH00A001S22R>0 & DPETBILP>0)
data_multiple_log <- lm(log(DDH00A001S22R)~log(DPETBILP)+DPSTBIFP+DPFPABILP,data = district_drop)

bptest(data_multiple_log)
## 
##  studentized Breusch-Pagan test
## 
## data:  data_multiple_log
## BP = 0.79281, df = 3, p-value = 0.8512
shapiro.test(data_multiple_log$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  data_multiple_log$residuals
## W = 0.86179, p-value < 2.2e-16
summary(data_multiple)
## 
## Call:
## lm(formula = DDH00A001S22R ~ DPETBILP + DPSTBIFP + DPFPABILP, 
##     data = district)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -74.914  -5.757   0.802   7.769  27.088 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 73.91439    0.47738 154.833  < 2e-16 ***
## DPETBILP    -0.22250    0.03451  -6.448 1.65e-10 ***
## DPSTBIFP     0.09047    0.07542   1.200    0.231    
## DPFPABILP    0.44518    0.34458   1.292    0.197    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.5 on 1193 degrees of freedom
##   (10 observations deleted due to missingness)
## Multiple R-squared:  0.04516,    Adjusted R-squared:  0.04276 
## F-statistic: 18.81 on 3 and 1193 DF,  p-value: 6.377e-12
#Assumptions met-
#The Homoscedasticity test shows that the model has a mostly straight line. 
#The model also meets the vif command. All the variables show a weak correlation of less than 5.
#The durbinWatsonTest shows a p-value greater than 0.05 and the D-W Statistic is very close to 2.
#The bptest shows that the model's p-value is more that 0.05 which is not significant, meaning the model is homoscedastic, not heteroscedastic.
#The rainbow test shows that the model is very linear with a p-value of 0.95.
#Assumptions violated-
#The Shapiro-Wilk test shows a significant p-value, meaning the residuals are not normal. A log transformation was performed for the normality of residuals.