Rationale

According to the cultivation theory, the amount of media that we consume changes people’s understandings on the world and how it works. The more we see a specific trend or theme on television or in any media, the more we believe it to be true.

For example, many television shows relate to criminal justice, medicine, and emergency services workers, mainly because these are viewed as high-stakes, fast-paced environments which are perfect for television. These employees encounter a problem, and have to act quickly to solve it, which may skew viewers beliefs about what regularly goes on in these fields.

Hypothesis

Participants who watch more television per week (“video”) will have a higher estimate of the percentage of the U.S. population employed in law enforcement, medicine, or emergency services (“pct”).

Variables & Method

A total of 400 participants tracked their weekly hours spent watching television content, recruited from a random sample of U.S. adults. Monitoring devices connected to participants’ television(s) and streaming devices automatically recorded the average number of hours viewed per week, over a period of six-months.

Following the six-month period of data collection, participants completed a questionnaire estimating the percentage of the U.S. population employed full-time in law enforcement, criminal justice, medicine, or emergency services. Each category was counted separately, and the percentages were summed.

A linear regression analysis was conducted to determine whether weekly television consumption and participants’ occupation estimates were statistically significant.

Results & Discussion

The graphs below show the distribution of the dependent variable (“pct”) and the independent variable (“video”). The regression plot, outliers table, regression analysis, and model fit statistics are also shown.

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

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

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

Row #	Leverage
Leverage estimates for 10 largest outliers
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

Term	Estimate	Std. Error	t	p-value
Regression Analysis Results
Coefficient Estimates
(Intercept)	23.2076	2.2026	10.5363	0.0000
IV	0.8440	0.0630	13.4056	0.0000

R-squared	Adj. R-squared	F-statistic	df (model)	df (residual)	Residual Std. Error
Model Fit Statistics
Overall Regression Performance
0.3111	0.3093	179.7107	1.0000	398.0000	9.7373

The results supported the hypothesis. While not necessarily seeing a strong correlation (R-squared = 31%), there is a trend that can be seen in the scatterplot. The scatterplot shows a positive, reasonably linear relationship between hours spent watching television and estimated percentage of U.S. law enforcement, medical, and emergency services workers. The regression analysis found the association to be statistically significant.

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$pct
mydata$IV <- mydata$video


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

Cultivation Analysis

Bradley Krpejs

2026-09-21

Rationale

Hypothesis

Variables & Method

Results & Discussion

Code