Handout9_Ch5_Vis_Reg

–Welcome to Lecture#8,Chapter V.

Let us start by finding our working directoty.

getwd()

[1] "C:/Users/antho/OneDrive/Documents/School/4.DataSecurity&Governance"

# 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)  
}

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)
head(za)

The file zeroaccess.csv was in my working directory folder so it uploaded We import the zeroaccess.csv using the read() function, notice the header is true(T) so the first row is considered headers. the values are separated by commas as should be in a comma separated file.

Notice above I used the head() function to view the fields, their data types, and weather they match.

# Load ggplot2 to create graphics
library(ggplot2)


# create a ggplot instance with zeroaccess data
gg <- ggplot(data=za, aes(x=long, y=lat)) 
# 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")
# simplify the theme for aesthetics
gg <- gg + theme_bw() 
# this may take a while, over 800,000 points plotted
print(gg)

Notice ggplot plots the points we provide it from our dataframe, and we can see that points cluster in the shape of what appears to be a map. That is because we are charting geographical cordinates latitude and longitude.

install.packages("mapproj")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/antho/AppData/Local/R/win-library/4.2’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.2/mapproj_1.2.8.zip'
Content type 'application/zip' length 59376 bytes (57 KB)
downloaded 57 KB

package ‘mapproj’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\antho\AppData\Local\Temp\RtmpIJe7Qo\downloaded_packages

We install the package mapproj above and load it to our notebook below.

library(mapproj)

Warning: package ‘mapproj’ was built under R version 4.2.1Loading required package: maps
Warning: package ‘maps’ was built under R version 4.2.1

# 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)

Notice by using the map_data() function from the “mapproj” package we are able to get the outline of the couninents and countries and not just the geographic lat/long points plotted on a scatter plot. Also notice we had to exclude some of the loaded map, antartica specifically which makes sense since we dont have records from there.

#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 = ",")

View(county.data)

After importing the countydataset.csv file from my working directory folder and using the view() function we can see we have 3,072 records 7 columns on infections, ip address, ufo sighting, population and area

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 

Above we test the summary function on 200 random numbers that are meant to have a mean of ~10. Interesting we can ensure that that we get the same random numbers by setting the seed to a specific number in the seed function, this is helpful when trying test thing with the same randomization everytime so that the values dont change inbetween tests.

The summary function outputs the Min, Max, Median and Mean, Not to mention the 1st and 3rd quantile values. We cab see that the mean did stay around 10

# 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)

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.

We use the variable we previously created with random number (input) to create a new variable also with 200 random numbers but this time their mean is twice that of the input. Then we use ggpplot again to graph these two points and a line using the the geo_smooth() function notice the “lm” method which creates the linear regression trend line.

We can add more information to the graph to convey more by adding a Title to it or naming the x and y variables. below i will give an example of renaming the chart.

print(gg + ggtitle("Random Input Vs Output-Mean Multiplied"))

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

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

Yes our p-value is 2.2e-16 which is almost 0 and less than 1% This is good. We also don’t have the worst r-squared at 76%

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

confint(model)

                2.5 %   97.5 %
(Intercept) -1.263895 1.808368
input        1.824502 2.129343

We can se that we have a 95% confidence level between the input coefficients 1.82 - 2.13

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

Given the output printed above, is ufo2010 significant at a 5% significance level?1%?

The p-value is significant since it is far below 1% sadly the r-squared value is very low. So ufo sighings do not predict infections well.

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

View(county.data)

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

Interpret the output. How would you proceed from now on in this handout given the results obtained above.

We can see that the p value for for population and ip address is high so they are least probable for predicting a high accuracy linear model. It would be in our interest to remove them(“pop” & “ipaddr”) from our model.

install.packages("carData")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/antho/AppData/Local/R/win-library/4.2’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.2/carData_3.0-5.zip'
Content type 'application/zip' length 1821225 bytes (1.7 MB)
downloaded 1.7 MB

package ‘carData’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\antho\AppData\Local\Temp\RtmpIJe7Qo\downloaded_packages

We install the package carData above and load the car package to our notebook below.

library(car) # for the vif() function

Warning: package ‘car’ was built under R version 4.2.1Loading required package: carData
Warning: package ‘carData’ was built under R version 4.2.1

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 accurate.

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 

We get the variance inflation factors for our predictors variables.

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

The null hypothesis is we would see that populations changes, while infections wouldn’t The alternative hypothesis we would see is that populations changes, as well as infections

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

Interpret the results obtained above.

We can see from the output above that the p-value is less than 1% so we have statistical significance but since we have such a low r-squared this has nearly no provability of predicting for infections.

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

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

Interpret the results obtained above.

Above we predict the population that will be infected of 6000000, of which the there is 95% confidence the value will fall between 296 and 462. The prediction of our fitted model is 379 which we wont value since the r-squared above showed our model was only accessing 1% of the populations variability.