Sasha-Sharp-PAD-6833-Homework-7
Sasha Sharp
2026-11-03
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.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 errorslibrary(readxl)
library(MASS)##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## selectlibrary(lmtest)## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(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':
##
## somegov_data <- read.csv("GOVSEMPTIMESERIES.GS00EMP01-2026-09-22T170141.csv")head(gov_data)## Year..time. Geographic.Area.Name..NAME.
## 1 2025 United States
## 2 2025 United States
## 3 2025 United States
## 4 2025 United States
## 5 2025 United States
## 6 2025 United States
## Meaning.of.Aggregate.Description..AGG_DESC_LABEL.
## 1 Total - All Government Employment Functions
## 2 Financial Administration
## 3 Other Government Administration
## 4 Judicial and Legal
## 5 Police Protection Total
## 6 Police Protection - Persons with Power of Arrest
## Meaning.of.Type.of.Government..GOVTYPE_LABEL. Full.Time.Employment..FT_EMP.
## 1 State and Local 15,790,344
## 2 State and Local 429,784
## 3 State and Local 270,830
## 4 State and Local 424,793
## 5 State and Local 936,007
## 6 State and Local 709,292
## Full.Time.Payroll..FT_PAY. Part.Time.Employment..PT_EMP.
## 1 106,743,344,645 4,524,442
## 2 2,992,697,403 53,004
## 3 1,888,513,707 180,646
## 4 3,161,227,083 37,611
## 5 7,972,096,680 87,506
## 6 6,477,500,863 33,677
## Part.Time.Payroll..PT_PAY. Part.Time.Hours..PT_HRS.
## 1 8,367,667,849 N
## 2 88,775,081 N
## 3 157,914,189 N
## 4 87,901,725 N
## 5 146,523,037 N
## 6 59,852,258 N
## Full.Time.Equivalent.Employment..FTE.
## 1 17,527,046
## 2 448,891
## 3 301,610
## 4 440,423
## 5 970,667
## 6 723,157
## Total.Full.Time.and.Part.Time.Employment..TOT_EMP.
## 1 20,314,786
## 2 482,788
## 3 451,476
## 4 462,404
## 5 1,023,513
## 6 742,969
## Total.Full.Time.and.Part.Time.Payroll..TOT_PAY.
## 1 115,111,012,494
## 2 3,081,472,484
## 3 2,046,427,896
## 4 3,249,128,808
## 5 8,118,619,717
## 6 6,537,353,121
## Full.Time.Employment.Coefficient.of.Variation..FT_EMP_CV.
## 1 0.38
## 2 1.66
## 3 2.55
## 4 0.65
## 5 1.93
## 6 2.35
## Full.Time.Payroll.Coefficient.of.Variation..FT_PAY_CV.
## 1 0.34
## 2 1.35
## 3 3.29
## 4 0.96
## 5 1.71
## 6 1.98
## Part.Time.Employment.Coefficient.of.Variation..PT_EMP_CV.
## 1 0.42
## 2 9.59
## 3 2.50
## 4 1.38
## 5 1.87
## 6 4.18
## Part.Time.Payroll.Coefficient.of.Variation..PT_PAY_CV.
## 1 1.42
## 2 6.14
## 3 8.68
## 4 1.77
## 5 4.41
## 6 10.59
## Part.Time.Hours.Coefficient.of.Variation..PT_HRS_CV.
## 1 N
## 2 N
## 3 N
## 4 N
## 5 N
## 6 N
## Full.Time.Equivalent.Employment.Coefficient.of.Variation..FTE_CV.
## 1 0.36
## 2 1.95
## 3 2.73
## 4 0.64
## 5 1.87
## 6 2.32
## Total.Full.Time.and.Part.Time.Employment.Coefficient.of.Variation..TOT_EMP_CV.
## 1 0.31
## 2 1.81
## 3 1.83
## 4 0.60
## 5 1.78
## 6 2.26
## Total.Full.Time.and.Part.Time.Payroll.Coefficient.of.Variation..TOT_PAY_CV.
## 1 0.33
## 2 1.32
## 3 3.11
## 4 0.94
## 5 1.69
## 6 1.96names(gov_data)## [1] "Year..time."
## [2] "Geographic.Area.Name..NAME."
## [3] "Meaning.of.Aggregate.Description..AGG_DESC_LABEL."
## [4] "Meaning.of.Type.of.Government..GOVTYPE_LABEL."
## [5] "Full.Time.Employment..FT_EMP."
## [6] "Full.Time.Payroll..FT_PAY."
## [7] "Part.Time.Employment..PT_EMP."
## [8] "Part.Time.Payroll..PT_PAY."
## [9] "Part.Time.Hours..PT_HRS."
## [10] "Full.Time.Equivalent.Employment..FTE."
## [11] "Total.Full.Time.and.Part.Time.Employment..TOT_EMP."
## [12] "Total.Full.Time.and.Part.Time.Payroll..TOT_PAY."
## [13] "Full.Time.Employment.Coefficient.of.Variation..FT_EMP_CV."
## [14] "Full.Time.Payroll.Coefficient.of.Variation..FT_PAY_CV."
## [15] "Part.Time.Employment.Coefficient.of.Variation..PT_EMP_CV."
## [16] "Part.Time.Payroll.Coefficient.of.Variation..PT_PAY_CV."
## [17] "Part.Time.Hours.Coefficient.of.Variation..PT_HRS_CV."
## [18] "Full.Time.Equivalent.Employment.Coefficient.of.Variation..FTE_CV."
## [19] "Total.Full.Time.and.Part.Time.Employment.Coefficient.of.Variation..TOT_EMP_CV."
## [20] "Total.Full.Time.and.Part.Time.Payroll.Coefficient.of.Variation..TOT_PAY_CV."gov_data <- gov_data %>%
rename(
PT_PAY = Part.Time.Payroll..PT_PAY.,
PT_EMP = Part.Time.Employment..PT_EMP.,
FT_EMP = Full.Time.Employment..FT_EMP.
)names(gov_data)## [1] "Year..time."
## [2] "Geographic.Area.Name..NAME."
## [3] "Meaning.of.Aggregate.Description..AGG_DESC_LABEL."
## [4] "Meaning.of.Type.of.Government..GOVTYPE_LABEL."
## [5] "FT_EMP"
## [6] "Full.Time.Payroll..FT_PAY."
## [7] "PT_EMP"
## [8] "PT_PAY"
## [9] "Part.Time.Hours..PT_HRS."
## [10] "Full.Time.Equivalent.Employment..FTE."
## [11] "Total.Full.Time.and.Part.Time.Employment..TOT_EMP."
## [12] "Total.Full.Time.and.Part.Time.Payroll..TOT_PAY."
## [13] "Full.Time.Employment.Coefficient.of.Variation..FT_EMP_CV."
## [14] "Full.Time.Payroll.Coefficient.of.Variation..FT_PAY_CV."
## [15] "Part.Time.Employment.Coefficient.of.Variation..PT_EMP_CV."
## [16] "Part.Time.Payroll.Coefficient.of.Variation..PT_PAY_CV."
## [17] "Part.Time.Hours.Coefficient.of.Variation..PT_HRS_CV."
## [18] "Full.Time.Equivalent.Employment.Coefficient.of.Variation..FTE_CV."
## [19] "Total.Full.Time.and.Part.Time.Employment.Coefficient.of.Variation..TOT_EMP_CV."
## [20] "Total.Full.Time.and.Part.Time.Payroll.Coefficient.of.Variation..TOT_PAY_CV."gov_data <- gov_data %>%
mutate(
PT_PAY = as.numeric(gsub(",", "", PT_PAY)),
PT_EMP = as.numeric(gsub(",", "", PT_EMP)),
FT_EMP = as.numeric(gsub(",", "", FT_EMP))
)str(gov_data)## 'data.frame': 1976 obs. of 20 variables:
## $ Year..time. : int 2025 2025 2025 2025 2025 2025 2025 2025 2025 2025 ...
## $ Geographic.Area.Name..NAME. : chr "United States" "United States" "United States" "United States" ...
## $ Meaning.of.Aggregate.Description..AGG_DESC_LABEL. : chr "Total - All Government Employment Functions" "Financial Administration" "Other Government Administration" "Judicial and Legal" ...
## $ Meaning.of.Type.of.Government..GOVTYPE_LABEL. : chr "State and Local" "State and Local" "State and Local" "State and Local" ...
## $ FT_EMP : num 15790344 429784 270830 424793 936007 ...
## $ Full.Time.Payroll..FT_PAY. : chr "106,743,344,645" "2,992,697,403" "1,888,513,707" "3,161,227,083" ...
## $ PT_EMP : num 4524442 53004 180646 37611 87506 ...
## $ PT_PAY : num 8.37e+09 8.88e+07 1.58e+08 8.79e+07 1.47e+08 ...
## $ Part.Time.Hours..PT_HRS. : chr "N" "N" "N" "N" ...
## $ Full.Time.Equivalent.Employment..FTE. : chr "17,527,046" "448,891" "301,610" "440,423" ...
## $ Total.Full.Time.and.Part.Time.Employment..TOT_EMP. : chr "20,314,786" "482,788" "451,476" "462,404" ...
## $ Total.Full.Time.and.Part.Time.Payroll..TOT_PAY. : chr "115,111,012,494" "3,081,472,484" "2,046,427,896" "3,249,128,808" ...
## $ Full.Time.Employment.Coefficient.of.Variation..FT_EMP_CV. : num 0.38 1.66 2.55 0.65 1.93 2.35 3.09 2.71 3.04 3.35 ...
## $ Full.Time.Payroll.Coefficient.of.Variation..FT_PAY_CV. : num 0.34 1.35 3.29 0.96 1.71 1.98 3.2 3.14 3.49 3.47 ...
## $ Part.Time.Employment.Coefficient.of.Variation..PT_EMP_CV. : num 0.42 9.59 2.5 1.38 1.87 ...
## $ Part.Time.Payroll.Coefficient.of.Variation..PT_PAY_CV. : num 1.42 6.14 8.68 1.77 4.41 ...
## $ Part.Time.Hours.Coefficient.of.Variation..PT_HRS_CV. : chr "N" "N" "N" "N" ...
## $ Full.Time.Equivalent.Employment.Coefficient.of.Variation..FTE_CV. : num 0.36 1.95 2.73 0.64 1.87 2.32 2.85 3.69 4.14 5.3 ...
## $ Total.Full.Time.and.Part.Time.Employment.Coefficient.of.Variation..TOT_EMP_CV.: chr "0.31" "1.81" "1.83" "0.60" ...
## $ Total.Full.Time.and.Part.Time.Payroll.Coefficient.of.Variation..TOT_PAY_CV. : chr "0.33" "1.32" "3.11" "0.94" ...gov_model <- lm(PT_PAY ~ PT_EMP + FT_EMP, data=gov_data)
summary(gov_model)##
## Call:
## lm(formula = PT_PAY ~ PT_EMP + FT_EMP, data = gov_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -472561885 -603608 206099 332698 536793227
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.327e+05 6.167e+05 -0.539 0.59
## PT_EMP 1.720e+03 1.488e+01 115.584 <2e-16 ***
## FT_EMP 4.152e+01 4.654e+00 8.920 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27240000 on 1973 degrees of freedom
## Multiple R-squared: 0.9893, Adjusted R-squared: 0.9893
## F-statistic: 9.119e+04 on 2 and 1973 DF, p-value: < 2.2e-16plot(gov_model, which=1)
raintest(gov_model)##
## Rainbow test
##
## data: gov_model
## Rain = 24.764, df1 = 988, df2 = 985, p-value < 2.2e-16Linearity: The linearity assumption was not met. The Rainbow test produced a p-value below 0.05 (p < 2.2e-16), indicating a significant departure from linearity. The Residuals vs. Fitted plot also showed a downward-sloping red line and several extreme observations. These results suggest that the relationship between the independent variables and the dependent variable may not be adequately represented by a linear model.
durbinWatsonTest(gov_model)## lag Autocorrelation D-W Statistic p-value
## 1 0.000473748 1.995788 0.72
## Alternative hypothesis: rho != 0Independence of Errors: The Durbin-Watson test did not detect significant autocorrelation in the original model. The test produced a Durbin-Watson statistic of 1.996 and a p-value of 0.714. Since the p-value is greater than 0.05 and the statistic is close to 2, there is insufficient evidence of autocorrelation. However, because the dataset contains geographic and government-level observations rather than a chronological time series, this is a result that should be interpreted cautiously.
plot(gov_model, which=3)
bptest(gov_model)##
## studentized Breusch-Pagan test
##
## data: gov_model
## BP = 226.27, df = 2, p-value < 2.2e-16Homoscedasticity: The homoscedasticity assumption was not met. The Breusch-Pagan test produced a p-value below 0.05 (p < 2.2e-16), indicating evidence of heteroscedasticity. The Scale-Location plot also showed a sharply increasing red line and an uneven distribution of residuals. These results suggest that the variance of the model’s errors is not constant across fitted values.
plot(gov_model, which=2)
shapiro.test(gov_model$residuals)##
## Shapiro-Wilk normality test
##
## data: gov_model$residuals
## W = 0.19564, p-value < 2.2e-16Normality of Residuals: The normality assumption was not met. The Shapiro-Wilk test produced a p-value below 0.05 (p < 2.2e-16), indicating that the residuals significantly differ from a normal distribution. The Q-Q plot also showed substantial deviations from the diagonal reference line, particularly at both ends of the distribution. These findings suggest that the model’s residuals are not normally distributed.
vif(gov_model)## PT_EMP FT_EMP
## 11.8334 11.8334cor(gov_data %>% dplyr::select(PT_EMP, FT_EMP),
use="complete.obs")## PT_EMP FT_EMP
## PT_EMP 1.0000000 0.9568142
## FT_EMP 0.9568142 1.0000000No Multicollinearity: The assumption of no multicollinearity was not met. The VIF results were 11.83 for both PT_EMP and FT_EMP, exceeding the recommended threshold of 10. The correlation matrix also showed a strong positive correlation of 0.957 between the two independent variables. These findings indicate substantial multicollinearity, suggesting that the independent variables are highly correlated with one another.
gov_model_revised <- lm(PT_PAY ~ PT_EMP, data=gov_data)
summary(gov_model_revised)##
## Call:
## lm(formula = PT_PAY ~ PT_EMP, data = gov_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -550605725 -454701 221939 310090 551110625
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.599e+05 6.288e+05 -0.413 0.679
## PT_EMP 1.847e+03 4.411e+00 418.713 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27780000 on 1974 degrees of freedom
## Multiple R-squared: 0.9889, Adjusted R-squared: 0.9889
## F-statistic: 1.753e+05 on 1 and 1974 DF, p-value: < 2.2e-16Multicollinearity Mitigation: To address multicollinearity, a revised regression model was created by removing FT_EMP and retaining PT_EMP as the independent variable. This decision was based on the strong correlation between the two original predictors and the relevance of part-time employment to the dependent variable, part-time payroll. Because the revised model contains only one independent variable, multicollinearity between predictors is no longer a concern. However, removing FT_EMP changes the model specification and does not automatically resolve the other regression assumption violations.
sapply(gov_data[, c("PT_PAY", "PT_EMP", "FT_EMP")],
function(x) sum(x <= 0, na.rm = TRUE))## PT_PAY PT_EMP FT_EMP
## 107 107 75gov_data_log <- gov_data %>%
filter(
PT_PAY > 0,
PT_EMP > 0,
FT_EMP > 0
)gov_model_log <- lm(
log(PT_PAY) ~ log(PT_EMP) + log(FT_EMP),
data = gov_data_log
)
summary(gov_model_log)##
## Call:
## lm(formula = log(PT_PAY) ~ log(PT_EMP) + log(FT_EMP), data = gov_data_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.39070 -0.28576 -0.00249 0.30819 1.81590
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.550988 0.054271 120.71 <2e-16 ***
## log(PT_EMP) 0.852762 0.008744 97.52 <2e-16 ***
## log(FT_EMP) 0.209633 0.010953 19.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5097 on 1865 degrees of freedom
## Multiple R-squared: 0.9618, Adjusted R-squared: 0.9618
## F-statistic: 2.35e+04 on 2 and 1865 DF, p-value: < 2.2e-16plot(gov_model_log, which = 1)
raintest(gov_model_log)##
## Rainbow test
##
## data: gov_model_log
## Rain = 0.99593, df1 = 934, df2 = 931, p-value = 0.5248Linearity After Transformation: The linearity assumption was met after applying a log transformation to the dependent and independent variables. The Rainbow test produced a p-value of 0.5248, which is greater than 0.05. The Residuals vs. Fitted plot also showed a relatively flat red line, suggesting that the log-transformed model better represents the linear relationship between the variables.
plot(gov_model_log, which = 3)
bptest(gov_model_log)##
## studentized Breusch-Pagan test
##
## data: gov_model_log
## BP = 89.484, df = 2, p-value < 2.2e-16Homoscedasticity After Transformation: The homoscedasticity assumption was not met after applying the log transformation. The Breusch-Pagan test produced a p-value below 0.05 (p < 2.2e-16), indicating evidence of heteroscedasticity. The Scale-Location plot also showed a downward-sloping red line, suggesting that the variance of the residuals is not constant across fitted values. Although the log transformation improved the linearity assumption, it did not fully correct heteroscedasticity.
plot(gov_model_log, which = 2)
shapiro.test(gov_model_log$residuals)##
## Shapiro-Wilk normality test
##
## data: gov_model_log$residuals
## W = 0.98599, p-value = 1.541e-12Normality After Transformation: The normality assumption was not fully met after applying the log transformation. The Shapiro-Wilk test produced a p-value of 1.541e-12, which is below 0.05, indicating that the residuals still differ significantly from a normal distribution. However, the W statistic improved from 0.19564 in the original model to 0.98599 in the transformed model. The Q-Q plot also showed that the residuals followed the reference line more closely, although deviations remained at the tails. This suggests that the log transformation substantially improved residual normality but did not completely resolve the issue.
durbinWatsonTest(gov_model_log)## lag Autocorrelation D-W Statistic p-value
## 1 0.2919066 1.415925 0
## Alternative hypothesis: rho != 0Independence of Errors After Transformation: The Durbin-Watson test indicated a violation of the independence of errors assumption. The test produced a Durbin-Watson statistic of 1.416 and a p-value below 0.001, indicating evidence of positive autocorrelation among the model’s residuals. However, because the dataset contains geographic and government-level observations rather than a chronological time series, the results should be interpreted cautiously.
vif(gov_model_log)## log(PT_EMP) log(FT_EMP)
## 3.601679 3.601679cor(gov_data_log %>% dplyr::select(PT_EMP, FT_EMP),
use = "complete.obs")## PT_EMP FT_EMP
## PT_EMP 1.0000000 0.9567837
## FT_EMP 0.9567837 1.0000000cor(
log(gov_data_log %>% dplyr::select(PT_EMP, FT_EMP)),
use = "complete.obs"
)## PT_EMP FT_EMP
## PT_EMP 1.0000000 0.8499128
## FT_EMP 0.8499128 1.0000000Multicollinearity After Transformation: The log-transformed model showed an improvement in multicollinearity. The VIF values decreased from 11.83 in the original model to 3.60 in the transformed model, falling below the recommended threshold of 5. The correlation between the transformed predictors was 0.850, indicating that a strong positive relationship remains. Overall, the results suggest that problematic multicollinearity was reduced.
gov_model_robust <- MASS::rlm(
log(PT_PAY) ~ log(PT_EMP) + log(FT_EMP),
data = gov_data_log
)
summary(gov_model_robust)##
## Call: rlm(formula = log(PT_PAY) ~ log(PT_EMP) + log(FT_EMP), data = gov_data_log)
## Residuals:
## Min 1Q Median 3Q Max
## -2.45071 -0.28638 -0.01034 0.29373 1.80282
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 6.6593 0.0503 132.4519
## log(PT_EMP) 0.8543 0.0081 105.4598
## log(FT_EMP) 0.1966 0.0101 19.3764
##
## Residual standard error: 0.4295 on 1865 degrees of freedomRobust Regression Mitigation: A robust regression model was created using the rlm() function from the MASS package. The model retained both log-transformed independent variables while reducing the influence of extreme observations. Robust regression provides an alternative approach that is less sensitive to outliers than ordinary least squares regression. However, it does not necessarily resolve the remaining issues with heteroscedasticity, autocorrelation, or non-normal residuals.
plot(gov_model_robust)
shapiro.test(residuals(gov_model_robust))##
## Shapiro-Wilk normality test
##
## data: residuals(gov_model_robust)
## W = 0.98396, p-value = 1.273e-13Normality After Robust Regression: The Shapiro-Wilk test produced a p-value of 1.273e-13, indicating that the residuals still differ significantly from a normal distribution. The W statistic of 0.98396 is substantially closer to 1 than the original model’s W statistic of 0.19564. However, the robust regression did not completely resolve the normality issue. Additionally, the Shapiro-Wilk test can be sensitive to small departures from normality in large samples.
Conclusion: The original regression model did not meet all of the regression assumptions. The results showed problems with linearity, homoscedasticity, normality of the residuals, and multicollinearity. After applying a log transformation, the model improved in linearity, normality, and multicollinearity, but heteroscedasticity was still present. The transformed model also showed evidence of autocorrelation, although this result should be considered carefully because of the type of data being used. I also created a robust regression model to reduce the influence of extreme observations. However, the residuals were still not completely normal. Overall, the changes helped improve the model, but they did not fix every issue with the regression assumptions.