##################################################
# 1. Install and load required packages
##################################################
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("gt")) install.packages("gt")
if (!require("gtExtras")) install.packages("gtExtras")

library(tidyverse)
library(gt)
library(gtExtras)

Introduction

Cultivation theory suggests that heavy television viewers come to see the real world the way television portrays it. Television programming overrepresents people who work in law enforcement, medicine, and emergency response, so heavy viewers may overestimate how many Americans hold those jobs.

Hypothesis: Average weekly hours of television viewing will positively predict estimates of the percentage of Americans employed in law enforcement/criminal justice, medicine, or emergency response.

Method

The study used 400 volunteer participants recruited from a random sample of U.S. adults. For six months, monitoring devices recorded the hours each participant personally spent watching television content, averaged into weekly viewing hours (video). At the end of the monitoring period, participants estimated the percentage of the U.S. population employed full time in each of the three occupational categories, and those estimates were summed (pct). A bivariate regression tested whether video predicted pct.

# Read the data from the web
FetchedData <- read.csv("https://github.com/drkblake/Data/raw/refs/heads/main/Cultivation.csv")
# Save the data on your computer
write.csv(FetchedData, "Cultivation.csv", row.names=FALSE)
# remove the data from the environment
rm (FetchedData)

##################################################
# 2. Read in the dataset
##################################################
mydata <- read.csv("Cultivation.csv")

##################################################
# 3. Define dependent variable (DV) and independent variable (IV)
##################################################
mydata$DV <- mydata$pct
mydata$IV <- mydata$video

Results

Distributions

##################################################
# 4. Explore distributions of DV and IV
##################################################
DVGraph <- ggplot(mydata, aes(x = DV)) + 
  geom_histogram(color = "black", fill = "#1f78b4")

IVGraph <- ggplot(mydata, aes(x = IV)) + 
  geom_histogram(color = "black", fill = "#1f78b4")

DVGraph

IVGraph

Regression

##################################################
# 5. Fit and summarize initial regression model
##################################################
options(scipen = 999)
myreg <- lm(DV ~ IV, data = mydata)
summary(myreg)
## 
## Call:
## lm(formula = DV ~ IV, data = mydata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -28.5917  -6.1607   0.6423   6.7433  23.8483 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 23.20761    2.20264   10.54 <0.0000000000000002 ***
## IV           0.84400    0.06296   13.41 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.737 on 398 degrees of freedom
## Multiple R-squared:  0.3111, Adjusted R-squared:  0.3093 
## F-statistic: 179.7 on 1 and 398 DF,  p-value: < 0.00000000000000022
##################################################
# 6. Visualize regression and check for bivariate outliers
##################################################
RegressionPlot <- ggplot(mydata, aes(x = IV, y = DV)) +
  geom_point(color = "#1f78b4") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(
    title = "Scatterplot of DV vs IV with Regression Line",
    x = "Independent Variable (IV)",
    y = "Dependent Variable (DV)"
  ) +
  theme_minimal()

RegressionPlot

Outlier check

##################################################
# 7. Check for potential outliers (high leverage points)
##################################################
hat_vals <- hatvalues(myreg)
threshold <- 2 * (length(coef(myreg)) / nrow(mydata))

outliers <- data.frame(
  Obs = 1:nrow(mydata),
  Leverage = hatvalues(myreg)
) %>%
  arrange(desc(Leverage)) %>%
  slice_head(n = 10)

outliers_table <- outliers %>%
  gt() %>%
  tab_header(title = "Leverage estimates for 10 largest outliers") %>%
  cols_label(Obs = "Row #", Leverage = "Leverage") %>%
  fmt_number(columns = Leverage, decimals = 4)

outliers_table
Leverage estimates for 10 largest outliers
Row # Leverage
164 0.0305
360 0.0207
359 0.0194
371 0.0174
72 0.0162
201 0.0159
265 0.0159
392 0.0148
44 0.0144
97 0.0144
 
threshold
## [1] 0.01

Regression tables

##################################################
# 8. Create nicely formatted regression results tables
##################################################
reg_results <- as.data.frame(coef(summary(myreg))) %>%
  tibble::rownames_to_column("Term") %>%
  rename(
    Estimate = Estimate,
    `Std. Error` = `Std. Error`,
    t = `t value`,
    `p-value` = `Pr(>|t|)`
  )

reg_table <- reg_results %>%
  gt() %>%
  tab_header(title = "Regression Analysis Results",
             subtitle = "Coefficient Estimates") %>%
  fmt_number(columns = c(Estimate, `Std. Error`, t, `p-value`), decimals = 4)

reg_summary <- summary(myreg)

fit_stats <- tibble::tibble(
  `R-squared` = reg_summary$r.squared,
  `Adj. R-squared` = reg_summary$adj.r.squared,
  `F-statistic` = reg_summary$fstatistic[1],
  `df (model)` = reg_summary$fstatistic[2],
  `df (residual)` = reg_summary$fstatistic[3],
  `Residual Std. Error` = reg_summary$sigma
)

fit_table <- fit_stats %>%
  gt() %>%
  tab_header(title = "Model Fit Statistics",
             subtitle = "Overall Regression Performance") %>%
  fmt_number(columns = everything(), decimals = 4)

reg_table
Regression Analysis Results
Coefficient Estimates
Term Estimate Std. Error t p-value
(Intercept) 23.2076 2.2026 10.5363 0.0000
IV 0.8440 0.0630 13.4056 0.0000
 
fit_table
Model Fit Statistics
Overall Regression Performance
R-squared Adj. R-squared F-statistic df (model) df (residual) Residual Std. Error
0.3111 0.3093 179.7107 1.0000 398.0000 9.7373