Handout10

–Welcome to Lecture#8,Chapter V.

Let us start by finding our working directoty.

getwd()

[1] "/Users/gretacapelletti/Downloads"

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

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

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

What is the output we have just obtained after running the chunk above.

The output is that the dataset zeroaccess.csv has been successfully loaded into R as a data frame named za. We can confirm this by running:

View(za)

I can clearly see the dataset in another window in R.

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

install.packages("mapproj")

trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-x86_64/contrib/4.4/mapproj_1.2.11.tgz'
Content type 'application/x-gzip' length 85498 bytes (83 KB)
==================================================
downloaded 83 KB

The downloaded binary packages are in
    /var/folders/g8/k40fgv1n08g8cmq06dchpspw0000gn/T//RtmpwA0CWp/downloaded_packages

library(mapproj)

Loading required package: maps

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

Interpret the output above.

The output is a world map with blue points representing geospatial data from the za dataset. Using the Mercator projection, it shows a high concentration of points in North America, Europe, and parts of Asia, with fewer in Africa and South America. Besides that, overall the distribution of points is uneven. #Let us proceed to define the file county.dataset found on Canvas.

# read up census data per county
county.data <- read.csv("countydataset.csv", header=T,sep = ",")

View(county.data)

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 

Interpret the synatx and the output above.

The code first sets a seed for reproducibility, ensuring that the same random numbers are generated each time the script is run. It then generates 200 random numbers from a normal distribution with a mean of 10 and the default standard deviation of 1. The summary(input) function provides key descriptive statistics for the dataset. The minimum value is 7.785, while the maximum is 12.402, showing the range of the data. The first quartile (9.386) and third quartile (10.613) indicate that the middle 50% of values are concentrated around the mean. The median (9.951) is very close to the mean (10.036), confirming that the distribution is approximately symmetrical 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.

The syntax used in the graph is based on ggplot2 in R, where geom_smooth() is applied with the formula y ~ x, indicating a linear regression model. The black dots represent individual data points, while the red line is the fitted regression line showing the trend between input and output.

To improve the visualization and make it more explanatory, I ensure that the points are more visible by modifying their size and transparency using alpha and size in geom_point(). The red regression line can be made bolder and complemented with a confidence interval (se = TRUE) for better interpretation. Axis labels and a title should be added for clarity.

# Generate output around a mean of 2 x input
output <- rnorm(200, mean = input * 2)

# Put into a data frame
our.data <- data.frame(input, output)

# Create plot with enhancements
gg <- ggplot(our.data, aes(input, output)) +
  geom_point(alpha = 0.6, size = 2, color = "blue") +  # Add transparency and color
  geom_smooth(method = "lm", se = TRUE, color = "red", linewidth = 1.2) +  # Add confidence interval
  labs(title = "Scatter Plot of Input vs. Output",
       x = "Input Values",
       y = "Output Values") +  # Add labels
  theme_minimal()  # Use a cleaner theme

print(gg)

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.9079 -0.6460 -0.0414  0.6721  3.7303 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.82971    0.82296  -1.008    0.315    
input        2.07851    0.08166  25.454   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.07 on 198 degrees of freedom
Multiple R-squared:  0.7659,    Adjusted R-squared:  0.7648 
F-statistic: 647.9 on 1 and 198 DF,  p-value: < 2.2e-16

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

The p-value for input is reported as < 2e-16, which is much smaller than both 0.05 and 0.01. Since the p-value is less than both thresholds, we reject the null hypothesis that the coefficient of input is zero. This means input is highly significant at both the 5% and 1% levels, strongly suggesting that it has a meaningful impact on output

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

confint(model)

                2.5 %    97.5 %
(Intercept) -2.452598 0.7931737
input        1.917481 2.2395382

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 for county.data$ufo2010 is 1.4e-12, which is much smaller than both 0.05 and 0.01. Since the p-value is significantly below both thresholds, we reject the null hypothesis that the coefficient of ufo2010 is zero. This means ufo2010 is highly significant at both the 5% and 1% levels, indicating a statistically significant (but likely weak) relationship between UFO sightings in 2010 and infections.

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. The linear regression model indicates that some predictors, such as income and ufo2010, are statistically significant in explaining the variation in infections, while others like population and ipaddr are not. The low R-squared value (0.01879) suggests that the model does not explain much of the variation in infections, indicating a weak fit. The significant relationship between income and infections suggests that higher income might be associated with an increase in infections, while the positive association with ufo2010 points to a potential link that warrants further exploration.

install.packages("carData")

trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-x86_64/contrib/4.4/carData_3.0-5.tgz'
Content type 'application/x-gzip' length 1827823 bytes (1.7 MB)
==================================================
downloaded 1.7 MB

The downloaded binary packages are in
    /var/folders/g8/k40fgv1n08g8cmq06dchpspw0000gn/T//RtmpKrtZZN/downloaded_packages

# Load the 'car' package (ensure it's installed first: install.packages("car"))
library(car)

# Use the vif() function to check variance inflation factor for multicollinearity
# Replace 'model' with your actual regression model (e.g., lm_model)
vif(model)

     pop   income   ipaddr  ufo2010 
4.689208 1.078413 1.094222 4.475392 

# VIF > 5 (or 10) indicates high multicollinearity, which can affect model accuracy

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.

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 

None of the VIF values are particularly high. Typically, a VIF value greater than 10 is considered problematic, indicating multicollinearity between predictors. Since all the VIF values are below 5, there is no strong multicollinearity in my model, so multicollinearity is unlikely to be a concern.

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

Interpret the results obtained above. The p-value for population is 3.67e-09, which is much smaller than 0.05, so we reject the null hypothesis and conclude that population significantly affects infections. The R-squared value of 0.01127 means population explains only about 1.13% of the variation in infections, suggesting it has a small effect overall.

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. The predicted number of infections for a population of 6 million is 378.71. The 95% confidence interval for this prediction is between 295.92 and 461.50, which means we are 95% confident that the true number of infections for a population of 6 million lies within this range.

Predict the number of infections based on the following values for the variable population: 1000000,2000000.

predict(pop.lm, data.frame(pop=c(1000000, 2000000)), interval="confidence")

       fit      lwr      upr
1 167.3069 153.9224 180.6915
2 209.5877 182.5947 236.5807

These predictions show how the number of infections increases as the population increases.

Recreate this solution building a model that depends on income. Predict the number of infections based on the following values for the variable income: 10000,20000,40000,90000.

income.lm <- lm(Infections ~ income, data=county.data)
predict(income.lm, data.frame(income=c(10000, 20000, 40000, 90000)), interval="confidence")

       fit       lwr      upr
1 100.2013  85.70806 114.6945
2 108.4966  97.66147 119.3318
3 125.0872 120.13584 130.0387
4 166.5638 148.35823 184.7694