Lab 1: Intro to R

library(tidyverse)
library(openintro)
library(statsr)
library(broom)
library(ggplot2)

Exercise 1

There are 1458 rows and 123 columns, each row represents a single country in a specific year.

data(hfi)
dim(hfi)

## [1] 1458  123

Exercise 2

hfi_2016 <- hfi %>%
  filter(year == 2016) %>%
  select(pf_score, pf_expression_control, hf_score, region, countries, ISO_code)
head(hfi_2016)

## # A tibble: 6 × 6
##   pf_score pf_expression_control hf_score region              countries ISO_code
##      <dbl>                 <dbl>    <dbl> <chr>               <chr>     <chr>   
## 1     7.60                  5.25     7.57 Eastern Europe      Albania   ALB     
## 2     5.28                  4        5.14 Middle East & Nort… Algeria   DZA     
## 3     6.11                  2.5      5.64 Sub-Saharan Africa  Angola    AGO     
## 4     8.10                  5.5      6.47 Latin America & th… Argentina ARG     
## 5     6.91                  4.25     7.24 Caucasus & Central… Armenia   ARM     
## 6     9.18                  7.75     8.58 Oceania             Australia AUS

Exercise 3

The form of the relationship is linear, the direction is positive, and the strength is strong, with a correlation of 0.845. There are no significant unusual observations, and the data points are fairly tightly clustered along the trend line, supporting the use of a linear model.

ggplot(hfi_2016, aes(x = pf_expression_control, y = pf_score)) +
  geom_point() +
  labs(title = "pf_score vs pf_expression_control")

hfi_2016 %>%
  summarise(correlation = cor(pf_expression_control, pf_score, use = "complete.obs"))

## # A tibble: 1 × 1
##   correlation
##         <dbl>
## 1       0.845

Exercise 4

The closer the fitted line is to the actual data points, the lower the sum of squared residuals (SSR).

plot_ss(x = pf_expression_control, y = pf_score, data = hfi_2016, showSquares = TRUE)

## Click two points to make a line.                                
## Call:
## lm(formula = y ~ x, data = pts)
## 
## Coefficients:
## (Intercept)            x  
##      4.2838       0.5418  
## 
## Sum of Squares:  102.213

Exercise 5

The regression model shows that for each 1-unit increase in expression control, the personal freedom score increases by about 0.542 points. This indicates a strong, positive linear relationship between the two variables.

m1 <- lm(pf_score ~ pf_expression_control, data = hfi_2016)
tidy(m1)

## # A tibble: 2 × 5
##   term                  estimate std.error statistic  p.value
##   <chr>                    <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)              4.28     0.149       28.8 4.23e-65
## 2 pf_expression_control    0.542    0.0271      20.0 2.31e-45

glance(m1)

## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.714         0.712 0.799      400. 2.31e-45     1  -193.  391.  400.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Exercise 6

The predicted personal freedom score for a country with a pf_expression_control of 3 is approximately 5.91. If the actual score is, for example, 5.20, then the residual is −0.71, indicating the model overestimated the true value. This helps us assess how accurate the model is for specific observations.

ggplot(hfi_2016, aes(x = pf_expression_control, y = pf_score)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

predicted <- predict(m1, newdata = data.frame(pf_expression_control = 3))
actual <- hfi_2016$pf_score[hfi_2016$pf_expression_control == 3]
residual <- actual - predicted

predict(m1, newdata = data.frame(pf_expression_control = 3))

##        1 
## 5.909351

print(actual)

## [1] 5.465632

print(residual)

##          1 
## -0.4437186

Exercise 7

Using the regression line, the predicted personal freedom score for a country with a pf_expression_control score of 3 is approximately 5.91. If the actual score is 5.20, the model overestimates the value by 0.71. This results in a residual of −0.71, which indicates that the model’s prediction was too high.

m1_aug <- augment(m1)

ggplot(m1_aug, aes(x = .fitted, y = .resid)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  xlab("Fitted values") +
  ylab("Residuals")

Exercise 8

There is no apparent skew or strong pattern in the histogram of residuals, and the shape appears roughly symmetric and bell-shaped. This suggests that the residuals are approximately normally distributed, which supports the validity of the linearity assumption in the regression model.

ggplot(data = m1_aug, aes(x = .resid)) +
  geom_histogram(binwidth = 0.25, fill = "skyblue") +
  labs(title = "Histogram of Residuals", x = "Residuals", y = "Count")

Exercise 9

The histogram looks roughly symmetric and bell-shaped, so the nearly normal residuals condition is satisfied.

Exercise 10

Based on the residuals vs. fitted values plot, there is no clear pattern or funnel shape, and the spread of residuals appears fairly consistent across all fitted values. This suggests that the constant variability condition is not violated, and the linear regression model satisfies the assumption of homoscedasticity.

m2 <- lm(pf_score ~ hf_score, data = hfi_2016)
glance(m2)$r.squared

## [1] 0.8978487