library(ggplot2)
library(janitor)
##
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
Animalcare <- read.csv("Desktop/Applied Quantitative Methods/Animal_Care_and_Control_Division_Annual_Statistics-1.csv")
head(Animalcare)
## Year Number.of.Employees Number.of.Division.Vehicles Annual.Budget
## 1 2004 15.150 3 782931
## 2 2005 16.600 3 760206
## 3 2006 16.600 3 1030661
## 4 2007 17.230 3 1062946
## 5 2008 17.725 3 1132507
## 6 2009 17.975 3 1164248
## Owner.Surrenders Strays Impounds.by.ACO..Added.in.2015.
## 1 2351 3257 NA
## 2 2104 3038 NA
## 3 2361 2800 NA
## 4 2294 2519 NA
## 5 2131 2690 NA
## 6 2203 2316 NA
## Total.Intake.of.Animals Adoptions Return.to.Owner Euthanized
## 1 5608 1896 558 2277
## 2 5142 1866 542 1724
## 3 5161 1737 554 1856
## 4 4813 1825 523 1753
## 5 4821 1849 544 1721
## 6 4591 1893 452 1612
## Transported.to.other.shelters.and.rescues Fostered.Animals Service.Calls
## 1 592 NA NA
## 2 670 380 1525
## 3 657 700 2900
## 4 539 800 2950
## 5 431 600 2700
## 6 275 712 2515
## Emergency.Call.Outs Grants.Received Annual.Adoption.Revenue
## 1 150 0 13146
## 2 150 18925 112649
## 3 150 5555 105401
## 4 150 0 100170
## 5 150 11215 106627
## 6 97 37498 112188
names(Animalcare)
## [1] "Year"
## [2] "Number.of.Employees"
## [3] "Number.of.Division.Vehicles"
## [4] "Annual.Budget"
## [5] "Owner.Surrenders"
## [6] "Strays"
## [7] "Impounds.by.ACO..Added.in.2015."
## [8] "Total.Intake.of.Animals"
## [9] "Adoptions"
## [10] "Return.to.Owner"
## [11] "Euthanized"
## [12] "Transported.to.other.shelters.and.rescues"
## [13] "Fostered.Animals"
## [14] "Service.Calls"
## [15] "Emergency.Call.Outs"
## [16] "Grants.Received"
## [17] "Annual.Adoption.Revenue"
Animalcare <- clean_names(Animalcare)
names(Animalcare)
## [1] "year"
## [2] "number_of_employees"
## [3] "number_of_division_vehicles"
## [4] "annual_budget"
## [5] "owner_surrenders"
## [6] "strays"
## [7] "impounds_by_aco_added_in_2015"
## [8] "total_intake_of_animals"
## [9] "adoptions"
## [10] "return_to_owner"
## [11] "euthanized"
## [12] "transported_to_other_shelters_and_rescues"
## [13] "fostered_animals"
## [14] "service_calls"
## [15] "emergency_call_outs"
## [16] "grants_received"
## [17] "annual_adoption_revenue"
model <- lm(annual_budget ~ number_of_employees + total_intake_of_animals + adoptions + annual_adoption_revenue,data = Animalcare)
summary(model)
##
## Call:
## lm(formula = annual_budget ~ number_of_employees + total_intake_of_animals +
## adoptions + annual_adoption_revenue, data = Animalcare)
##
## Residuals:
## Min 1Q Median 3Q Max
## -146586 -71386 7507 35198 233012
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.397e+06 7.254e+05 -3.304 0.00419 **
## number_of_employees 2.070e+05 2.589e+04 7.998 3.66e-07 ***
## total_intake_of_animals -1.259e+02 5.277e+01 -2.385 0.02897 *
## adoptions 4.290e+02 1.170e+02 3.668 0.00191 **
## annual_adoption_revenue -3.207e+00 1.238e+00 -2.592 0.01900 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 103700 on 17 degrees of freedom
## Multiple R-squared: 0.9488, Adjusted R-squared: 0.9368
## F-statistic: 78.78 on 4 and 17 DF, p-value: 9.662e-11
plot(model, which = 1)
raintest(model)
##
## Rainbow test
##
## data: model
## Rain = 2.698, df1 = 11, df2 = 6, p-value = 0.1172
durbinWatsonTest(model)
## lag Autocorrelation D-W Statistic p-value
## 1 0.2829313 1.430979 0.042
## Alternative hypothesis: rho != 0
plot(model, which = 3)
bptest(model)
##
## studentized Breusch-Pagan test
##
## data: model
## BP = 2.8005, df = 4, p-value = 0.5917
plot(model, which = 2)
shapiro.test(residuals(model))
##
## Shapiro-Wilk normality test
##
## data: residuals(model)
## W = 0.94025, p-value = 0.2003
vif(model)
## number_of_employees total_intake_of_animals adoptions
## 2.515897 3.189513 1.642012
## annual_adoption_revenue
## 1.248496
cor(Animalcare[, c("number_of_employees", "total_intake_of_animals",
"adoptions", "annual_adoption_revenue")])
## number_of_employees total_intake_of_animals adoptions
## number_of_employees 1.0000000 -0.7505948 0.3924993
## total_intake_of_animals -0.7505948 1.0000000 -0.6162405
## adoptions 0.3924993 -0.6162405 1.0000000
## annual_adoption_revenue 0.4430212 -0.3672911 0.1989204
## annual_adoption_revenue
## number_of_employees 0.4430212
## total_intake_of_animals -0.3672911
## adoptions 0.1989204
## annual_adoption_revenue 1.0000000