Rationale

Cultivation theory says that constant exposure to television can affect viewers’ perceptions of social reality. Television may cause certain occupations, events, or circumstances to appear more normal than they actually are. This analysis shows whether the amount of television participants watch is associated with their estimates of the percentage of the U.S. population employed full-time in either law enforcement/criminal justice, medicine, or emergency response services.

Hypothesis

Participants who watch more television per week will estimate a higher percentage of workers in law enforcement/criminal justice, medicine, or emergency response services in the U.S. population.

Variables & method

A total of 400 volunteer participants were recruited from a random sample of U.S. adults. Participants’ television viewing was monitored for six months, and they completed a questionnaire afterward. The independent variable: video, measured each participant’s average weekly hours of television watched during the six-month period. The dependent variable: pct, measured each participant’s estimated percentage of the U.S. population that is employed in law enforcement/criminal justice, medicine, or emergency response services.

Both variables were continuous measurements. A bivariate linear regression was done to see whether weekly television viewing was associated with participants’ estimates of the percentage of the population employed in these occupations.

Results & discussion

The scatterplot below shows a positive association between weekly television viewing and participants’ estimates of the percentage of the U.S. population that is employed in law enforcement/criminal justice, medicine, or emergency response services. The regression table showed a statistically significant relationship between the dependent and independent variables (p < .001). The fit table shows a 31.1% variation in the dependent variable (R² = .311).  

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

The results supported the hypothesis. Participants who watched more television tended to estimate a higher percentage of the population worked in these occupations. This finding is consistent with cultivation theory, suggesting that more television exposure can be associated with perceptions of social reality.

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

NOTE: If you see what looks like an outlier or two, you can return to this part of the script:

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

Note: The code below will download the (made-up) data used in the analysis.

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