Version 2
Melanie Loaiza
2024-04-12
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, unionlibrary(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectdata <- read.csv('data.csv' ,
header = TRUE ,
sep=',') options(width=120)
head(data)## fYear month areaResponsability state demographic citizenship encounterCount
## 1 2021 DEC Big Bend Sector TX FMUA CUBA 2
## 2 2021 DEC Big Bend Sector TX FMUA CUBA 2
## 3 2021 DEC Big Bend Sector TX FMUA ECUADOR 2
## 4 2021 DEC Big Bend Sector TX FMUA EL SALVADOR 2
## 5 2021 DEC Big Bend Sector TX FMUA GUATEMALA 3
## 6 2021 DEC Big Bend Sector TX FMUA HONDURAS 15Convert into factors
Area of responsability
is.factor(data\(areaResponsability) data\)areaResponsability = factor(data$areaResponsability)
State
is.factor(data\(state) data\)state = factor(data$state)
Demographic is.factor(data\(demographic) data\)demographic = factor(data$demographic)
Citizenship
is.factor(data\(citizenship) data\)citizenship = factor(data$citizenship)
Anova raw data
model <- aov(data$encounterCount~data$fYear+data$month+data$areaResponsability+data$demographic+data$citizenship, data=data)
summary(model) ## Df Sum Sq Mean Sq F value Pr(>F)
## data$fYear 1 4.193e+07 41926653 30.285 3.78e-08 ***
## data$month 11 6.319e+07 5744172 4.149 3.81e-06 ***
## data$areaResponsability 8 7.012e+08 87648713 63.311 < 2e-16 ***
## data$demographic 2 5.378e+08 268887148 194.226 < 2e-16 ***
## data$citizenship 21 2.269e+09 108051507 78.049 < 2e-16 ***
## Residuals 17783 2.462e+10 1384405
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here, I created the areasize variable, which contains the area of each sector in square kilometers. Importantly, for the Tucson sector, the patrol consists of only 262 border miles. Therefore, I estimated its area by calculating the average of the Yuma Sector and El Paso, considering geographic similarities.
##################### RATE
data <- data %>%
mutate(areasize = case_when(
areaResponsability == "Big Bend Sector" ~ 427747.21,
areaResponsability == "Del Rio Sector" ~ 137432.64,
areaResponsability == "El Centro Sector" ~ 1061.90,
areaResponsability == "El Paso Sector" ~ 325043.75,
areaResponsability == "Laredo Sector" ~ 217665.35,
areaResponsability == "Rio Grande Valley Sector" ~ 88059.66,
areaResponsability == "San Diego Sector" ~ 18129.93,
areaResponsability == "Tucson Sector" ~ 650485.70,
areaResponsability == "Yuma Sector" ~ 470523.48,
TRUE ~ NA_real_ # This is the default case if none of the above conditions are met
))head(data)## fYear month areaResponsability state demographic citizenship encounterCount areasize
## 1 2021 DEC Big Bend Sector TX FMUA CUBA 2 427747.2
## 2 2021 DEC Big Bend Sector TX FMUA CUBA 2 427747.2
## 3 2021 DEC Big Bend Sector TX FMUA ECUADOR 2 427747.2
## 4 2021 DEC Big Bend Sector TX FMUA EL SALVADOR 2 427747.2
## 5 2021 DEC Big Bend Sector TX FMUA GUATEMALA 3 427747.2
## 6 2021 DEC Big Bend Sector TX FMUA HONDURAS 15 427747.2count <- sum(is.na(data$encounters.areasize))
print(count)## [1] 0The Encounters.areasize variable has been created and transformed into logencounters.areasize. However, the distribution remains strongly right-skewed even after the log transformation.
options(scipen = 999)
data <- data %>%
mutate(encounters.areasize = encounterCount / areasize)
head(data)## fYear month areaResponsability state demographic citizenship encounterCount areasize encounters.areasize
## 1 2021 DEC Big Bend Sector TX FMUA CUBA 2 427747.2 0.000004675659
## 2 2021 DEC Big Bend Sector TX FMUA CUBA 2 427747.2 0.000004675659
## 3 2021 DEC Big Bend Sector TX FMUA ECUADOR 2 427747.2 0.000004675659
## 4 2021 DEC Big Bend Sector TX FMUA EL SALVADOR 2 427747.2 0.000004675659
## 5 2021 DEC Big Bend Sector TX FMUA GUATEMALA 3 427747.2 0.000007013488
## 6 2021 DEC Big Bend Sector TX FMUA HONDURAS 15 427747.2 0.000035067441#log transformation
data$logencounters.areasize <- log(data$encounters.areasize + 1)
hist(data$logencounters.areasize, main = "Histogram of Log-Transformed Encounter Counts", xlab = "Log(Encounter Count + 1)")
boxplot(data$logencounters.areasize, main = "Boxplot for Log-Transformed Encounter Counts") 
Model 2
model2 <- aov(data$logencounters.areasize~data$fYear+data$month+data$areaResponsability+data$demographic+data$citizenship, data=data)
summary(model2)## Df Sum Sq Mean Sq F value Pr(>F)
## data$fYear 1 0.01 0.0081 2.072 0.150
## data$month 11 0.03 0.0026 0.665 0.773
## data$areaResponsability 8 6.78 0.8473 217.420 <0.0000000000000002 ***
## data$demographic 2 0.56 0.2807 72.039 <0.0000000000000002 ***
## data$citizenship 21 2.55 0.1216 31.207 <0.0000000000000002 ***
## Residuals 17783 69.30 0.0039
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Assumptions model 2 - log transformation
qqnorm(model2$residuals)
qqline(model2$residuals, col = "red") 
plot(model2, which = 4) 
plot(model2$fitted.values, model2$residuals)
abline(h = 0, col = "red")
hist(model2$residuals, breaks = "FD", main = "Histogram of Residuals", xlab = "Residuals")
Option 1 :
Create indicator variables for year as 2021 reference group
## year
is.factor(data$fYear) ## [1] FALSEdata$fYear = factor(data$fYear)
data$g2022 <- ifelse(data$fYear=='2022', 1, 0)
data$g2023 <- ifelse(data$fYear=='2023', 1, 0)
data$g2024 <- ifelse(data$fYear=='2024', 1, 0)
model3 <- lm(data$logencounters.areasize~data$g2022+data$g2023+data$g2024)
summary(model3)##
## Call:
## lm(formula = data$logencounters.areasize ~ data$g2022 + data$g2023 +
## data$g2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.01176 -0.01121 -0.00936 -0.00809 1.73705
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0090284 0.0009298 9.710 <0.0000000000000002 ***
## data$g2022 0.0027298 0.0012811 2.131 0.0331 *
## data$g2023 0.0022114 0.0013069 1.692 0.0906 .
## data$g2024 0.0017250 0.0018668 0.924 0.3555
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06667 on 17823 degrees of freedom
## Multiple R-squared: 0.0002815, Adjusted R-squared: 0.0001132
## F-statistic: 1.673 on 3 and 17823 DF, p-value: 0.1705Create indicator variables for month . Overall model as factor.month not significant
#month
is.factor(data$month) ## [1] FALSEdata$month = factor(data$month)
model4 <- aov(data$logencounters.areasize~data$month, data=data)
summary(model4) ## Df Sum Sq Mean Sq F value Pr(>F)
## data$month 11 0.02 0.002142 0.482 0.916
## Residuals 17815 79.21 0.004446data$gFEB <- ifelse(data$month=='FEB', 1, 0)
data$gMAR <- ifelse(data$month=='MAR', 1, 0)
data$gAPR <- ifelse(data$month=='APR', 1, 0)
data$gJUN <- ifelse(data$month=='JUN', 1, 0)
data$gJUL <- ifelse(data$month=='JUL', 1, 0)
data$gAUG <- ifelse(data$month=='AUG', 1, 0)
data$gSEP <- ifelse(data$month=='SEP', 1, 0)
data$gOCT <- ifelse(data$month=='OCT', 1, 0)
data$gNOV <- ifelse(data$month=='NOV', 1, 0)
data$gDEC <- ifelse(data$month=='DEC', 1, 0)
model5 <- lm(data$logencounters.areasize~data$gFEB+data$gDEC)
summary(model5)##
## Call:
## lm(formula = data$logencounters.areasize ~ data$gFEB + data$gDEC)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.01171 -0.01079 -0.01042 -0.00807 1.73526
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0108214 0.0005555 19.481 <0.0000000000000002 ***
## data$gFEB -0.0019538 0.0017115 -1.142 0.254
## data$gDEC 0.0008941 0.0016982 0.527 0.599
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06667 on 17824 degrees of freedom
## Multiple R-squared: 9.691e-05, Adjusted R-squared: -1.529e-05
## F-statistic: 0.8637 on 2 and 17824 DF, p-value: 0.4216Option 2 :
datatime <- data.frame(
fYear = c(rep(2021:2023, each=12), rep(2024, each=5)), # Full years 2021 to 2023, partial for 2024
month = c(rep(c("JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP", "OCT", "NOV", "DEC"), times=3),
c("JAN", "FEB", "OCT" , "NOV" , "DEC"))
)
head(datatime)## fYear month
## 1 2021 JAN
## 2 2021 FEB
## 3 2021 MAR
## 4 2021 APR
## 5 2021 MAY
## 6 2021 JUNdatatime$format <- as.Date(paste(datatime$fYear, datatime$month, "01", sep="-"), format="%Y-%b-%d")
refdate <- as.Date("2021-01-01")
datatime$timesequence <- as.numeric(format(datatime$format, "%Y%m")) - as.numeric(format(refdate, "%Y%m"))
data <- merge(data, datatime[, c("fYear", "month", "timesequence")],
by = c("fYear", "month"), all.x = TRUE)
head(data)## fYear month areaResponsability state demographic citizenship encounterCount areasize encounters.areasize
## 1 2021 APR Big Bend Sector TX FMUA ECUADOR 10 427747.2 0.000023378294
## 2 2021 APR Big Bend Sector TX FMUA EL SALVADOR 20 427747.2 0.000046756588
## 3 2021 APR Big Bend Sector TX FMUA GUATEMALA 15 427747.2 0.000035067441
## 4 2021 APR Big Bend Sector TX FMUA HONDURAS 49 427747.2 0.000114553640
## 5 2021 APR Big Bend Sector TX FMUA HONDURAS 3 427747.2 0.000007013488
## 6 2021 APR Big Bend Sector TX FMUA MEXICO 47 427747.2 0.000109877981
## logencounters.areasize g2022 g2023 g2024 gFEB gMAR gAPR gJUN gJUL gAUG gSEP gOCT gNOV gDEC timesequence
## 1 0.000023378021 0 0 0 0 0 1 0 0 0 0 0 0 0 3
## 2 0.000046755495 0 0 0 0 0 1 0 0 0 0 0 0 0 3
## 3 0.000035066826 0 0 0 0 0 1 0 0 0 0 0 0 0 3
## 4 0.000114547079 0 0 0 0 0 1 0 0 0 0 0 0 0 3
## 5 0.000007013464 0 0 0 0 0 1 0 0 0 0 0 0 0 3
## 6 0.000109871945 0 0 0 0 0 1 0 0 0 0 0 0 0 3model6 <- aov(data$logencounters.areasize~data$timesequence+data$demographic+data$citizenship, data=data)
summary(model6) ## Df Sum Sq Mean Sq F value Pr(>F)
## data$timesequence 1 0.01 0.00866 2.023 0.155
## data$demographic 2 0.57 0.28634 66.892 <0.0000000000000002 ***
## data$citizenship 21 2.45 0.11659 27.238 <0.0000000000000002 ***
## Residuals 17802 76.20 0.00428
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model7 <- lm(data$logencounters.areasize~data$timesequence, data=data)
summary(model7) ##
## Call:
## lm(formula = data$logencounters.areasize ~ data$timesequence,
## data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.01206 -0.01059 -0.00984 -0.00826 1.73623
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.009816776 0.000818535 11.993 <0.0000000000000002 ***
## data$timesequence 0.000007229 0.000005179 1.396 0.163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06667 on 17825 degrees of freedom
## Multiple R-squared: 0.0001093, Adjusted R-squared: 5.32e-05
## F-statistic: 1.948 on 1 and 17825 DF, p-value: 0.1628qqnorm(model6$residuals)
qqline(model6$residuals, col = "red") 
plot(model6, which = 4) 
plot(model6$fitted.values, model6$residuals)
abline(h = 0, col = "red")
hist(model6$residuals, breaks = "FD", main = "Histogram of Residuals", xlab = "Residuals")