Rationale

A regression analysis was used to examine whether there is a relationship between the estimated percentage of the U.S. population employed in law enforcement/criminal justice, medicine, or emergency response services and the amount of television content participants watch each week. Specifically, the analysis examines whether changes in pct are associated with changes in video.

In this regression analysis, pct is used as the independent variable to predict video, the dependent variable. The regression model allows the direction and strength of the relationship between the two quantitative variables to be measured and determines whether pct is statistically significant predictor of weekly television viewing.

Hypothesis

The estimated percentage of the U.S. population employed in law enforcement/criminal justice, medicine, or emergency response services is linearly related to average weekly television viewing hours.

Variables & method

The dependent variable in the analysis was video, which measured the average number of hours per week each participant spent watching television content during the six-month monitoring period. The independent variable was pct, which measured each participant’s estimated percentage of the U.S. population employed in law enforcement/criminal justice, medicine, or emergency response services.

A sample of 400 participants was used in the analysis. Television viewing was measured using monitoring devices that recorded the amount of time spent watching television content. At the end of the six-month period, participants completed a questionnaire that included their estimates of the employment percentages.

A linear regression analysis was used to test for a statistically significant relationship between participants’ estimated employment percentage and their weekly television-viewing hours. A scatterplot and histograms were also used to examine the distributions and relationship between the two variables.

Results & discussion

The scatterplot shows a positive relationship between pct and video. This means that as the estimated percentage of people employed in law enforcement/criminal justice, medicine, or emergency response services increases, weekly television viewing hours also tend to increase.

The regression equation was:

video=14.9525+ 0.3686 (pct)

The slope of 0.3686 shows that for every 1% increase in pct, weekly television viewing increases by about 0.37 hours.

The regression was statistically significant, with a p-value less than 0.001. The R squared value was 0.3111, meaning that about 31.1% of the variation in weekly television viewing is related to pct.

Overall, the results show a significant positive linear relationship between the two variables, supporting the hypothesis.

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

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

Leverage estimates for 10 largest outliers
Row # Leverage
359 0.0275
346 0.0236
108 0.0212
198 0.0168
236 0.0168
191 0.0158
388 0.0158
39 0.0130
333 0.0130
392 0.0130
Regression Analysis Results
Coefficient Estimates
Term Estimate Std. Error t p-value
(Intercept) 14.9525 1.4656 10.2026 0.0000
IV 0.3686 0.0275 13.4056 0.0000
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 6.4347

Code

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


##################################################
# 2. Read in the dataset
##################################################
# Replace "YOURFILENAME.csv" with the actual filename
mydata <- read.csv("Cultivation.csv")


# ################################################
# # (Optional) 2b. Remove specific cases by row number
# ################################################
# # Example: remove rows 10 and 25
# rows_to_remove <- c(10, 25) # Edit and uncomment this line
# mydata <- mydata[-rows_to_remove, ] # Uncomment this line


##################################################
# 3. Define dependent variable (DV) and independent variable (IV)
##################################################
# Replace YOURDVNAME and YOURIVNAME with actual column names
mydata$DV <- mydata$video
mydata$IV <- mydata$pct


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

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


##################################################
# 5. Fit and summarize initial regression model
##################################################
# Suppress scientific notation
options(scipen = 999)

# Fit model
myreg <- lm(DV ~ IV, data = mydata)

# Model summary
summary(myreg)


##################################################
# 6. Visualize regression and check for bivariate outliers
##################################################
# Create scatterplot with regression line as a ggplot object
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()


##################################################
# 7. Check for potential outliers (high leverage points)
##################################################
# Calculate leverage values
hat_vals <- hatvalues(myreg)

# Rule of thumb: leverage > 2 * (number of predictors + 1) / n may be influential
threshold <- 2 * (length(coef(myreg)) / nrow(mydata))

# Create table showing 10 largest leverage values
outliers <- data.frame(
  Obs = 1:nrow(mydata),
  Leverage = hatvalues(myreg)
) %>%
  arrange(desc(Leverage)) %>%
  slice_head(n = 10)

# Format as a gt table
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
  )


##################################################
# 8. Create nicely formatted regression results tables
##################################################
# --- Coefficient-level results ---
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
  )


# --- Model fit statistics ---
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
  )


##################################################
# 9. Final print of key graphics and tables
##################################################
DVGraph
IVGraph
RegressionPlot
outliers_table
reg_table
fit_table