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
#2. Create a linear model "lm()" from the variables, with a continuous dependent variable as the outcome.
data_multiple <- lm(DDH00A001S22R~DPETBILP+DPSTBIFP+DPFPABILP,data=district)
#3. Check the following 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.528
## 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 multicolinarity (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
#4. Does your model meet those assumptions?
#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.
#5. If your model violates an assumption, which one?
#The Shapiro-Wilk test shows a significant p-value, meaning the residuals are not normal.
#6. What would you do to mitigate this assumption?
#A log transformation can be performed for the 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