–Welcome to Lecture#8,Chapter V.

Let us start by finding our working directoty.

getwd() # My working directory is the documents
## [1] "C:/Users/jonat/OneDrive/Documents"
# 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, repos = "https://cloud.r-project.org/")  
}
## Installing package into 'C:/Users/jonat/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)
## Warning: 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
# a different approach to installing the packages

Let us now make sure that the file zeroaccess is uploaded into your new project. Then let us define the variable za.

# read the CSV with headers
za <- read.csv("zeroaccess.csv", header=T,sep ="," )
#View(za)
# Always comment out the view function because when publishing this will cause an error. If you do not comment out the view function then you will have a popup error.

What is the output we have just obtained after running the chunk above. If you comment out the view function you will not get the pop up The output we received is the file being read and a view window showing the data

# Load ggplot2 to create graphics
install.packages('ggplot2', repos = "https://cloud.r-project.org/")
## Installing package into 'C:/Users/jonat/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)
## package 'ggplot2' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\jonat\AppData\Local\Temp\Rtmp0qxp5S\downloaded_packages
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.2
# create a ggplot instance with zeroaccess data
gg <- ggplot(data=za, aes(x=long, y=lat)) # call the data, and call the x and y 

# add the points, set transparency to 1/40th 
gg <- gg + geom_point(size=1, color="#000099", alpha=1/40) 

# add axes labels
gg <- gg + xlab("Longitude") + ylab("Latitude") # this adds labels to the plot

# simplify the theme for aesthetics
gg <- gg + theme_bw() 

# this may take a while, over 800,000 points plotted
print(gg)

# This R code utilizes the ggplot2 package to create a scatter plot of geographical data
#install.packages("mapproj", repos = "http://cran.us.r-project.org")
#library(mapproj)
# requires package : ggplot2
# requires object: za (5-1)
# the "maps" and "mapproj" packages are used by ggplot
# load map data of the world
world <- map_data("world")
#Remove Antarctica
world <- subset(world, world$region!="Antarctica")
# load world data into ggplot object
gg <- ggplot(data=world, aes(x=long, y=lat))
# trace along the lat/long coords by group (countries)
gg <- gg + geom_path(aes(group=group), colour="gray70")
# now project using the mercator projection
# try different projections with ?mapproject
gg <- gg + coord_map("mercator", xlim=c(-200, 200))
# load up the ZeroAccess points, overiding the default data set
gg <- gg + geom_point(data=za, aes(long, lat), 
                      colour="#000099", alpha=1/40, size=1)
# remove text, axes ticks, grid lines and do gray border on white
gg <- gg + theme(text=element_blank(), 
                 axis.ticks=element_blank(),
                 panel.grid=element_blank(),
                 panel.background=element_rect(color="gray50",
                                               fill="white"))
print(gg)

# This code creates a plot showing a map of the world with country boundaries traced in light gray and overlays it with data points

Interpret the output above.

#Let us proceed to define the file county.dataset found on blackboard.

# read up census data per county
county.data <- read.csv("countydataset.csv", header=T,sep = ",")
# Header = T reads the first row as the columns names
#View(county.data)
# this allows you to view the data in full, always comment out
set.seed(1)
# generate 200 random numbers around 10
input <- rnorm(200, mean=10)
summary(input)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   7.785   9.386   9.951  10.036  10.613  12.402
# This code generates a set of 200 random numbers drawn from a normal distribution with a mean of 10 and a default standard deviation of 1

Interpret the synatx and the output above.It suggests that the numbers are roughly centered around the mean of 10, with a spread from approximately 7.785 to 12.402.

# requires objects: input (5-16)
# generate output around a mean of 2 x input
output <- rnorm(200, mean=input*2)
# put into data frame to plot it
our.data <- data.frame(input, output)
gg <- ggplot(our.data, aes(input, output))
gg <- gg + geom_point()
gg <- gg + geom_smooth(method = "lm", se=F, color="red")
gg <- gg + theme_bw()
print(gg)
## `geom_smooth()` using formula = 'y ~ x'

# this code generates a scatter plot of input v output and adds a line for linear regression

Make comments about both the syntax and the output of the task executed above. Is it possible to customize this graph to make it more explanatory? If so explain and run the code. The graph and output above is very well distributed and easily understandable in its current state

Let us now run output vs input.

model <- lm(output ~ input)
summary(model)
## 
## Call:
## lm(formula = output ~ input)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.93275 -0.54273 -0.02523  0.66833  2.58615 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.27224    0.77896   0.349    0.727    
## input        1.97692    0.07729  25.577   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.013 on 198 degrees of freedom
## Multiple R-squared:  0.7677, Adjusted R-squared:  0.7665 
## F-statistic: 654.2 on 1 and 198 DF,  p-value: < 2.2e-16
# this code fits a linear regression model of output vs input variables and prints a summary of the model

Is the input significant at a 5% significance level?1%?

Yes because the P-value is less than 0.05

#If it was 1% significance then it would be a P-value of less than 0.01

Let us now proceed to get the confidence interval for this model.

confint(model) # provides confidence level
##                 2.5 %   97.5 %
## (Intercept) -1.263895 1.808368
## input        1.824502 2.129343
# if it does not go through 0, it is significant

Using the county.data dataset, let us run a model of Infections vs ufo2010(Aliens Visits according to the UFO).

summary(lm(county.data$Infections ~ county.data$ufo2010, data= county.data))
## 
## Call:
## lm(formula = county.data$Infections ~ county.data$ufo2010, data = county.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -407.50  -70.35   -5.78   60.05 2736.39 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         124.77740    2.38569  52.302  < 2e-16 ***
## county.data$ufo2010   0.56899    0.07998   7.114  1.4e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 127.5 on 3070 degrees of freedom
## Multiple R-squared:  0.01622,    Adjusted R-squared:  0.0159 
## F-statistic: 50.61 on 1 and 3070 DF,  p-value: 1.398e-12
# This R code is fitting a linear regression model to the relationship between the number of infections and a predictor variable using the lm function
# In the different us regions,there is not enough evidence that more reports of UFO's affects the number of infections in computer stations

Given the output printed above, is ufo2010 significant at a 5% significance level?1%? Yes, it is significant at a 5% significance level and at 1% because the P-value is less than 0.05 and less than 0.01

Let us now run a model of Infections vs every single quantitative variable that is included in the dataset.

#View(county.data)
# this views the dataset and always comment out the view function
summary(lm(Infections ~ pop + income + ipaddr + ufo2010, 
           data=county.data))
## 
## Call:
## lm(formula = Infections ~ pop + income + ipaddr + ufo2010, data = county.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -342.94  -70.34   -5.40   58.85 2745.36 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.000e+02  9.426e+00  10.610  < 2e-16 ***
## pop         -1.615e-05  1.543e-05  -1.047  0.29510    
## income       5.670e-04  2.066e-04   2.745  0.00609 ** 
## ipaddr      -9.865e-08  5.148e-07  -0.192  0.84803    
## ufo2010      6.804e-01  1.691e-01   4.024 5.85e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 127.3 on 3067 degrees of freedom
## Multiple R-squared:  0.01879,    Adjusted R-squared:  0.01751 
## F-statistic: 14.69 on 4 and 3067 DF,  p-value: 6.915e-12
# The R code you provided fits a multiple linear regression model to predict the number of infections based on several variables
# by looking at the R squared, it is not significant

Interpret the output. How would you proceed from now on in this handout given the results obtained above. The population and ipaddr are not significant because they are above 0.05 at a 5% significance level

install.packages("car", repos = "http://cran.us.r-project.org")
## Installing package into 'C:/Users/jonat/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)
## package 'car' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\jonat\AppData\Local\Temp\Rtmp0qxp5S\downloaded_packages
library(car) # for the vif() function
## Warning: package 'car' was built under R version 4.3.2
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.2

Let us just explore the variance inflation factor(VIF) of the model to see if there is a chance of high correlation between my predictors. I remind you that a strong correlation between two of my predictors will likely end up in heteroskedasticity, and therefore our model would not be accuarate. That is high correlation between values

model <- lm(Infections ~ pop + income + ipaddr + ufo2010, 
            data=county.data)
sqrt(vif(model))
##      pop   income   ipaddr  ufo2010 
## 2.165458 1.038467 1.046051 2.115512
# The code provided calculates the square root of the variance inflation factor (VIF) for each predictor variable in the multiple linear regression model
# Values greater than 5 are typically considered indicative of significant multicollinearity.

Let us see if population affects the number of infections in this dataset. Write the null and alternative hypothesis you would use to test this relationship.

summary(lm(Infections ~ pop, data=county.data))
## 
## Call:
## lm(formula = Infections ~ pop, data = county.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -365.26  -70.35   -5.70   59.59 2724.00 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.250e+02  2.416e+00  51.755  < 2e-16 ***
## pop         4.228e-05  7.147e-06   5.916 3.67e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 127.8 on 3070 degrees of freedom
## Multiple R-squared:  0.01127,    Adjusted R-squared:  0.01095 
## F-statistic:    35 on 1 and 3070 DF,  p-value: 3.668e-09
# overall goodness of fit of the simple linear regression model.

Interpret the results obtained above. pop is significant because it is less than 0.05

Now let us define the regression of Infections vs pop as pop.lm and predict the number of infections based on the variable population.

# population of 6000000
pop.lm <- lm(Infections ~ pop, data=county.data)
predict(pop.lm, data.frame(pop=6000000), interval="confidence")
##        fit      lwr      upr
## 1 378.7109 295.9209 461.5009
# lwr is the confidence level and fit is the amount of infections
# the range of values from this output is quite large
# population of 1000000
pop.lm <- lm(Infections ~ pop, data=county.data)
predict(pop.lm, data.frame(pop=1000000), interval="confidence")
##        fit      lwr      upr
## 1 167.3069 153.9224 180.6915
# these values have a mid size range
# population of 2000000
pop.lm <- lm(Infections ~ pop, data=county.data)
predict(pop.lm, data.frame(pop=2000000), interval="confidence")
##        fit      lwr      upr
## 1 209.5877 182.5947 236.5807
income.lm <- lm(Infections ~ income, data=county.data)
predict(income.lm, data.frame(income=10000), interval="confidence")
##        fit      lwr      upr
## 1 100.2013 85.70806 114.6945
# very large range between upper and lower
income.lm <- lm(Infections ~ income, data=county.data)
predict(income.lm, data.frame(income=20000), interval="confidence")
##        fit      lwr      upr
## 1 108.4966 97.66147 119.3318
# the range of this is small but not as well fit at  the other values
income.lm <- lm(Infections ~ income, data=county.data)
predict(income.lm, data.frame(income=40000), interval="confidence")
##        fit      lwr      upr
## 1 125.0872 120.1358 130.0387
# Based on 40,000 in income, it provides a an upper and lower level that has a small range. This shows that this may be the best possible fit for the correlation between the number of infections and income. 
# the range for these values is very small showing that this may be the best interval range
income.lm <- lm(Infections ~ income, data=county.data)
predict(income.lm, data.frame(income=90000), interval="confidence")
##        fit      lwr      upr
## 1 166.5638 148.3582 184.7694
# the range for these is not too large, but it is not the best fit

Interpret the results obtained above.