Handout 11

#get the working directory 
getwd()

[1] "/cloud/project"

# make sure the packages for this chapter
# are installed, install if necessary
pkg <- c("ggplot2", "scales", "maptools",
              "sp", "maps", "grid", "car" )
new.pkg <- pkg[!(pkg %in% installed.packages())]
if (length(new.pkg)) {
  install.packages(new.pkg)  
}

Installing packages into ‘/cloud/lib/x86_64-pc-linux-gnu-library/4.3’
(as ‘lib’ is unspecified)

Warning in install.packages :
  package ‘maptools’ is not available for this version of R

A version of this package for your version of R might be available elsewhere,
see the ideas at
https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packages

also installing the dependencies ‘rbibutils’, ‘cowplot’, ‘Deriv’, ‘microbenchmark’, ‘Rdpack’, ‘numDeriv’, ‘doBy’, ‘SparseM’, ‘MatrixModels’, ‘minqa’, ‘nloptr’, ‘reformulas’, ‘Rcpp’, ‘RcppEigen’, ‘carData’, ‘abind’, ‘Formula’, ‘pbkrtest’, ‘quantreg’, ‘lme4’

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/rbibutils_2.3.tar.gz'
Content type 'application/x-gzip' length 1137301 bytes (1.1 MB)
==================================================
downloaded 1.1 MB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/cowplot_1.1.3.tar.gz'
Content type 'application/x-gzip' length 1376043 bytes (1.3 MB)
==================================================
downloaded 1.3 MB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/Deriv_4.1.6.tar.gz'
Content type 'application/x-gzip' length 149977 bytes (146 KB)
==================================================
downloaded 146 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/microbenchmark_1.5.0.tar.gz'
Content type 'application/x-gzip' length 64450 bytes (62 KB)
==================================================
downloaded 62 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/Rdpack_2.6.2.tar.gz'
Content type 'application/x-gzip' length 747892 bytes (730 KB)
==================================================
downloaded 730 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/numDeriv_2016.8-1.1.tar.gz'
Content type 'application/x-gzip' length 112843 bytes (110 KB)
==================================================
downloaded 110 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/doBy_4.6.25.tar.gz'
Content type 'application/x-gzip' length 4831866 bytes (4.6 MB)
==================================================
downloaded 4.6 MB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/SparseM_1.84-2.tar.gz'
Content type 'application/x-gzip' length 879962 bytes (859 KB)
==================================================
downloaded 859 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/MatrixModels_0.5-3.tar.gz'
Content type 'application/x-gzip' length 413182 bytes (403 KB)
==================================================
downloaded 403 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/minqa_1.2.8.tar.gz'
Content type 'application/x-gzip' length 121451 bytes (118 KB)
==================================================
downloaded 118 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/nloptr_2.1.1.tar.gz'
Content type 'application/x-gzip' length 554121 bytes (541 KB)
==================================================
downloaded 541 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/reformulas_0.4.0.tar.gz'
Content type 'application/x-gzip' length 89712 bytes (87 KB)
==================================================
downloaded 87 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/Rcpp_1.0.14.tar.gz'
Content type 'application/x-gzip' length 2176986 bytes (2.1 MB)
==================================================
downloaded 2.1 MB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/RcppEigen_0.3.4.0.2.tar.gz'
Content type 'application/x-gzip' length 1845612 bytes (1.8 MB)
==================================================
downloaded 1.8 MB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/carData_3.0-5.tar.gz'
Content type 'application/x-gzip' length 1820881 bytes (1.7 MB)
==================================================
downloaded 1.7 MB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/abind_1.4-8.tar.gz'
Content type 'application/x-gzip' length 64323 bytes (62 KB)
==================================================
downloaded 62 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/Formula_1.2-5.tar.gz'
Content type 'application/x-gzip' length 158399 bytes (154 KB)
==================================================
downloaded 154 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/pbkrtest_0.5.3.tar.gz'
Content type 'application/x-gzip' length 174550 bytes (170 KB)
==================================================
downloaded 170 KB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/quantreg_6.00.tar.gz'
Content type 'application/x-gzip' length 1446791 bytes (1.4 MB)
==================================================
downloaded 1.4 MB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/lme4_1.1-36.tar.gz'
Content type 'application/x-gzip' length 4235318 bytes (4.0 MB)
==================================================
downloaded 4.0 MB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/maps_3.4.2.1.tar.gz'
Content type 'application/x-gzip' length 3096278 bytes (3.0 MB)
==================================================
downloaded 3.0 MB

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/car_3.1-3.tar.gz'
Content type 'application/x-gzip' length 1512749 bytes (1.4 MB)
==================================================
downloaded 1.4 MB

* installing *binary* package ‘rbibutils’ ...
* DONE (rbibutils)
* installing *binary* package ‘cowplot’ ...
* DONE (cowplot)
* installing *binary* package ‘Deriv’ ...
* DONE (Deriv)
* installing *binary* package ‘microbenchmark’ ...
* DONE (microbenchmark)
* installing *binary* package ‘numDeriv’ ...
* DONE (numDeriv)
* installing *binary* package ‘SparseM’ ...
* DONE (SparseM)
* installing *binary* package ‘MatrixModels’ ...
* DONE (MatrixModels)
* installing *binary* package ‘nloptr’ ...
* DONE (nloptr)
* installing *binary* package ‘Rcpp’ ...
* DONE (Rcpp)
* installing *binary* package ‘carData’ ...
* DONE (carData)
* installing *binary* package ‘abind’ ...
* DONE (abind)
* installing *binary* package ‘Formula’ ...
* DONE (Formula)
* installing *binary* package ‘maps’ ...
* DONE (maps)
* installing *binary* package ‘Rdpack’ ...
* DONE (Rdpack)
* installing *binary* package ‘doBy’ ...
* DONE (doBy)
* installing *binary* package ‘minqa’ ...
* DONE (minqa)
* installing *binary* package ‘RcppEigen’ ...
* DONE (RcppEigen)
* installing *binary* package ‘quantreg’ ...
* DONE (quantreg)
* installing *binary* package ‘reformulas’ ...
* DONE (reformulas)
* installing *binary* package ‘lme4’ ...
* DONE (lme4)
* installing *binary* package ‘pbkrtest’ ...
* DONE (pbkrtest)
* installing *binary* package ‘car’ ...
* DONE (car)

The downloaded source packages are in
    ‘/tmp/RtmpEDqZ4I/downloaded_packages’

# read the CSV with headers

regression1<-read.csv("/cloud/project/incidents (5).csv", header=T,sep =",")

View(regression1)

summary(regression1)

     area               zone            population          incidents     
 Length:16          Length:16          Length:16          Min.   : 103.0  
 Class :character   Class :character   Class :character   1st Qu.: 277.8  
 Mode  :character   Mode  :character   Mode  :character   Median : 654.0  
                                                          Mean   : 695.2  
                                                          3rd Qu.: 853.0  
                                                          Max.   :2072.0  

str(regression1)

'data.frame':   16 obs. of  4 variables:
 $ area      : chr  "Boulder" "California-lexington" "Huntsville" "Seattle" ...
 $ zone      : chr  "west" "east" "east" "west" ...
 $ population: chr  "107,353" "326,534" "444,752" "750,000" ...
 $ incidents : int  605 103 161 1703 1003 527 721 704 105 403 ...

regression1$population <- as.numeric(gsub(",","",regression1$population))
regression1$population

 [1]  107353  326534  444752  750000   64403 2744878 1600000 2333000 1572816  712091 6900000 2700000
[13] 4900000 4200000 5200000 7100000

str(regression1$population)

 num [1:16] 107353 326534 444752 750000 64403 ...

regression2<-regression1[,-1]#new data frame with the deletion of column 1 

head(regression2)

reg.fit1<-lm(regression1$incidents ~ regression1$population)

summary(reg.fit1)

Call:
lm(formula = regression1$incidents ~ regression1$population)

Residuals:
   Min     1Q Median     3Q    Max 
-684.5 -363.5 -156.2  133.9 1164.7 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)  
(Intercept)            4.749e+02  2.018e+02   2.353   0.0337 *
regression1$population 8.462e-05  5.804e-05   1.458   0.1669  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 534.9 on 14 degrees of freedom
Multiple R-squared:  0.1318,    Adjusted R-squared:  0.0698 
F-statistic: 2.126 on 1 and 14 DF,  p-value: 0.1669

#This output is from fitting a linear regression model where the dependent variable (regression1$incidents) is predicted based on the independent variable (regression1$population). lm(formula = regression1$incidents ~ regression1$population)  shows that the model is trying to predict incidents (the dependent variable) using population as the predictor variable.Residuals are the differences between the observed values and the predicted values. A quick look at these statistics shows that the residuals range from -684.5 to 1164.7.
#The median residual is -156.2, meaning that the model’s predictions tend to slightly under-predict the observed values. Intercept (Estimate = 474.9): When the population is 0, the model predicts about 474.9 incidents. While this may not be realistic in practical terms (since a population of 0 is unlikely), it’s still the baseline value when the population is extremely small.

#Population Coefficient (Estimate = 8.462e-05): For every 1-unit increase in population, the model predicts an increase of approximately 0.0000846 incidents.

#However, the p-value (0.1669) for this coefficient is greater than 0.05, which means this predictor (population) is not statistically significant at the 5% significance level. In other words, there’s not enough evidence to conclude that population has a meaningful effect on the number of incidents.
#The residual standard error of 534.9 is the typical distance that the data points are from the regression line, on average.
#14 degrees of freedom suggests the model was fit with a small sample size (likely 16 observations: 15 data points and 1 degree of freedom for the regression model).

#Multiple R-squared (0.1318) indicates that approximately 13.18% of the variance in the number of incidents is explained by the population variable. This is quite low, meaning the population is not a strong predictor of the number of incidents.
#Adjusted R-squared (0.0698) adjusts for the number of predictors and sample size. This also indicates that the model doesn’t do a great job explaining the variance.
##The p-value (0.1669) is above 0.05, indicating that the model is not statistically significant at the 5% significance level. In simpler terms, the relationship between population and incidents is not statistically meaningful in this model.

reg.fit2<-lm(incidents ~ zone+population, data = regression1)

summary(reg.fit2)

Call:
lm(formula = incidents ~ zone + population, data = regression1)

Residuals:
    Min      1Q  Median      3Q     Max 
-537.21 -273.14  -57.89  188.17  766.03 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) 1.612e+02  1.675e+02   0.962  0.35363   
zonewest    7.266e+02  1.938e+02   3.749  0.00243 **
population  6.557e-05  4.206e-05   1.559  0.14300   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 384.8 on 13 degrees of freedom
Multiple R-squared:  0.5828,    Adjusted R-squared:  0.5186 
F-statistic: 9.081 on 2 and 13 DF,  p-value: 0.003404

#Zone is statistically significant at the 5% significance level, but population is not.
#The Adjusted R-squared of the model is 0.5186, indicating that about 51.86% of the variance in the number of incidents is explained by the model.

regression1$zone <- ifelse(regression1$zone == "west", 1, 0)#Please explain the syntax and the output

View(regression1)

str(regression1)

'data.frame':   16 obs. of  4 variables:
 $ area      : chr  "Boulder" "California-lexington" "Huntsville" "Seattle" ...
 $ zone      : num  1 0 0 1 1 0 1 1 0 0 ...
 $ population: num  107353 326534 444752 750000 64403 ...
 $ incidents : int  605 103 161 1703 1003 527 721 704 105 403 ...

regression1$zone<-as.integer((regression1$zone),replace=TRUE) was not necessary

Error: unexpected symbol in "regression1$zone<-as.integer((regression1$zone),replace=TRUE) was"

interaction<-regression1$zone*regression1$population#Explain the syntax

reg.fit3<-lm(regression1$incidents~interaction+regression1$population+regression1$zone)

summary(reg.fit3)

Call:
lm(formula = regression1$incidents ~ interaction + regression1$population + 
    regression1$zone)

Residuals:
    Min      1Q  Median      3Q     Max 
-540.91 -270.93  -59.56  187.99  767.99 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)  
(Intercept)            1.659e+02  2.313e+02   0.717   0.4869  
interaction            2.974e-06  9.469e-05   0.031   0.9755  
regression1$population 6.352e-05  7.868e-05   0.807   0.4352  
regression1$zone       7.192e+02  3.108e+02   2.314   0.0392 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 400.5 on 12 degrees of freedom
Multiple R-squared:  0.5829,    Adjusted R-squared:  0.4786 
F-statistic: 5.589 on 3 and 12 DF,  p-value: 0.01237

#1. Is Population significant at a 5% significance level?
Population has a p-value of 0.4352 in this model. Since this is greater than 0.05, population is not statistically significant at the 5% significance level.
#2. Is Zone significant at a 5% significance level?
#Zone has a p-value of 0.0392 in this model. Since this is less than 0.05, zone is statistically significant at the 5% significance level. This indicates that the zone has a significant effect on the number of incidents.
#3. Is the interaction term significant at a 5% significance level?
#The interaction term has a p-value of 0.9755. Since this is much greater than 0.05, the interaction term is not statistically significant at the 5% significance level. This suggests that the interaction between the variables has no significant effect on the number of incidents.
#4. What is the Adjusted R-squared of the model?
#The Adjusted R-squared of the model is 0.4786. This means that approximately 47.86% of the variance in the number of incidents is explained by the predictors (population, zone, and interaction).

reg.fit4<-lm(regression1$incidents~interaction)

summary(reg.fit4)

Call:
lm(formula = regression1$incidents ~ interaction)

Residuals:
    Min      1Q  Median      3Q     Max 
-650.28 -301.09  -83.71  123.23 1103.76 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) 4.951e+02  1.320e+02   3.751  0.00215 **
interaction 1.389e-04  4.737e-05   2.932  0.01093 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 451.9 on 14 degrees of freedom
Multiple R-squared:  0.3804,    Adjusted R-squared:  0.3361 
F-statistic: 8.595 on 1 and 14 DF,  p-value: 0.01093

#Decision:
#Model 1 (incidents ~ zone + population) appears to be the most reasonable choice for making predictions despite the population variable not being statistically significant. This model explains more of the variance (47.86%) compared to Model 3, which only explains 33.61%.

#Model 3 only includes the interaction term, which is statistically significant, but the model’s lower Adjusted R-squared (0.3361) and higher residual errors suggest it is less effective at explaining the variance in incidents.

#Model 1 provides more context with two predictors (zone and population) and has a better overall fit. Even though population is not statistically significant, zone is significant, which likely makes the model more reliable than Model 3. For prediction purposes, Model 1 (incidents ~ zone + population) is preferred because it explains more of the variance in the data and has a better overall fit, despite the lack of significance of population. Model 3, with only the interaction term, has lower explanatory power and a higher error rate.