Punishment Study 2.1.2 - Final Analysis
DGK
January 2026
Pre-registered Hypotheses: H1: Punitiveness positively correlated with all correlate measures H2: Harshness cluster > Crime Concerns cluster in predicting punitiveness H3: Most correlate measures positively intercorrelated
PHASE 0: SETUP AND PACKAGE LOADING
# Clear environment
rm(list = ls())
# Required packages - install if needed
required_packages <- c(
"tidyverse", # Data manipulation and visualization
"psych", # Cronbach's alpha, correlations, factor analysis
"corrplot", # Correlation heat maps
"Hmisc", # rcorr for correlation matrices with p-values
"cocor", # Steiger's Z-test for comparing correlations
"lme4", # Mixed-effects models
"lmerTest", # P-values for mixed models
"car", # ANOVA, diagnostics
"effectsize", # Effect size calculations
"ggcorrplot", # Alternative correlation plots
"RColorBrewer", # Color palettes
"knitr", # Tables
"broom", # Tidy model outputs
"scales", # Percentage formatting
"performance", # ICC and model diagnostics
"lavaan" # CFA
)
# Install missing packages
new_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)
# Load packages
lapply(required_packages, library, character.only = TRUE)## [[1]]
## [1] "lubridate" "forcats" "stringr" "dplyr" "purrr" "readr"
## [7] "tidyr" "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "psych" "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [7] "readr" "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [13] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "corrplot" "psych" "lubridate" "forcats" "stringr" "dplyr"
## [7] "purrr" "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [13] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [19] "base"
##
## [[4]]
## [1] "Hmisc" "corrplot" "psych" "lubridate" "forcats" "stringr"
## [7] "dplyr" "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [13] "tidyverse" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"
##
## [[5]]
## [1] "cocor" "Hmisc" "corrplot" "psych" "lubridate" "forcats"
## [7] "stringr" "dplyr" "purrr" "readr" "tidyr" "tibble"
## [13] "ggplot2" "tidyverse" "stats" "graphics" "grDevices" "utils"
## [19] "datasets" "methods" "base"
##
## [[6]]
## [1] "lme4" "Matrix" "cocor" "Hmisc" "corrplot" "psych"
## [7] "lubridate" "forcats" "stringr" "dplyr" "purrr" "readr"
## [13] "tidyr" "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [19] "grDevices" "utils" "datasets" "methods" "base"
##
## [[7]]
## [1] "lmerTest" "lme4" "Matrix" "cocor" "Hmisc" "corrplot"
## [7] "psych" "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [13] "readr" "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [19] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[8]]
## [1] "car" "carData" "lmerTest" "lme4" "Matrix" "cocor"
## [7] "Hmisc" "corrplot" "psych" "lubridate" "forcats" "stringr"
## [13] "dplyr" "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [19] "tidyverse" "stats" "graphics" "grDevices" "utils" "datasets"
## [25] "methods" "base"
##
## [[9]]
## [1] "effectsize" "car" "carData" "lmerTest" "lme4"
## [6] "Matrix" "cocor" "Hmisc" "corrplot" "psych"
## [11] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [16] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
##
## [[10]]
## [1] "ggcorrplot" "effectsize" "car" "carData" "lmerTest"
## [6] "lme4" "Matrix" "cocor" "Hmisc" "corrplot"
## [11] "psych" "lubridate" "forcats" "stringr" "dplyr"
## [16] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [21] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [26] "datasets" "methods" "base"
##
## [[11]]
## [1] "RColorBrewer" "ggcorrplot" "effectsize" "car" "carData"
## [6] "lmerTest" "lme4" "Matrix" "cocor" "Hmisc"
## [11] "corrplot" "psych" "lubridate" "forcats" "stringr"
## [16] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [21] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [26] "utils" "datasets" "methods" "base"
##
## [[12]]
## [1] "knitr" "RColorBrewer" "ggcorrplot" "effectsize" "car"
## [6] "carData" "lmerTest" "lme4" "Matrix" "cocor"
## [11] "Hmisc" "corrplot" "psych" "lubridate" "forcats"
## [16] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [21] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [26] "grDevices" "utils" "datasets" "methods" "base"
##
## [[13]]
## [1] "broom" "knitr" "RColorBrewer" "ggcorrplot" "effectsize"
## [6] "car" "carData" "lmerTest" "lme4" "Matrix"
## [11] "cocor" "Hmisc" "corrplot" "psych" "lubridate"
## [16] "forcats" "stringr" "dplyr" "purrr" "readr"
## [21] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [26] "graphics" "grDevices" "utils" "datasets" "methods"
## [31] "base"
##
## [[14]]
## [1] "scales" "broom" "knitr" "RColorBrewer" "ggcorrplot"
## [6] "effectsize" "car" "carData" "lmerTest" "lme4"
## [11] "Matrix" "cocor" "Hmisc" "corrplot" "psych"
## [16] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [21] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [26] "stats" "graphics" "grDevices" "utils" "datasets"
## [31] "methods" "base"
##
## [[15]]
## [1] "performance" "scales" "broom" "knitr" "RColorBrewer"
## [6] "ggcorrplot" "effectsize" "car" "carData" "lmerTest"
## [11] "lme4" "Matrix" "cocor" "Hmisc" "corrplot"
## [16] "psych" "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [26] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [31] "datasets" "methods" "base"
##
## [[16]]
## [1] "lavaan" "performance" "scales" "broom" "knitr"
## [6] "RColorBrewer" "ggcorrplot" "effectsize" "car" "carData"
## [11] "lmerTest" "lme4" "Matrix" "cocor" "Hmisc"
## [16] "corrplot" "psych" "lubridate" "forcats" "stringr"
## [21] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [26] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [31] "utils" "datasets" "methods" "base"# Set options
options(scipen = 999) # Avoid scientific notation
options(digits = 4) # Default decimal placesPHASE 1: DATA LOADING AND CLEANING
1.1 Load Data
# Read the data file
# UPDATE THIS PATH TO YOUR LOCAL FILE LOCATION
raw_data <- read.csv("/Users/dgkamper/Library/CloudStorage/[email protected]/My Drive/DGK Lab/Collaborations/Dan Simon/Punishment/Data/Data_Final.csv", stringsAsFactors = FALSE)
cat("Raw data loaded:", nrow(raw_data), "rows,", ncol(raw_data), "columns\n\n")## Raw data loaded: 514 rows, 95 columns1.2 Data Filtering
# Track exclusions at each step
n_raw <- nrow(raw_data)
cat("=== DATA FILTERING ===\n\n")## === DATA FILTERING ===cat("Starting N:", n_raw, "\n")## Starting N: 514# Filter 1: Progress = 100 (completed survey)
df <- raw_data %>% filter(Progress == 100)
n_after_progress <- nrow(df)
cat("After Progress = 100:", n_after_progress,
"(excluded", n_raw - n_after_progress, ")\n")## After Progress = 100: 499 (excluded 15 )# Filter 2: Finished = 1
df <- df %>% filter(Finished == 1)
n_after_finished <- nrow(df)
cat("After Finished = 1:", n_after_finished,
"(excluded", n_after_progress - n_after_finished, ")\n")## After Finished = 1: 499 (excluded 0 )# Filter 3: Attention Check 1 = 9 ("Strongly agree")
df <- df %>% filter(AttnCheck1 == 9)
n_after_attn1 <- nrow(df)
cat("After AttnCheck1 = 9:", n_after_attn1,
"(excluded", n_after_finished - n_after_attn1, ")\n")## After AttnCheck1 = 9: 497 (excluded 2 )# Filter 4: Attention Check 2 = 4 (8+7=15, option 4)
df <- df %>% filter(AttnCheck2 == 4)
n_after_attn2 <- nrow(df)
cat("After AttnCheck2 = 4:", n_after_attn2,
"(excluded", n_after_attn1 - n_after_attn2, ")\n")## After AttnCheck2 = 4: 496 (excluded 1 )# Final sample
cat("\n*** FINAL SAMPLE: N =", nrow(df), "***\n")##
## *** FINAL SAMPLE: N = 496 ***cat("Total excluded:", n_raw - nrow(df), "\n\n")## Total excluded: 181.3 Check for Missing Data
# Key variables to check
key_vars <- c("Sentence_1", "punishmore_1", "punishmore_2", "vignette",
"hate_1", "revenge_1", "rwa_1", "sdo_1")
missing_summary <- df %>%
select(all_of(key_vars)) %>%
summarise(across(everything(), ~sum(is.na(.)))) %>%
pivot_longer(everything(), names_to = "Variable", values_to = "N_Missing")
cat("=== MISSING DATA CHECK ===\n\n")## === MISSING DATA CHECK ===print(as.data.frame(missing_summary))## Variable N_Missing
## 1 Sentence_1 0
## 2 punishmore_1 0
## 3 punishmore_2 0
## 4 vignette 0
## 5 hate_1 0
## 6 revenge_1 0
## 7 rwa_1 0
## 8 sdo_1 0cat("\n")PHASE 2: REVERSE CODING
# Create working dataframe
df_clean <- df2.1 Punitiveness Items (reverse code so higher = more punitive)
# punishmore_1: "Current level of punishment is too high" - REVERSE
# (Agreeing means LESS punitive, so reverse)
df_clean$punishmore_1_R <- 8 - df_clean$punishmore_1
# parsimony_1: "Criminal punishments should be as light as possible" - REVERSE
df_clean$parsimony_1_R <- 8 - df_clean$parsimony_1
# parsimony_2: "The state should not punish more than necessary" - REVERSE
df_clean$parsimony_2_R <- 8 - df_clean$parsimony_2
cat("Punitiveness items reverse coded: punishmore_1, parsimony_1, parsimony_2\n")## Punitiveness items reverse coded: punishmore_1, parsimony_1, parsimony_22.2 Emotion Items
# hate_2: "I am sympathetic towards criminals" - REVERSE
# (Sympathy = opposite of hatred)
df_clean$hate_2_R <- 8 - df_clean$hate_2
cat("Emotion items reverse coded: hate_2 (sympathy)\n")## Emotion items reverse coded: hate_2 (sympathy)2.3 Harshness/Hostile Aggression Items
# degrad_1: "Prisoners deserve dignity and respect" - REVERSE
df_clean$degrad_1_R <- 8 - df_clean$degrad_1
# prisonvi_1: "I feel bad when a prisoner is assaulted" - REVERSE
# (Feeling bad = NOT tolerating violence)
df_clean$prisonvi_1_R <- 8 - df_clean$prisonvi_1
# prisonvi_2: "Sexual assault of prisoners is a serious problem" - REVERSE
# (Viewing it as serious = NOT tolerating violence)
df_clean$prisonvi_2_R <- 8 - df_clean$prisonvi_2
# harsh_2: "Prisoners deserve decent living conditions" - REVERSE
df_clean$harsh_2_R <- 8 - df_clean$harsh_2
cat("Harshness items reverse coded: degrad_1, prisonvi_1, prisonvi_2, harsh_2\n")## Harshness items reverse coded: degrad_1, prisonvi_1, prisonvi_2, harsh_22.4 Crime Concerns Items
# fear_1: "How safe do you feel walking at night" - REVERSE
# (Feeling safe = LOW fear)
df_clean$fear_1_R <- 8 - df_clean$fear_1
cat("Crime concerns items reverse coded: fear_1 (safety)\n")## Crime concerns items reverse coded: fear_1 (safety)2.5 Personality Items - RWA
# rwa_3: "Nothing wrong with premarital sex" - REVERSE
df_clean$rwa_3_R <- 8 - df_clean$rwa_3
cat("RWA items reverse coded: rwa_3\n")## RWA items reverse coded: rwa_32.6 Personality Items - SDO
# sdo_3: "No one group should dominate" - REVERSE (egalitarian)
df_clean$sdo_3_R <- 8 - df_clean$sdo_3
# sdo_4: "Bottom groups are just as deserving" - REVERSE (egalitarian)
df_clean$sdo_4_R <- 8 - df_clean$sdo_4
# sdo_7: "Should equalize conditions for groups" - REVERSE (egalitarian)
df_clean$sdo_7_R <- 8 - df_clean$sdo_7
# sdo_8: "Should give all groups equal chance" - REVERSE (egalitarian)
df_clean$sdo_8_R <- 8 - df_clean$sdo_8
cat("SDO items reverse coded: sdo_3, sdo_4, sdo_7, sdo_8\n")## SDO items reverse coded: sdo_3, sdo_4, sdo_7, sdo_82.7 Personality Items - Vengefulness
# venge_1: "Not worth my time to pay back someone" - REVERSE
df_clean$venge_1_R <- 8 - df_clean$venge_1
# venge_5: "I am not a vengeful person" - REVERSE
df_clean$venge_5_R <- 8 - df_clean$venge_5
cat("Vengefulness items reverse coded: venge_1, venge_5\n")## Vengefulness items reverse coded: venge_1, venge_52.8 Personality Items - Racial Resentment
# Raceresent_1: "Slavery/discrimination created difficult conditions" - REVERSE
# (Acknowledging systemic issues = LOW racial resentment)
df_clean$Raceresent_1_R <- 8 - df_clean$Raceresent_1
# Raceresent_4: "Blacks have gotten less than they deserve" - REVERSE
df_clean$Raceresent_4_R <- 8 - df_clean$Raceresent_4
cat("Racial Resentment items reverse coded: Raceresent_1, Raceresent_4\n")## Racial Resentment items reverse coded: Raceresent_1, Raceresent_42.9 Due Process Items (Exploratory - collected but cluster dropped)
# dueprocess_2: "Every defendant deserves due process" - REVERSE
# (Supporting due process = NOT tolerating violations)
df_clean$dueprocess_2_R <- 8 - df_clean$dueprocess_2
cat("Due Process items reverse coded: dueprocess_2\n")## Due Process items reverse coded: dueprocess_2PHASE 3: Z-SCORE SENTENCES WITHIN VIGNETTE
# Check vignette distribution
cat("=== VIGNETTE DISTRIBUTION ===\n\n")## === VIGNETTE DISTRIBUTION ===vignette_table <- table(df_clean$vignette, useNA = "ifany")
print(vignette_table)##
## 1 2 3
## 168 176 152cat("\nVignette Labels:\n")##
## Vignette Labels:cat(" 1 = Stranger Felony-Murder (purse snatching)\n")## 1 = Stranger Felony-Murder (purse snatching)cat(" 2 = Domestic Violence Murder (restraining order violation)\n")## 2 = Domestic Violence Murder (restraining order violation)cat(" 3 = Organized Crime Murder (car theft operation)\n\n")## 3 = Organized Crime Murder (car theft operation)# Descriptive stats by vignette BEFORE z-scoring
cat("=== SENTENCE DESCRIPTIVES BY VIGNETTE (Raw) ===\n\n")## === SENTENCE DESCRIPTIVES BY VIGNETTE (Raw) ===sentence_by_vignette <- df_clean %>%
group_by(vignette) %>%
summarise(
n = n(),
mean = mean(Sentence_1, na.rm = TRUE),
sd = sd(Sentence_1, na.rm = TRUE),
median = median(Sentence_1, na.rm = TRUE),
min = min(Sentence_1, na.rm = TRUE),
max = max(Sentence_1, na.rm = TRUE)
)
print(as.data.frame(sentence_by_vignette))## vignette n mean sd median min max
## 1 1 168 34.41 12.59 30 3 50
## 2 2 176 37.46 12.60 40 2 50
## 3 3 152 33.41 11.99 30 0 50cat("\n")# Z-score sentences within each vignette condition
# This removes vignette-specific variance while preserving individual differences
df_clean <- df_clean %>%
group_by(vignette) %>%
mutate(Sentence_z = scale(Sentence_1)[,1]) %>%
ungroup()
cat("Sentence_z created: z-scored within vignette condition\n")## Sentence_z created: z-scored within vignette condition# Verify z-scoring worked (means should be ~0, SDs should be ~1)
cat("\n=== SENTENCE DESCRIPTIVES BY VIGNETTE (Z-scored) ===\n\n")##
## === SENTENCE DESCRIPTIVES BY VIGNETTE (Z-scored) ===sentence_z_check <- df_clean %>%
group_by(vignette) %>%
summarise(
n = n(),
mean = mean(Sentence_z, na.rm = TRUE),
sd = sd(Sentence_z, na.rm = TRUE)
)
print(as.data.frame(sentence_z_check))## vignette n mean sd
## 1 1 168 0.000000000000000047581 1
## 2 2 176 -0.000000000000000008831 1
## 3 3 152 0.000000000000000150464 1cat("\n")PHASE 4: CREATE COMPOSITE SCORES WITH RELIABILITY
# Function to compute composite and report alpha
compute_composite <- function(data, items, composite_name, print_output = TRUE) {
# Select items
item_data <- data %>% select(all_of(items))
# Compute alpha
alpha_result <- psych::alpha(item_data, check.keys = FALSE)
alpha_value <- alpha_result$total$raw_alpha
# Compute composite (mean of items)
composite <- rowMeans(item_data, na.rm = TRUE)
if(print_output) {
cat(sprintf("%-30s: α = %.3f (k = %d)\n",
composite_name, alpha_value, length(items)))
}
return(list(composite = composite, alpha = alpha_value, n_items = length(items)))
}
# Store all alphas for later reporting
alpha_table <- data.frame(
Construct = character(),
Alpha = numeric(),
N_Items = integer(),
stringsAsFactors = FALSE
)4.1 PUNITIVENESS CONSTRUCT (DV)
cat("=== PUNITIVENESS (DV) ===\n\n")## === PUNITIVENESS (DV) ===# Punish More (2 items)
punishmore_items <- c("punishmore_1_R", "punishmore_2")
result <- compute_composite(df_clean, punishmore_items, "Punish More")## Punish More : α = 0.683 (k = 2)df_clean$punishmore_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Punish More",
Alpha = result$alpha, N_Items = result$n_items))
# Parsimony (2 items) - already reversed
parsimony_items <- c("parsimony_1_R", "parsimony_2_R")
result <- compute_composite(df_clean, parsimony_items, "Parsimony (R)")## Parsimony (R) : α = 0.477 (k = 2)df_clean$parsimony_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Parsimony",
Alpha = result$alpha, N_Items = result$n_items))
# Three Strikes (2 items)
threestrikes_items <- c("threestrikes_1", "threestrikes_2")
result <- compute_composite(df_clean, threestrikes_items, "Three Strikes")## Three Strikes : α = 0.751 (k = 2)df_clean$threestrikes_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Three Strikes",
Alpha = result$alpha, N_Items = result$n_items))
# Single items (report as single-item measures)
cat(sprintf("%-30s: Single item\n", "LWOP"))## LWOP : Single itemcat(sprintf("%-30s: Single item\n", "Death Penalty"))## Death Penalty : Single itemcat(sprintf("%-30s: Z-scored within vignette\n", "Sentence"))## Sentence : Z-scored within vignette# Punitiveness Aggregate (all items)
# Include: punishmore, parsimony, threestrikes, LWOP, deathpenalty, Sentence_z
punitiveness_items <- c("punishmore_1_R", "punishmore_2",
"parsimony_1_R", "parsimony_2_R",
"threestrikes_1", "threestrikes_2",
"LWOP", "deathpenalty")
# For aggregate, also include z-scored sentence
result <- compute_composite(df_clean, punitiveness_items, "Punitiveness (8 items)")## Punitiveness (8 items) : α = 0.842 (k = 8)alpha_table <- rbind(alpha_table, data.frame(Construct = "Punitiveness (8 items)",
Alpha = result$alpha, N_Items = result$n_items))
# Create punitiveness aggregate WITH sentence_z
# First standardize the 8-item composite, then average with Sentence_z
df_clean$punitiveness_8item <- result$composite
df_clean$punitiveness_8item_z <- scale(df_clean$punitiveness_8item)[,1]
df_clean$punitiveness_agg <- rowMeans(cbind(df_clean$punitiveness_8item_z,
df_clean$Sentence_z), na.rm = TRUE)
cat(sprintf("%-30s: Standardized composite + Sentence_z\n\n", "Punitiveness Aggregate"))## Punitiveness Aggregate : Standardized composite + Sentence_z4.2 CRIME CONCERNS CLUSTER
cat("=== CRIME CONCERNS CLUSTER ===\n\n")## === CRIME CONCERNS CLUSTER ===# Perceived Crime Rates (2 items)
rates_items <- c("rates_1", "rates_2")
result <- compute_composite(df_clean, rates_items, "Perceived Crime Rates")## Perceived Crime Rates : α = 0.845 (k = 2)df_clean$crime_rates_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Perceived Crime Rates",
Alpha = result$alpha, N_Items = result$n_items))
# Fear of Crime (3 items)
fear_items <- c("fear_1_R", "fear_2", "fear_3")
result <- compute_composite(df_clean, fear_items, "Fear of Crime")## Fear of Crime : α = 0.780 (k = 3)df_clean$fear_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Fear of Crime",
Alpha = result$alpha, N_Items = result$n_items))
# Crime Concerns Aggregate
crime_concerns_items <- c("rates_1", "rates_2", "fear_1_R", "fear_2", "fear_3")
result <- compute_composite(df_clean, crime_concerns_items, "Crime Concerns Cluster")## Crime Concerns Cluster : α = 0.782 (k = 5)df_clean$crime_concerns_agg <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Crime Concerns Cluster",
Alpha = result$alpha, N_Items = result$n_items))
cat("\n")4.3 EMOTIONS CLUSTER
cat("=== EMOTIONS CLUSTER ===\n\n")## === EMOTIONS CLUSTER ===# Hatred (3 items)
hatred_items <- c("hate_1", "hate_2_R", "hate_3")
result <- compute_composite(df_clean, hatred_items, "Hatred")## Hatred : α = 0.778 (k = 3)df_clean$hatred_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Hatred",
Alpha = result$alpha, N_Items = result$n_items))
# Anger (2 items)
anger_items <- c("anger_1", "anger_2")
result <- compute_composite(df_clean, anger_items, "Anger")## Anger : α = 0.828 (k = 2)df_clean$anger_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Anger",
Alpha = result$alpha, N_Items = result$n_items))
# Emotions Aggregate
emotions_items <- c("hate_1", "hate_2_R", "hate_3", "anger_1", "anger_2")
result <- compute_composite(df_clean, emotions_items, "Emotions Cluster")## Emotions Cluster : α = 0.865 (k = 5)df_clean$emotions_agg <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Emotions Cluster",
Alpha = result$alpha, N_Items = result$n_items))
cat("\n")4.4 HARSHNESS / HOSTILE AGGRESSION CLUSTER
cat("=== HARSHNESS / HOSTILE AGGRESSION CLUSTER ===\n\n")## === HARSHNESS / HOSTILE AGGRESSION CLUSTER ===# Social Exclusion (3 items)
exclusion_items <- c("exclusion_1", "exclusion_2", "exclusion_3")
result <- compute_composite(df_clean, exclusion_items, "Social Exclusion")## Social Exclusion : α = 0.761 (k = 3)df_clean$exclusion_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Social Exclusion",
Alpha = result$alpha, N_Items = result$n_items))
# Degradation (3 items)
degradation_items <- c("degrad_1_R", "degrad_2", "degrad_3")
result <- compute_composite(df_clean, degradation_items, "Degradation")## Degradation : α = 0.670 (k = 3)df_clean$degradation_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Degradation",
Alpha = result$alpha, N_Items = result$n_items))
# Infliction of Suffering (2 items)
suffering_items <- c("suffer_1", "suffer_2")
result <- compute_composite(df_clean, suffering_items, "Infliction of Suffering")## Infliction of Suffering : α = 0.794 (k = 2)df_clean$suffering_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Infliction of Suffering",
Alpha = result$alpha, N_Items = result$n_items))
# Prison Violence Tolerance (2 items) - reversed items
prisonvi_items <- c("prisonvi_1_R", "prisonvi_2_R")
result <- compute_composite(df_clean, prisonvi_items, "Prison Violence Tolerance")## Prison Violence Tolerance : α = 0.559 (k = 2)df_clean$prisonvi_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Prison Violence Tolerance",
Alpha = result$alpha, N_Items = result$n_items))
# Harsh Prison Conditions (3 items)
harsh_items <- c("harsh_1", "harsh_2_R", "harsh_3")
result <- compute_composite(df_clean, harsh_items, "Harsh Prison Conditions")## Harsh Prison Conditions : α = 0.818 (k = 3)df_clean$harsh_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Harsh Prison Conditions",
Alpha = result$alpha, N_Items = result$n_items))
# Revenge (3 items)
revenge_items <- c("revenge_1", "revenge_2", "revenge_3")
result <- compute_composite(df_clean, revenge_items, "Revenge")## Revenge : α = 0.665 (k = 3)df_clean$revenge_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Revenge",
Alpha = result$alpha, N_Items = result$n_items))
# Hostile Aggression Aggregate
hostile_agg_items <- c("exclusion_1", "exclusion_2", "exclusion_3",
"degrad_1_R", "degrad_2", "degrad_3",
"suffer_1", "suffer_2",
"prisonvi_1_R", "prisonvi_2_R",
"harsh_1", "harsh_2_R", "harsh_3",
"revenge_1", "revenge_2", "revenge_3")
result <- compute_composite(df_clean, hostile_agg_items, "Harshness Cluster")## Harshness Cluster : α = 0.921 (k = 16)df_clean$hostile_agg <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Harshness Cluster",
Alpha = result$alpha, N_Items = result$n_items))
cat("\n")4.5 PERSONALITY CLUSTER
cat("=== PERSONALITY CLUSTER ===\n\n")## === PERSONALITY CLUSTER ===# Right-Wing Authoritarianism (5 items)
rwa_items <- c("rwa_1", "rwa_2", "rwa_3_R", "rwa_4", "rwa_5")
result <- compute_composite(df_clean, rwa_items, "RWA")## RWA : α = 0.869 (k = 5)df_clean$rwa_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "RWA",
Alpha = result$alpha, N_Items = result$n_items))
# Social Dominance Orientation (8 items)
sdo_items <- c("sdo_1", "sdo_2", "sdo_3_R", "sdo_4_R",
"sdo_5", "sdo_6", "sdo_7_R", "sdo_8_R")
result <- compute_composite(df_clean, sdo_items, "SDO")## SDO : α = 0.923 (k = 8)df_clean$sdo_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "SDO",
Alpha = result$alpha, N_Items = result$n_items))
# Vengefulness (5 items)
venge_items <- c("venge_1_R", "venge_2", "venge_3", "venge_4", "venge_5_R")
result <- compute_composite(df_clean, venge_items, "Vengefulness")## Vengefulness : α = 0.879 (k = 5)df_clean$venge_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Vengefulness",
Alpha = result$alpha, N_Items = result$n_items))
# Violence Proneness (4 items)
vprone_items <- c("vprone_1", "vprone_2", "vprone_3", "vprone_4")
result <- compute_composite(df_clean, vprone_items, "Violence Proneness")## Violence Proneness : α = 0.729 (k = 4)df_clean$vprone_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Violence Proneness",
Alpha = result$alpha, N_Items = result$n_items))
# Racial Resentment (4 items)
raceresent_items <- c("Raceresent_1_R", "Raceresent_2", "Raceresent_3", "Raceresent_4_R")
result <- compute_composite(df_clean, raceresent_items, "Racial Resentment")## Racial Resentment : α = 0.904 (k = 4)df_clean$raceresent_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Racial Resentment",
Alpha = result$alpha, N_Items = result$n_items))
# Blood Sports (4 items: Boxing, UFC, Hunting, Wrestling)
# TV_4 = Boxing, TV_6 = UFC, TV_8 = Hunting, TV_9 = Wrestling
bloodsports_items <- c("TV_4", "TV_6", "TV_8", "TV_9")
result <- compute_composite(df_clean, bloodsports_items, "Blood Sports")## Blood Sports : α = 0.837 (k = 4)df_clean$bloodsports_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Blood Sports",
Alpha = result$alpha, N_Items = result$n_items))
# Personality Aggregate (all personality measures)
personality_items <- c(rwa_items, sdo_items, venge_items, vprone_items,
raceresent_items, bloodsports_items)
result <- compute_composite(df_clean, personality_items, "Personality Cluster")## Personality Cluster : α = 0.924 (k = 30)df_clean$personality_agg <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Personality Cluster",
Alpha = result$alpha, N_Items = result$n_items))
cat("\n")4.6 PROCESS VIOLATIONS (EXPLORATORY ONLY)
cat("=== PROCESS VIOLATIONS (Exploratory Only) ===\n\n")## === PROCESS VIOLATIONS (Exploratory Only) ===# Due Process (3 items)
dueprocess_items <- c("dueprocess_1", "dueprocess_2_R", "dueprocess_3")
result <- compute_composite(df_clean, dueprocess_items, "Due Process Violations")## Due Process Violations : α = 0.598 (k = 3)df_clean$dueprocess_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Due Process Violations*",
Alpha = result$alpha, N_Items = result$n_items))
# Uncertain Evidence (4 items)
uncertain_items <- c("uncertain_1", "uncertain_2", "uncertain_3", "uncertain_4")
result <- compute_composite(df_clean, uncertain_items, "Uncertain Evidence")## Uncertain Evidence : α = 0.667 (k = 4)df_clean$uncertain_comp <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Uncertain Evidence*",
Alpha = result$alpha, N_Items = result$n_items))
# Process Violations Aggregate
process_items <- c(dueprocess_items, uncertain_items)
result <- compute_composite(df_clean, process_items, "Process Violations Cluster")## Process Violations Cluster : α = 0.758 (k = 7)df_clean$process_agg <- result$composite
alpha_table <- rbind(alpha_table, data.frame(Construct = "Process Violations Cluster*",
Alpha = result$alpha, N_Items = result$n_items))
cat("\n* = Exploratory only (cluster dropped from confirmatory analyses per pre-registration)\n\n")##
## * = Exploratory only (cluster dropped from confirmatory analyses per pre-registration)4.7 Print Complete Alpha Table
cat("\n=== COMPLETE RELIABILITY TABLE ===\n\n")##
## === COMPLETE RELIABILITY TABLE ===print(as.data.frame(alpha_table))## Construct Alpha N_Items
## 1 Punish More 0.6832 2
## 2 Parsimony 0.4773 2
## 3 Three Strikes 0.7510 2
## 4 Punitiveness (8 items) 0.8421 8
## 5 Perceived Crime Rates 0.8453 2
## 6 Fear of Crime 0.7802 3
## 7 Crime Concerns Cluster 0.7824 5
## 8 Hatred 0.7784 3
## 9 Anger 0.8280 2
## 10 Emotions Cluster 0.8649 5
## 11 Social Exclusion 0.7615 3
## 12 Degradation 0.6703 3
## 13 Infliction of Suffering 0.7943 2
## 14 Prison Violence Tolerance 0.5594 2
## 15 Harsh Prison Conditions 0.8182 3
## 16 Revenge 0.6653 3
## 17 Harshness Cluster 0.9205 16
## 18 RWA 0.8693 5
## 19 SDO 0.9225 8
## 20 Vengefulness 0.8793 5
## 21 Violence Proneness 0.7291 4
## 22 Racial Resentment 0.9041 4
## 23 Blood Sports 0.8365 4
## 24 Personality Cluster 0.9237 30
## 25 Due Process Violations* 0.5976 3
## 26 Uncertain Evidence* 0.6669 4
## 27 Process Violations Cluster* 0.7581 7# Flag low reliability composites
low_alpha <- alpha_table %>% filter(Alpha < 0.60)
if(nrow(low_alpha) > 0) {
cat("\n WARNING: Composites with α < .60 (interpret with caution): \n")
print(as.data.frame(low_alpha))
}##
## WARNING: Composites with α < .60 (interpret with caution):
## Construct Alpha N_Items
## 1 Parsimony 0.4773 2
## 2 Prison Violence Tolerance 0.5594 2
## 3 Due Process Violations* 0.5976 3cat("\n")PHASE 5: DESCRIPTIVE STATISTICS
5.1 Sample Demographics
cat("=== SAMPLE DEMOGRAPHICS ===\n\n")## === SAMPLE DEMOGRAPHICS ===# Age
cat("AGE:\n")## AGE:cat(sprintf(" Mean: %.1f years\n", mean(df_clean$age_1, na.rm = TRUE)))## Mean: 46.2 yearscat(sprintf(" SD: %.1f\n", sd(df_clean$age_1, na.rm = TRUE)))## SD: 16.2cat(sprintf(" Range: %d - %d\n\n",
min(df_clean$age_1, na.rm = TRUE),
max(df_clean$age_1, na.rm = TRUE)))## Range: 19 - 85# Gender (1=Male, 2=Female, 3=Other)
cat("GENDER:\n")## GENDER:gender_table <- df_clean %>%
count(gender) %>%
mutate(
Gender_Label = case_when(
gender == 1 ~ "Male",
gender == 2 ~ "Female",
gender == 3 ~ "Other",
TRUE ~ "Missing"
),
Percent = sprintf("%.1f%%", n/sum(n)*100)
) %>%
select(Gender_Label, n, Percent)
print(as.data.frame(gender_table))## Gender_Label n Percent
## 1 Male 244 49.2%
## 2 Female 249 50.2%
## 3 Other 3 0.6%cat("\n")# Race/Ethnicity
cat("RACE/ETHNICITY:\n")## RACE/ETHNICITY:race_table <- df_clean %>%
count(race) %>%
mutate(
Race_Label = case_when(
race == 1 ~ "White, non-Hispanic",
race == 2 ~ "Black/African-American",
race == 3 ~ "Asian",
race == 4 ~ "Hispanic/Latino",
race == 5 ~ "Middle Eastern",
race == 6 ~ "Native American",
race == 7 ~ "Pacific Islander",
TRUE ~ "Missing/Other"
),
Percent = sprintf("%.1f%%", n/sum(n)*100)
) %>%
select(Race_Label, n, Percent) %>%
arrange(desc(n))
print(as.data.frame(race_table))## Race_Label n Percent
## 1 White, non-Hispanic 324 65.3%
## 2 Black/African-American 66 13.3%
## 3 Hispanic/Latino 53 10.7%
## 4 Asian 39 7.9%
## 5 Native American 8 1.6%
## 6 Middle Eastern 4 0.8%
## 7 Pacific Islander 2 0.4%cat("\n")# Education
cat("EDUCATION:\n")## EDUCATION:education_table <- df_clean %>%
count(education) %>%
mutate(
Education_Label = case_when(
education == 1 ~ "Did not graduate high school",
education == 2 ~ "High school only",
education == 3 ~ "1-2 years post HS",
education == 4 ~ "3-4 years post HS",
education == 5 ~ "5-6 years post HS",
education == 6 ~ "More than 6 years post HS",
TRUE ~ "Missing"
),
Percent = sprintf("%.1f%%", n/sum(n)*100)
) %>%
select(Education_Label, n, Percent)
print(as.data.frame(education_table))## Education_Label n Percent
## 1 Did not graduate high school 4 0.8%
## 2 High school only 74 14.9%
## 3 1-2 years post HS 74 14.9%
## 4 3-4 years post HS 164 33.1%
## 5 5-6 years post HS 58 11.7%
## 6 More than 6 years post HS 122 24.6%cat("\n")# Political Orientation (1=Strongly liberal to 7=Strongly conservative)
cat("POLITICAL ORIENTATION:\n")## POLITICAL ORIENTATION:cat(sprintf(" Mean: %.2f (1=Strongly liberal, 7=Strongly conservative)\n",
mean(df_clean$politid, na.rm = TRUE)))## Mean: 4.00 (1=Strongly liberal, 7=Strongly conservative)cat(sprintf(" SD: %.2f\n", sd(df_clean$politid, na.rm = TRUE)))## SD: 1.88cat(sprintf(" Median: %.0f\n\n", median(df_clean$politid, na.rm = TRUE)))## Median: 45.2 Construct Descriptive Statistics
cat("=== CONSTRUCT DESCRIPTIVE STATISTICS ===\n\n")## === CONSTRUCT DESCRIPTIVE STATISTICS ===# List of all composites for descriptives
construct_vars_desc <- c(
# Punitiveness
"punishmore_comp", "parsimony_comp", "threestrikes_comp",
"LWOP", "deathpenalty", "Sentence_1", "punitiveness_agg",
# Crime Concerns
"crime_rates_comp", "fear_comp", "crime_concerns_agg",
# Emotions
"hatred_comp", "anger_comp", "emotions_agg",
# Harshness
"exclusion_comp", "degradation_comp", "suffering_comp",
"prisonvi_comp", "harsh_comp", "revenge_comp", "hostile_agg",
# Personality
"rwa_comp", "sdo_comp", "venge_comp", "vprone_comp",
"raceresent_comp", "bloodsports_comp", "personality_agg",
# Process (Exploratory)
"dueprocess_comp", "uncertain_comp", "process_agg"
)
# Function to compute midpoint percentages
compute_midpoint_pcts <- function(x, midpoint) {
n_valid <- sum(!is.na(x))
c(
Pct_Below = round(sum(x < midpoint, na.rm = TRUE) / n_valid, 2),
Pct_At = round(sum(x == midpoint, na.rm = TRUE) / n_valid, 2),
Pct_Above = round(sum(x > midpoint, na.rm = TRUE) / n_valid, 2)
)
}
# Compute descriptives with midpoint analysis
descriptives_list <- lapply(construct_vars_desc, function(var) {
x <- df_clean[[var]]
# Determine midpoint based on variable
if (var == "Sentence_1") {
midpoint <- 25 # 0-50 scale
} else if (var == "punitiveness_agg") {
midpoint <- 0 # z-scored, midpoint is 0
} else if (var %in% c("crime_concerns_agg", "emotions_agg", "hostile_agg",
"personality_agg", "process_agg")) {
midpoint <- 4 # aggregates of 1-7 scales
} else {
midpoint <- 4 # standard 1-7 scale
}
# Basic stats
stats <- data.frame(
Variable = var,
Mean = round(mean(x, na.rm = TRUE), 2),
SD = round(sd(x, na.rm = TRUE), 2),
Median = round(median(x, na.rm = TRUE), 2),
Min = round(min(x, na.rm = TRUE), 2),
Max = round(max(x, na.rm = TRUE), 2),
N = sum(!is.na(x))
)
# Midpoint percentages
pcts <- compute_midpoint_pcts(x, midpoint)
stats$Pct_Below <- pcts["Pct_Below"]
stats$Pct_At <- pcts["Pct_At"]
stats$Pct_Above <- pcts["Pct_Above"]
return(stats)
})
descriptives_table <- do.call(rbind, descriptives_list)
rownames(descriptives_table) <- NULL
cat("Descriptive Statistics with Midpoint Analysis:\n")## Descriptive Statistics with Midpoint Analysis:cat("(Midpoint = 4 for 1-7 scales; 25 for Sentence; 0 for z-scored punitiveness_agg)\n")## (Midpoint = 4 for 1-7 scales; 25 for Sentence; 0 for z-scored punitiveness_agg)cat("(Pct columns show proportions below/at/above midpoint)\n\n")## (Pct columns show proportions below/at/above midpoint)print(as.data.frame(descriptives_table), row.names = FALSE)## Variable Mean SD Median Min Max N Pct_Below Pct_At
## punishmore_comp 4.12 1.63 4.00 1.00 7.00 496 0.40 0.13
## parsimony_comp 4.21 1.31 4.00 1.00 7.00 496 0.33 0.21
## threestrikes_comp 3.20 1.69 3.00 1.00 7.00 496 0.62 0.12
## LWOP 4.99 1.78 5.00 1.00 7.00 496 0.20 0.16
## deathpenalty 4.31 2.07 5.00 1.00 7.00 496 0.33 0.15
## Sentence_1 35.19 12.51 35.00 0.00 50.00 496 0.15 0.15
## punitiveness_agg 0.00 0.82 -0.01 -2.18 1.74 496 0.50 0.00
## crime_rates_comp 5.00 1.42 5.00 1.00 7.00 496 0.19 0.09
## fear_comp 3.02 1.34 2.67 1.00 7.00 496 0.70 0.09
## crime_concerns_agg 3.81 1.15 3.80 1.00 7.00 496 0.54 0.08
## hatred_comp 4.34 1.44 4.33 1.00 7.00 496 0.38 0.08
## anger_comp 4.84 1.51 5.00 1.00 7.00 496 0.22 0.12
## emotions_agg 4.54 1.36 4.60 1.00 7.00 496 0.32 0.06
## exclusion_comp 4.46 1.45 4.67 1.00 7.00 496 0.32 0.09
## degradation_comp 3.66 1.37 3.67 1.00 7.00 496 0.58 0.09
## suffering_comp 3.82 1.68 4.00 1.00 7.00 496 0.45 0.12
## prisonvi_comp 3.04 1.38 3.00 1.00 7.00 496 0.69 0.14
## harsh_comp 3.30 1.55 3.33 1.00 7.00 496 0.64 0.08
## revenge_comp 4.11 1.43 4.00 1.00 7.00 496 0.40 0.11
## hostile_agg 3.77 1.21 3.84 1.00 7.00 496 0.56 0.03
## rwa_comp 3.72 1.65 3.80 1.00 7.00 496 0.51 0.05
## sdo_comp 2.58 1.44 2.38 1.00 7.00 496 0.81 0.03
## venge_comp 2.97 1.41 2.80 1.00 7.00 496 0.74 0.04
## vprone_comp 3.29 1.41 3.25 1.00 7.00 496 0.66 0.05
## raceresent_comp 3.52 1.76 3.50 1.00 7.00 496 0.57 0.10
## bloodsports_comp 3.17 1.60 3.25 1.00 7.00 496 0.64 0.05
## personality_agg 3.13 1.05 3.17 1.00 6.33 496 0.77 0.01
## dueprocess_comp 3.02 1.18 3.00 1.00 7.00 496 0.72 0.10
## uncertain_comp 3.17 1.12 3.25 1.00 6.50 496 0.73 0.07
## process_agg 3.11 1.01 3.14 1.00 6.29 496 0.77 0.04
## Pct_Above
## 0.47
## 0.46
## 0.26
## 0.64
## 0.52
## 0.70
## 0.50
## 0.72
## 0.21
## 0.38
## 0.54
## 0.66
## 0.62
## 0.59
## 0.33
## 0.43
## 0.18
## 0.28
## 0.49
## 0.42
## 0.43
## 0.15
## 0.22
## 0.29
## 0.33
## 0.30
## 0.22
## 0.18
## 0.19
## 0.18cat("\n")5.3 Sentence Distribution
cat("=== SENTENCE DISTRIBUTION ===\n\n")## === SENTENCE DISTRIBUTION ===cat("Overall:\n")## Overall:cat(sprintf(" Mean: %.1f years\n", mean(df_clean$Sentence_1, na.rm = TRUE)))## Mean: 35.2 yearscat(sprintf(" SD: %.1f\n", sd(df_clean$Sentence_1, na.rm = TRUE)))## SD: 12.5cat(sprintf(" Median: %.0f\n", median(df_clean$Sentence_1, na.rm = TRUE)))## Median: 35# Percentage at key values
cat(sprintf(" %% at maximum (50 years): %.1f%%\n",
mean(df_clean$Sentence_1 == 50, na.rm = TRUE) * 100))## % at maximum (50 years): 32.3%cat(sprintf(" %% at 0 years: %.1f%%\n",
mean(df_clean$Sentence_1 == 0, na.rm = TRUE) * 100))## % at 0 years: 0.2%cat("\nBy Vignette:\n")##
## By Vignette:print(as.data.frame(sentence_by_vignette))## vignette n mean sd median min max
## 1 1 168 34.41 12.59 30 3 50
## 2 2 176 37.46 12.60 40 2 50
## 3 3 152 33.41 11.99 30 0 50cat("\n")5.4 Sentence Distribution Visualizations
cat("=== SENTENCE DISTRIBUTION HISTOGRAMS ===\n\n")## === SENTENCE DISTRIBUTION HISTOGRAMS ===# Overall histogram
png("histogram_sentence.png", width = 800, height = 600, res = 150)
ggplot(df_clean, aes(x = Sentence_1)) +
geom_histogram(binwidth = 5, fill = "steelblue", color = "white") +
geom_vline(xintercept = mean(df_clean$Sentence_1, na.rm = TRUE),
color = "red", linetype = "dashed", size = 1) +
labs(title = "Distribution of Sentence Recommendations",
subtitle = sprintf("M = %.1f, SD = %.1f, Median = %.0f",
mean(df_clean$Sentence_1, na.rm = TRUE),
sd(df_clean$Sentence_1, na.rm = TRUE),
median(df_clean$Sentence_1, na.rm = TRUE)),
x = "Sentence (Years)", y = "Count") +
theme_minimal() +
scale_x_continuous(breaks = seq(0, 50, 10))
dev.off()## quartz_off_screen
## 2cat("Saved: histogram_sentence.png\n")## Saved: histogram_sentence.png# By vignette
png("histogram_sentence_by_vignette.png", width = 1000, height = 400, res = 150)
ggplot(df_clean, aes(x = Sentence_1, fill = factor(vignette))) +
geom_histogram(binwidth = 5, color = "white") +
facet_wrap(~vignette, labeller = labeller(vignette = c(
"1" = "Stranger Felony-Murder",
"2" = "Domestic Violence",
"3" = "Organized Crime"
))) +
labs(title = "Sentence Distribution by Vignette",
x = "Sentence (Years)", y = "Count") +
theme_minimal() +
theme(legend.position = "none") +
scale_fill_brewer(palette = "Set2")
dev.off()## quartz_off_screen
## 2cat("Saved: histogram_sentence_by_vignette.png\n\n")## Saved: histogram_sentence_by_vignette.pngPHASE 6: CONFIRMATORY HYPOTHESIS TESTS
6.2 H2: Harshness > Crime Concerns
cat("=== H2: Harshness > Crime Concerns in predicting punitiveness ===\n\n")## === H2: Harshness > Crime Concerns in predicting punitiveness ===# Correlations
r_hostile <- cor(df_clean$punitiveness_agg, df_clean$hostile_agg,
use = "pairwise.complete.obs")
r_crime <- cor(df_clean$punitiveness_agg, df_clean$crime_concerns_agg,
use = "pairwise.complete.obs")
r_hostile_crime <- cor(df_clean$hostile_agg, df_clean$crime_concerns_agg,
use = "pairwise.complete.obs")
cat(sprintf("r(Punitiveness, Harshness) = %.3f\n", r_hostile))## r(Punitiveness, Harshness) = 0.586cat(sprintf("r(Punitiveness, Crime Concerns) = %.3f\n", r_crime))## r(Punitiveness, Crime Concerns) = 0.322cat(sprintf("r(Harshness, Crime Concerns) = %.3f\n\n", r_hostile_crime))## r(Harshness, Crime Concerns) = 0.380# Steiger's Z-test for dependent correlations
# Using cocor package
n_steiger <- sum(complete.cases(df_clean[, c("punitiveness_agg", "hostile_agg", "crime_concerns_agg")]))
steiger_result <- cocor.dep.groups.overlap(
r.jk = r_hostile, # r(punitiveness, hostile)
r.jh = r_crime, # r(punitiveness, crime_concerns)
r.kh = r_hostile_crime, # r(hostile, crime_concerns)
n = n_steiger,
test = "steiger1980"
)
cat("Steiger's Z-test for Dependent Correlations:\n")## Steiger's Z-test for Dependent Correlations:cat(sprintf(" N = %d\n", n_steiger))## N = 496cat(sprintf(" Difference: %.3f - %.3f = %.3f\n", r_hostile, r_crime, r_hostile - r_crime))## Difference: 0.586 - 0.322 = 0.264print(steiger_result)##
## Results of a comparison of two overlapping correlations based on dependent groups
##
## Comparison between r.jk = 0.5858 and r.jh = 0.3216
## Difference: r.jk - r.jh = 0.2642
## Related correlation: r.kh = 0.3805
## Group size: n = 496
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
##
## steiger1980: Steiger's (1980) modification of Dunn and Clark's z (1969) using average correlations
## z = 6.2614, p-value = 0.0000
## Null hypothesis rejected# Extract p-value safely
h2_pvalue <- steiger_result@steiger1980$p.value
cat(sprintf("\n*** H2 SUMMARY ***\n"))##
## *** H2 SUMMARY ***cat(sprintf(" Harshness r = %.3f\n", r_hostile))## Harshness r = 0.586cat(sprintf(" Crime Concerns r = %.3f\n", r_crime))## Crime Concerns r = 0.322cat(sprintf(" Difference = %.3f\n", r_hostile - r_crime))## Difference = 0.264cat(sprintf(" Steiger's Z p-value = %.4f\n", h2_pvalue))## Steiger's Z p-value = 0.0000cat(sprintf(" H2 SUPPORTED: %s\n\n",
ifelse(r_hostile > r_crime & h2_pvalue < .05,
"YES (significantly greater)", "CHECK P-VALUE")))## H2 SUPPORTED: YES (significantly greater)6.2b: Extended H2 - All Clusters vs Crime Concerns
cat("=== EXTENDED H2: All Clusters vs Crime Concerns ===\n\n")## === EXTENDED H2: All Clusters vs Crime Concerns ===clusters_to_compare <- c("hostile_agg", "emotions_agg", "personality_agg", "process_agg")
cluster_labels_h2 <- c("Harshness", "Emotions", "Personality", "Process Violations")
for(i in 1:length(clusters_to_compare)) {
r_other <- cor(df_clean$punitiveness_agg, df_clean[[clusters_to_compare[i]]],
use = "pairwise.complete.obs")
r_between <- cor(df_clean[[clusters_to_compare[i]]], df_clean$crime_concerns_agg,
use = "pairwise.complete.obs")
steiger <- cocor.dep.groups.overlap(
r.jk = r_other,
r.jh = r_crime,
r.kh = r_between,
n = n_steiger,
test = "steiger1980"
)
cat(sprintf("%s vs Crime Concerns:\n", cluster_labels_h2[i]))
cat(sprintf(" r(%s) = %.3f vs r(Crime) = %.3f\n", cluster_labels_h2[i], r_other, r_crime))
cat(sprintf(" Difference = %.3f, p = %.4f\n", r_other - r_crime,
steiger@steiger1980$p.value))
cat(sprintf(" %s > Crime Concerns: %s\n\n", cluster_labels_h2[i],
ifelse(steiger@steiger1980$p.value < .05 & r_other > r_crime, "YES", "NO")))
}## Harshness vs Crime Concerns:
## r(Harshness) = 0.586 vs r(Crime) = 0.322
## Difference = 0.264, p = 0.0000
## Harshness > Crime Concerns: YES
##
## Emotions vs Crime Concerns:
## r(Emotions) = 0.535 vs r(Crime) = 0.322
## Difference = 0.213, p = 0.0000
## Emotions > Crime Concerns: YES
##
## Personality vs Crime Concerns:
## r(Personality) = 0.475 vs r(Crime) = 0.322
## Difference = 0.153, p = 0.0012
## Personality > Crime Concerns: YES
##
## Process Violations vs Crime Concerns:
## r(Process Violations) = 0.439 vs r(Crime) = 0.322
## Difference = 0.117, p = 0.0091
## Process Violations > Crime Concerns: YES6.3 H3: Intercorrelations Among Correlates
cat("=== H3: Correlates positively intercorrelated ===\n\n")## === H3: Correlates positively intercorrelated ===# Create matrix of all correlate constructs (excluding punitiveness)
correlate_vars <- c(
"crime_rates_comp", "fear_comp",
"dueprocess_comp", "uncertain_comp",
"hatred_comp", "anger_comp",
"exclusion_comp", "degradation_comp", "suffering_comp",
"prisonvi_comp", "harsh_comp", "revenge_comp",
"rwa_comp", "sdo_comp", "venge_comp", "vprone_comp",
"raceresent_comp", "bloodsports_comp"
)
correlate_data <- df_clean %>% select(all_of(correlate_vars)) %>% na.omit()
correlate_cor <- corr.test(correlate_data, use = "pairwise", method = "pearson")
# Apply clean labels
colnames(correlate_cor$r) <- rownames(correlate_cor$r) <- construct_labels[colnames(correlate_cor$r)]
colnames(correlate_cor$p) <- rownames(correlate_cor$p) <- construct_labels[colnames(correlate_cor$p)]
# Count off-diagonal correlations
n_vars <- length(correlate_vars)
h3_n_correlations <- (n_vars * (n_vars - 1)) / 2 # Off-diagonal
# Extract upper triangle
upper_r <- correlate_cor$r[upper.tri(correlate_cor$r)]
upper_p <- correlate_cor$p[upper.tri(correlate_cor$p)]
h3_n_positive <- sum(upper_r > 0)
h3_n_sig_positive <- sum(upper_r > 0 & upper_p < .05)
h3_n_sig <- sum(upper_p < .05)
cat(sprintf("Total off-diagonal correlations: %d\n", h3_n_correlations))## Total off-diagonal correlations: 153cat(sprintf("Positive correlations: %d (%.1f%%)\n",
h3_n_positive, h3_n_positive/h3_n_correlations*100))## Positive correlations: 153 (100.0%)cat(sprintf("Significant (p < .05): %d (%.1f%%)\n",
h3_n_sig, h3_n_sig/h3_n_correlations*100))## Significant (p < .05): 144 (94.1%)cat(sprintf("Significant AND positive: %d (%.1f%%)\n",
h3_n_sig_positive, h3_n_sig_positive/h3_n_correlations*100))## Significant AND positive: 144 (94.1%)cat(sprintf("\n*** H3 SUMMARY ***\n"))##
## *** H3 SUMMARY ***cat(sprintf(" H3 SUPPORTED: %s\n\n",
ifelse(h3_n_positive/h3_n_correlations > 0.5, "YES (majority positive)", "NO")))## H3 SUPPORTED: YES (majority positive)# Full intercorrelation heat map
png("heatmap_intercorrelations.png", width = 1600, height = 1600, res = 150)
corrplot(correlate_cor$r,
method = "color",
type = "upper",
addCoef.col = "black",
number.cex = 0.45,
tl.col = "black",
tl.srt = 45,
tl.cex = 0.55,
p.mat = correlate_cor$p,
sig.level = 0.05,
insig = "n",
title = "Intercorrelations Among Correlate Constructs",
mar = c(0,0,2,0))
dev.off()## quartz_off_screen
## 2cat("Saved: heatmap_intercorrelations.png\n\n")## Saved: heatmap_intercorrelations.png6.4 Apply FDR Correction
cat("=== FDR CORRECTION FOR MULTIPLE COMPARISONS ===\n\n")## === FDR CORRECTION FOR MULTIPLE COMPARISONS ===# Collect all p-values from H1 construct-level tests
all_pvalues <- punitiveness_cors$p
fdr_adjusted <- p.adjust(all_pvalues, method = "fdr")
fdr_table <- punitiveness_cors %>%
mutate(
p_fdr = round(fdr_adjusted, 4),
sig_fdr = case_when(
p_fdr < .001 ~ "***",
p_fdr < .01 ~ "**",
p_fdr < .05 ~ "*",
TRUE ~ ""
)
) %>%
select(Construct, r, p, sig, p_fdr, sig_fdr)
cat("H1 Results with FDR Correction:\n")## H1 Results with FDR Correction:print(as.data.frame(fdr_table), row.names = FALSE)## Construct r p sig p_fdr
## Harsh Conditions 0.56 0.00000000000000000000000000000000000000146 *** 0.0000
## Exclusion 0.54 0.00000000000000000000000000000000000088367 *** 0.0000
## Hatred 0.53 0.00000000000000000000000000000000003408773 *** 0.0000
## Suffering 0.52 0.00000000000000000000000000000000692216355 *** 0.0000
## RWA 0.50 0.00000000000000000000000000000255468690156 *** 0.0000
## Racial Resentment 0.46 0.00000000000000000000000007576270968897927 *** 0.0000
## Anger 0.45 0.00000000000000000000000774325910360865235 *** 0.0000
## Degradation 0.44 0.00000000000000000000002751006908269931110 *** 0.0000
## Due Process 0.43 0.00000000000000000000050490585795064224687 *** 0.0000
## Revenge 0.43 0.00000000000000000000089296391036317302573 *** 0.0000
## Violence Proneness 0.41 0.00000000000000000004492528069221440832656 *** 0.0000
## Crime Rates 0.40 0.00000000000000000258996473346146174903954 *** 0.0000
## Uncertain Evidence 0.35 0.00000000000004045873073833809418642079117 *** 0.0000
## Prison Violence 0.35 0.00000000000006968182401471611736698009688 *** 0.0000
## SDO 0.32 0.00000000001179328365267163918690994028362 *** 0.0000
## Fear 0.18 0.00132377313759396300976711735586377471918 ** 0.0015
## Vengefulness 0.12 0.07367435580997361288524416522704996168613 0.0780
## Blood Sports 0.11 0.08017318731704685397243537181566352955997 0.0802
## sig_fdr
## ***
## ***
## ***
## ***
## ***
## ***
## ***
## ***
## ***
## ***
## ***
## ***
## ***
## ***
## ***
## **
##
## h1_n_sig_fdr <- sum(fdr_adjusted < .05 & punitiveness_cors$r > 0)
cat(sprintf("\nSignificant after FDR correction: %d/%d\n\n", h1_n_sig_fdr, nrow(fdr_table)))##
## Significant after FDR correction: 16/18PHASE 7: CORRELATION TABLES (2 decimal places)
# Function to create formatted correlation table
format_cor_table <- function(cor_matrix, p_matrix, decimals = 2) {
# Round correlations
r_formatted <- round(cor_matrix, decimals)
# Add significance stars
sig_stars <- ifelse(p_matrix < .001, "***",
ifelse(p_matrix < .01, "**",
ifelse(p_matrix < .05, "*", "")))
# Combine
result <- matrix(paste0(format(r_formatted, nsmall = decimals), sig_stars),
nrow = nrow(cor_matrix),
dimnames = dimnames(cor_matrix))
# Set diagonal to "-"
diag(result) <- "-"
return(result)
}
# --- Cluster-Level Table ---
cat("=== CLUSTER-LEVEL CORRELATION TABLE ===\n\n")## === CLUSTER-LEVEL CORRELATION TABLE ===cluster_table <- format_cor_table(cluster_cor$r, cluster_cor$p)
print(noquote(cluster_table))## Punitiveness Crime Concerns Process Violations Emotions
## Punitiveness - 0.32*** 0.44*** 0.53***
## Crime Concerns 0.32*** - 0.38*** 0.43***
## Process Violations 0.44*** 0.38*** - 0.48***
## Emotions 0.53*** 0.43*** 0.48*** -
## Harshness 0.59*** 0.38*** 0.59*** 0.69***
## Personality 0.47*** 0.29*** 0.62*** 0.53***
## Harshness Personality
## Punitiveness 0.59*** 0.47***
## Crime Concerns 0.38*** 0.29***
## Process Violations 0.59*** 0.62***
## Emotions 0.69*** 0.53***
## Harshness - 0.67***
## Personality 0.67*** -cat("\nNote: * p < .05, ** p < .01, *** p < .001\n\n")##
## Note: * p < .05, ** p < .01, *** p < .001# --- Construct-Level Table ---
cat("=== CONSTRUCT-LEVEL CORRELATION TABLE ===\n\n")## === CONSTRUCT-LEVEL CORRELATION TABLE ===cat("(See exported CSV for full table)\n\n")## (See exported CSV for full table)# Export full correlation matrix to CSV
construct_cor_export <- round(construct_cor$r, 2)
write.csv(construct_cor_export, "correlation_table_construct_level.csv")
cat("Saved: correlation_table_construct_level.csv\n")## Saved: correlation_table_construct_level.csv# Export p-values
write.csv(round(construct_cor$p, 4), "pvalues_construct_level.csv")
cat("Saved: pvalues_construct_level.csv\n\n")## Saved: pvalues_construct_level.csv# --- Intercorrelation Table ---
cat("=== INTERCORRELATION TABLE (Correlates Only) ===\n\n")## === INTERCORRELATION TABLE (Correlates Only) ===correlate_table <- format_cor_table(correlate_cor$r, correlate_cor$p)
write.csv(correlate_table, "correlation_table_intercorrelations.csv")
cat("Saved: correlation_table_intercorrelations.csv\n\n")## Saved: correlation_table_intercorrelations.csvPHASE 8: SENSITIVITY ANALYSES
8.1 Vignette Effects on Sentence
cat("=== VIGNETTE EFFECTS ON SENTENCE ===\n\n")## === VIGNETTE EFFECTS ON SENTENCE ===# One-way ANOVA
anova_model <- aov(Sentence_1 ~ factor(vignette), data = df_clean)
cat("One-way ANOVA: Sentence ~ Vignette\n")## One-way ANOVA: Sentence ~ Vignetteprint(summary(anova_model))## Df Sum Sq Mean Sq F value Pr(>F)
## factor(vignette) 2 1492 746 4.84 0.0083 **
## Residuals 493 75951 154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Effect size
eta_sq <- effectsize::eta_squared(anova_model)
cat("\nEffect size:\n")##
## Effect size:print(as.data.frame(eta_sq))## Parameter Eta2 CI CI_low CI_high
## 1 factor(vignette) 0.01926 0.95 0.002889 1cat("\n")# Post-hoc tests
cat("Tukey HSD Post-hoc Comparisons:\n")## Tukey HSD Post-hoc Comparisons:print(TukeyHSD(anova_model))## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sentence_1 ~ factor(vignette), data = df_clean)
##
## $`factor(vignette)`
## diff lwr upr p adj
## 2-1 3.050 -0.09773 6.1968 0.0599
## 3-1 -1.003 -4.26915 2.2635 0.7507
## 3-2 -4.052 -7.28321 -0.8215 0.0094cat("\n")8.2 Correlation Stability Across Vignettes
cat("=== CORRELATION STABILITY ACROSS VIGNETTES ===\n\n")## === CORRELATION STABILITY ACROSS VIGNETTES ===# Function to compute punitiveness correlations for each vignette
compute_vignette_cors <- function(data, vignette_num) {
subset_data <- data %>% filter(vignette == vignette_num)
cors <- sapply(correlate_vars, function(var) {
cor(subset_data$punitiveness_agg, subset_data[[var]], use = "pairwise.complete.obs")
})
return(cors)
}
# Compute for each vignette
vignette_cors <- data.frame(
Construct = construct_labels[correlate_vars],
Vignette_1 = compute_vignette_cors(df_clean, 1),
Vignette_2 = compute_vignette_cors(df_clean, 2),
Vignette_3 = compute_vignette_cors(df_clean, 3)
)
# Add overall correlation
vignette_cors$Overall <- sapply(correlate_vars, function(var) {
cor(df_clean$punitiveness_agg, df_clean[[var]], use = "pairwise.complete.obs")
})
# Round for display
vignette_cors <- vignette_cors %>%
mutate(across(where(is.numeric), ~round(., 2)))
cat("Punitiveness Correlations by Vignette:\n")## Punitiveness Correlations by Vignette:print(as.data.frame(vignette_cors), row.names = FALSE)## Construct Vignette_1 Vignette_2 Vignette_3 Overall
## Crime Rates 0.41 0.42 0.36 0.40
## Fear 0.20 0.30 0.03 0.18
## Due Process 0.45 0.45 0.39 0.43
## Uncertain Evidence 0.29 0.37 0.42 0.35
## Hatred 0.60 0.49 0.49 0.53
## Anger 0.55 0.39 0.38 0.45
## Exclusion 0.64 0.46 0.51 0.54
## Degradation 0.52 0.41 0.38 0.44
## Suffering 0.59 0.51 0.43 0.52
## Prison Violence 0.31 0.39 0.35 0.35
## Harsh Conditions 0.62 0.54 0.49 0.56
## Revenge 0.51 0.40 0.37 0.43
## RWA 0.60 0.43 0.46 0.50
## SDO 0.36 0.32 0.28 0.32
## Vengefulness 0.21 0.05 0.09 0.12
## Violence Proneness 0.47 0.33 0.43 0.41
## Racial Resentment 0.51 0.40 0.49 0.46
## Blood Sports 0.12 0.10 0.12 0.11# Compute range across vignettes
vignette_cors$Range <- apply(vignette_cors[, c("Vignette_1", "Vignette_2", "Vignette_3")],
1, function(x) max(x) - min(x))
cat("\n\nLargest Range (least stable correlations):\n")##
##
## Largest Range (least stable correlations):print(as.data.frame(vignette_cors %>% arrange(desc(Range)) %>% head(5) %>%
select(Construct, Vignette_1, Vignette_2, Vignette_3, Range)), row.names = FALSE)## Construct Vignette_1 Vignette_2 Vignette_3 Range
## Fear 0.20 0.30 0.03 0.27
## Exclusion 0.64 0.46 0.51 0.18
## Anger 0.55 0.39 0.38 0.17
## RWA 0.60 0.43 0.46 0.17
## Suffering 0.59 0.51 0.43 0.16cat("\n")8.3 Mixed-Effects Model
cat("=== MIXED-EFFECTS MODEL (Vignette as Random Intercept) ===\n\n")## === MIXED-EFFECTS MODEL (Vignette as Random Intercept) ===# Test key relationship with random intercept for vignette
# Predicting punitiveness from harshness with vignette random effect
mixed_model <- lmer(punitiveness_agg ~ hostile_agg + (1|vignette), data = df_clean)
cat("Model: punitiveness_agg ~ hostile_agg + (1|vignette)\n\n")## Model: punitiveness_agg ~ hostile_agg + (1|vignette)print(summary(mixed_model))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: punitiveness_agg ~ hostile_agg + (1 | vignette)
## Data: df_clean
##
## REML criterion at convergence: 1008
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9325 -0.6919 0.0299 0.7368 2.6031
##
## Random effects:
## Groups Name Variance Std.Dev.
## vignette (Intercept) 0.000 0.000
## Residual 0.439 0.662
## Number of obs: 496, groups: vignette, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.4958 0.0977 494.0000 -15.3 <0.0000000000000002 ***
## hostile_agg 0.3967 0.0247 494.0000 16.1 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## hostile_agg -0.953
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')# Compare to model without random effect
fixed_model <- lm(punitiveness_agg ~ hostile_agg, data = df_clean)
cat("\n\nComparison: Fixed vs Random Effect Models\n")##
##
## Comparison: Fixed vs Random Effect Modelscat(sprintf("Fixed model R²: %.3f\n", summary(fixed_model)$r.squared))## Fixed model R²: 0.343# ICC - with error handling for singularity
icc_result <- tryCatch({
performance::icc(mixed_model)
}, warning = function(w) NULL, error = function(e) NULL)
if(is.null(icc_result) || is.atomic(icc_result)) {
cat("ICC: Near zero (vignette accounts for negligible variance)\n")
cat("(Model showed singularity - random effect variance ≈ 0)\n\n")
} else if(is.list(icc_result) && !is.null(icc_result$ICC_adjusted)) {
cat(sprintf("ICC (Vignette): %.4f\n", icc_result$ICC_adjusted))
cat("(Low ICC suggests vignette accounts for minimal variance)\n\n")
} else {
cat("ICC: Could not be computed (likely near zero)\n\n")
}## ICC: Near zero (vignette accounts for negligible variance)
## (Model showed singularity - random effect variance ≈ 0)PHASE 9: EXPLORATORY ANALYSES
9.1 Process Violations (Dropped Cluster)
cat("=== PROCESS VIOLATIONS (Exploratory) ===\n\n")## === PROCESS VIOLATIONS (Exploratory) ===# Correlations with punitiveness
process_cor_due <- cor.test(df_clean$punitiveness_agg, df_clean$dueprocess_comp)
process_cor_unc <- cor.test(df_clean$punitiveness_agg, df_clean$uncertain_comp)
process_cor_agg <- cor.test(df_clean$punitiveness_agg, df_clean$process_agg)
cat("Punitiveness correlations with Process Violations:\n")## Punitiveness correlations with Process Violations:cat(sprintf(" Due Process Violations: r = %.2f, p %s\n",
process_cor_due$estimate,
ifelse(process_cor_due$p.value < .001, "< .001",
sprintf("= %.3f", process_cor_due$p.value))))## Due Process Violations: r = 0.43, p < .001cat(sprintf(" Uncertain Evidence: r = %.2f, p %s\n",
process_cor_unc$estimate,
ifelse(process_cor_unc$p.value < .001, "< .001",
sprintf("= %.3f", process_cor_unc$p.value))))## Uncertain Evidence: r = 0.35, p < .001cat(sprintf(" Process Violations Aggregate: r = %.2f, p %s\n\n",
process_cor_agg$estimate,
ifelse(process_cor_agg$p.value < .001, "< .001",
sprintf("= %.3f", process_cor_agg$p.value))))## Process Violations Aggregate: r = 0.44, p < .001# Correlations with other clusters
cat("Process Violations correlations with other clusters:\n")## Process Violations correlations with other clusters:process_cluster_vars <- c("crime_concerns_agg", "emotions_agg", "hostile_agg", "personality_agg")
process_cluster_labels_map <- c("crime_concerns_agg" = "Crime Concerns", "emotions_agg" = "Emotions",
"hostile_agg" = "Harshness", "personality_agg" = "Personality")
for(var in process_cluster_vars) {
r <- cor(df_clean$process_agg, df_clean[[var]], use = "pairwise.complete.obs")
cat(sprintf(" %s: r = %.2f\n", process_cluster_labels_map[var], r))
}## Crime Concerns: r = 0.38
## Emotions: r = 0.48
## Harshness: r = 0.59
## Personality: r = 0.62cat("\n")9.2 Blood Sports Analysis
cat("=== BLOOD SPORTS ANALYSIS ===\n\n")## === BLOOD SPORTS ANALYSIS ===# Individual sport correlations with punitiveness
tv_sports <- c("TV_4", "TV_6", "TV_8", "TV_9")
tv_labels <- c("Boxing", "MMA/UFC", "Hunting", "Wrestling")
cat("Individual Blood Sports correlations with Punitiveness:\n")## Individual Blood Sports correlations with Punitiveness:for(i in 1:length(tv_sports)) {
r <- cor(df_clean$punitiveness_agg, df_clean[[tv_sports[i]]],
use = "pairwise.complete.obs")
test <- cor.test(df_clean$punitiveness_agg, df_clean[[tv_sports[i]]])
cat(sprintf(" %s: r = %.2f, p %s\n",
tv_labels[i], r,
ifelse(test$p.value < .001, "< .001", sprintf("= %.3f", test$p.value))))
}## Boxing: r = 0.05, p = 0.247
## MMA/UFC: r = 0.11, p = 0.017
## Hunting: r = 0.15, p < .001
## Wrestling: r = 0.06, p = 0.159cat(sprintf("\n Blood Sports Composite: r = %.2f\n",
cor(df_clean$punitiveness_agg, df_clean$bloodsports_comp,
use = "pairwise.complete.obs")))##
## Blood Sports Composite: r = 0.11cat("\n")9.3 Revenge vs. Vengefulness Distinction
cat("=== REVENGE (Beliefs) vs. VENGEFULNESS (Personality Trait) ===\n\n")## === REVENGE (Beliefs) vs. VENGEFULNESS (Personality Trait) ===# Correlation between revenge and vengefulness
rev_veng_test <- cor.test(df_clean$revenge_comp, df_clean$venge_comp)
cat(sprintf("Correlation: r = %.3f, p %s\n",
rev_veng_test$estimate,
ifelse(rev_veng_test$p.value < .001, "< .001",
sprintf("= %.3f", rev_veng_test$p.value))))## Correlation: r = 0.458, p < .001cat(sprintf("95%% CI: [%.3f, %.3f]\n\n",
rev_veng_test$conf.int[1], rev_veng_test$conf.int[2]))## 95% CI: [0.386, 0.525]# Both with punitiveness for comparison
rev_pun_test <- cor.test(df_clean$punitiveness_agg, df_clean$revenge_comp)
veng_pun_test <- cor.test(df_clean$punitiveness_agg, df_clean$venge_comp)
cat("Correlations with Punitiveness:\n")## Correlations with Punitiveness:cat(sprintf(" Revenge (beliefs/attitudes): r = %.3f, p %s\n",
rev_pun_test$estimate,
ifelse(rev_pun_test$p.value < .001, "< .001",
sprintf("= %.3f", rev_pun_test$p.value))))## Revenge (beliefs/attitudes): r = 0.430, p < .001cat(sprintf(" Vengefulness (personality): r = %.3f, p %s\n\n",
veng_pun_test$estimate,
ifelse(veng_pun_test$p.value < .001, "< .001",
sprintf("= %.3f", veng_pun_test$p.value))))## Vengefulness (personality): r = 0.119, p = 0.008# Steiger's Z comparing the two
r_rev_veng <- cor(df_clean$revenge_comp, df_clean$venge_comp, use = "pairwise.complete.obs")
n_rv <- sum(complete.cases(df_clean[, c("punitiveness_agg", "revenge_comp", "venge_comp")]))
steiger_rv <- cocor.dep.groups.overlap(
r.jk = rev_pun_test$estimate,
r.jh = veng_pun_test$estimate,
r.kh = r_rev_veng,
n = n_rv,
test = "steiger1980"
)
cat("Steiger's Z-test (Revenge vs Vengefulness predicting Punitiveness):\n")## Steiger's Z-test (Revenge vs Vengefulness predicting Punitiveness):cat(sprintf(" Difference = %.3f, p = %.4f\n",
rev_pun_test$estimate - veng_pun_test$estimate,
steiger_rv@steiger1980$p.value))## Difference = 0.311, p = 0.0000cat(sprintf(" Revenge significantly stronger: %s\n\n",
ifelse(steiger_rv@steiger1980$p.value < .05 &
rev_pun_test$estimate > veng_pun_test$estimate, "YES", "NO")))## Revenge significantly stronger: YES9.4 Exploratory Factor Analysis
cat("=== EXPLORATORY FACTOR ANALYSIS ===\n\n")## === EXPLORATORY FACTOR ANALYSIS ===# Select all correlate items for EFA
# Using construct composites for cleaner interpretation
efa_vars <- correlate_vars
efa_data <- df_clean %>% select(all_of(efa_vars)) %>% na.omit()
# Check sample size adequacy
cat(sprintf("Sample size for EFA: N = %d\n", nrow(efa_data)))## Sample size for EFA: N = 496cat(sprintf("Variables: k = %d\n", ncol(efa_data)))## Variables: k = 18cat(sprintf("Subjects per variable: %.1f\n\n", nrow(efa_data)/ncol(efa_data)))## Subjects per variable: 27.6# KMO and Bartlett's test
kmo_result <- KMO(efa_data)
cat(sprintf("KMO Measure of Sampling Adequacy: %.3f\n", kmo_result$MSA))## KMO Measure of Sampling Adequacy: 0.910cat("(>.90 = marvelous, >.80 = meritorious, >.70 = middling, >.60 = mediocre)\n\n")## (>.90 = marvelous, >.80 = meritorious, >.70 = middling, >.60 = mediocre)# Parallel analysis to determine number of factors
cat("Parallel Analysis (determining number of factors):\n")## Parallel Analysis (determining number of factors):parallel_result <- fa.parallel(efa_data, fa = "fa", n.iter = 100, show.legend = FALSE)
## Parallel analysis suggests that the number of factors = 5 and the number of components = NAcat("\n")# Extract suggested number of factors
n_factors <- parallel_result$nfact
# Run EFA
cat(sprintf("Running EFA with %d factors...\n\n", n_factors))## Running EFA with 5 factors...efa_result <- fa(efa_data, nfactors = n_factors, rotate = "varimax", fm = "ml")
# Print loadings
cat("Factor Loadings (sorted by factor):\n")## Factor Loadings (sorted by factor):print(efa_result$loadings, cutoff = 0.3, sort = TRUE)##
## Loadings:
## ML2 ML1 ML3 ML4 ML5
## exclusion_comp 0.719 0.357
## degradation_comp 0.760
## suffering_comp 0.645 0.308
## prisonvi_comp 0.513 0.369
## harsh_comp 0.658 0.304 0.340
## revenge_comp 0.607 0.377
## hatred_comp 0.510 0.572
## anger_comp 0.955
## sdo_comp 0.704
## vprone_comp 0.505 0.397
## raceresent_comp 0.775
## rwa_comp 0.342 0.386 0.613
## bloodsports_comp 0.565
## venge_comp 0.623
## crime_rates_comp 0.475
## fear_comp
## dueprocess_comp 0.393 0.452
## uncertain_comp 0.309 0.389
##
## ML2 ML1 ML3 ML4 ML5
## SS loadings 3.505 2.236 2.187 1.663 0.935
## Proportion Var 0.195 0.124 0.122 0.092 0.052
## Cumulative Var 0.195 0.319 0.440 0.533 0.585# Variance explained
cat("\nVariance Explained:\n")##
## Variance Explained:print(efa_result$Vaccounted)## ML2 ML1 ML3 ML4 ML5
## SS loadings 3.5055 2.2363 2.1871 1.66341 0.93482
## Proportion Var 0.1947 0.1242 0.1215 0.09241 0.05193
## Cumulative Var 0.1947 0.3190 0.4405 0.53291 0.58484
## Proportion Explained 0.3330 0.2124 0.2078 0.15801 0.08880
## Cumulative Proportion 0.3330 0.5454 0.7532 0.91120 1.00000cat("\n")# Save EFA results
png("efa_scree_plot.png", width = 800, height = 600, res = 150)
plot(efa_result$values, type = "b", main = "Scree Plot",
xlab = "Factor", ylab = "Eigenvalue")
abline(h = 1, lty = 2, col = "red")
dev.off()## quartz_off_screen
## 2cat("Saved: efa_scree_plot.png\n\n")## Saved: efa_scree_plot.png9.5 Confirmatory Factor Analysis (Pre-registered)
cat("=== CONFIRMATORY FACTOR ANALYSIS ===\n\n")## === CONFIRMATORY FACTOR ANALYSIS ===# Theoretical cluster model (confirmatory + Process exploratory)
cfa_model <- '
CrimeConcerns =~ crime_rates_comp + fear_comp
ProcessViolations =~ dueprocess_comp + uncertain_comp
Emotions =~ hatred_comp + anger_comp
Harshness =~ exclusion_comp + degradation_comp + suffering_comp +
prisonvi_comp + harsh_comp + revenge_comp
Personality =~ rwa_comp + sdo_comp + venge_comp + vprone_comp +
raceresent_comp + bloodsports_comp
'
cfa_fit <- tryCatch({
cfa(cfa_model, data = df_clean, std.lv = TRUE)
}, error = function(e) {
cat("CFA model did not converge. Trying with different estimator...\n")
cfa(cfa_model, data = df_clean, std.lv = TRUE, estimator = "MLR")
})
cat("CFA Model Fit Indices:\n")## CFA Model Fit Indices:fit_indices <- fitMeasures(cfa_fit, c("cfi", "tli", "rmsea", "rmsea.ci.lower",
"rmsea.ci.upper", "srmr", "chisq", "df", "pvalue"))
print(round(fit_indices, 3))## cfi tli rmsea rmsea.ci.lower rmsea.ci.upper
## 0.858 0.826 0.105 0.098 0.112
## srmr chisq df pvalue
## 0.064 805.808 125.000 0.000cat("\nInterpretation guidelines:\n")##
## Interpretation guidelines:cat(" CFI/TLI > .90 acceptable, > .95 good\n")## CFI/TLI > .90 acceptable, > .95 goodcat(" RMSEA < .08 acceptable, < .05 good\n")## RMSEA < .08 acceptable, < .05 goodcat(" SRMR < .08 acceptable\n\n")## SRMR < .08 acceptable# Factor loadings
cat("Standardized Factor Loadings:\n")## Standardized Factor Loadings:std_solution <- standardizedSolution(cfa_fit) %>%
filter(op == "=~") %>%
select(Factor = lhs, Indicator = rhs, Loading = est.std, p = pvalue) %>%
mutate(Loading = round(Loading, 3), p = round(p, 4))
print(as.data.frame(std_solution), row.names = FALSE)## Factor Indicator Loading p
## CrimeConcerns crime_rates_comp 0.739 0
## CrimeConcerns fear_comp 0.531 0
## ProcessViolations dueprocess_comp 0.814 0
## ProcessViolations uncertain_comp 0.680 0
## Emotions hatred_comp 0.910 0
## Emotions anger_comp 0.788 0
## Harshness exclusion_comp 0.809 0
## Harshness degradation_comp 0.830 0
## Harshness suffering_comp 0.831 0
## Harshness prisonvi_comp 0.554 0
## Harshness harsh_comp 0.843 0
## Harshness revenge_comp 0.745 0
## Personality rwa_comp 0.724 0
## Personality sdo_comp 0.674 0
## Personality venge_comp 0.380 0
## Personality vprone_comp 0.774 0
## Personality raceresent_comp 0.683 0
## Personality bloodsports_comp 0.389 0cat("\n")# Factor correlations
cat("Factor Correlations:\n")## Factor Correlations:factor_cors <- standardizedSolution(cfa_fit) %>%
filter(op == "~~", lhs != rhs) %>%
select(Factor1 = lhs, Factor2 = rhs, r = est.std) %>%
mutate(r = round(r, 3))
print(as.data.frame(factor_cors), row.names = FALSE)## Factor1 Factor2 r
## CrimeConcerns ProcessViolations 0.589
## CrimeConcerns Emotions 0.611
## CrimeConcerns Harshness 0.540
## CrimeConcerns Personality 0.514
## ProcessViolations Emotions 0.626
## ProcessViolations Harshness 0.755
## ProcessViolations Personality 0.794
## Emotions Harshness 0.802
## Emotions Personality 0.649
## Harshness Personality 0.786cat("\n")9.6 Demographic Moderators
cat("=== DEMOGRAPHIC MODERATORS (Exploratory) ===\n\n")## === DEMOGRAPHIC MODERATORS (Exploratory) ===# Political Orientation
cat("POLITICAL ORIENTATION:\n")## POLITICAL ORIENTATION:# Split into liberal (1-3), moderate (4), conservative (5-7)
df_clean <- df_clean %>%
mutate(political_group = case_when(
politid <= 3 ~ "Liberal",
politid == 4 ~ "Moderate",
politid >= 5 ~ "Conservative"
))
# Correlations by political group
political_cors <- df_clean %>%
group_by(political_group) %>%
summarise(
n = n(),
r_hostile = cor(punitiveness_agg, hostile_agg, use = "pairwise.complete.obs"),
r_crime = cor(punitiveness_agg, crime_concerns_agg, use = "pairwise.complete.obs"),
r_personality = cor(punitiveness_agg, personality_agg, use = "pairwise.complete.obs")
) %>%
mutate(across(starts_with("r_"), ~round(., 2)))
print(as.data.frame(political_cors))## political_group n r_hostile r_crime r_personality
## 1 Conservative 190 0.56 0.20 0.39
## 2 Liberal 190 0.57 0.35 0.41
## 3 Moderate 116 0.47 0.31 0.34cat("\n")# Gender
cat("GENDER:\n")## GENDER:gender_cors <- df_clean %>%
filter(gender %in% c(1, 2)) %>%
mutate(gender_label = ifelse(gender == 1, "Male", "Female")) %>%
group_by(gender_label) %>%
summarise(
n = n(),
mean_punitiveness = mean(punitiveness_agg, na.rm = TRUE),
r_hostile = cor(punitiveness_agg, hostile_agg, use = "pairwise.complete.obs")
) %>%
mutate(across(where(is.numeric) & !n, ~round(., 2)))
print(as.data.frame(gender_cors))## gender_label n mean_punitiveness r_hostile
## 1 Female 249 0.09 0.59
## 2 Male 244 -0.10 0.60# T-test for gender difference in punitiveness
gender_ttest <- t.test(punitiveness_agg ~ gender,
data = df_clean %>% filter(gender %in% c(1, 2)))
cat(sprintf("\nGender difference in punitiveness: t = %.2f, p = %.3f\n",
gender_ttest$statistic, gender_ttest$p.value))##
## Gender difference in punitiveness: t = -2.71, p = 0.007# Effect size for gender difference
gender_d <- effectsize::cohens_d(punitiveness_agg ~ gender,
data = df_clean %>% filter(gender %in% c(1, 2)))
cat(sprintf("Cohen's d = %.2f\n\n", gender_d$Cohens_d))## Cohen's d = -0.24PHASE 10: MULTIPLE REGRESSION
10.1 All Clusters Predicting Punitiveness
cat("=== ALL CLUSTERS PREDICTING PUNITIVENESS ===\n\n")## === ALL CLUSTERS PREDICTING PUNITIVENESS ===# Model with all cluster aggregates
cluster_reg <- lm(punitiveness_agg ~ crime_concerns_agg + process_agg + emotions_agg +
hostile_agg + personality_agg, data = df_clean)
cat("Model: Punitiveness ~ Crime Concerns + Process Violations + Emotions + Harshness + Personality\n\n")## Model: Punitiveness ~ Crime Concerns + Process Violations + Emotions + Harshness + Personalityprint(summary(cluster_reg))##
## Call:
## lm(formula = punitiveness_agg ~ crime_concerns_agg + process_agg +
## emotions_agg + hostile_agg + personality_agg, data = df_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8230 -0.4588 0.0146 0.4993 1.4830
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.9303 0.1241 -15.55 < 0.0000000000000002 ***
## crime_concerns_agg 0.0406 0.0285 1.43 0.154
## process_agg 0.0593 0.0389 1.53 0.128
## emotions_agg 0.1246 0.0306 4.08 0.0000534 ***
## hostile_agg 0.2127 0.0394 5.40 0.0000001 ***
## personality_agg 0.0717 0.0398 1.80 0.072 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.641 on 490 degrees of freedom
## Multiple R-squared: 0.39, Adjusted R-squared: 0.383
## F-statistic: 62.6 on 5 and 490 DF, p-value: <0.0000000000000002# Standardized coefficients
cat("\nStandardized Coefficients (Beta):\n")##
## Standardized Coefficients (Beta):cluster_reg_std <- lm(scale(punitiveness_agg) ~ scale(crime_concerns_agg) +
scale(process_agg) + scale(emotions_agg) +
scale(hostile_agg) + scale(personality_agg), data = df_clean)
print(round(coef(cluster_reg_std)[-1], 3))## scale(crime_concerns_agg) scale(process_agg) scale(emotions_agg)
## 0.057 0.073 0.208
## scale(hostile_agg) scale(personality_agg)
## 0.314 0.092cat("\n")# Unique variance (squared semi-partial correlations)
cat("Unique Variance Explained (sr²):\n")## Unique Variance Explained (sr²):# Using car package for Type III SS
anova_type3 <- Anova(cluster_reg, type = 3)
ss_total <- sum(anova_type3$`Sum Sq`)
unique_var <- anova_type3$`Sum Sq` / ss_total
names(unique_var) <- rownames(anova_type3)
print(round(unique_var[-1], 4)) # Exclude intercept## crime_concerns_agg process_agg emotions_agg hostile_agg
## 0.0026 0.0030 0.0212 0.0371
## personality_agg Residuals
## 0.0041 0.6239cat("\n")10.2 All Constructs Predicting Punitiveness
cat("=== ALL CONSTRUCTS PREDICTING PUNITIVENESS ===\n\n")## === ALL CONSTRUCTS PREDICTING PUNITIVENESS ===construct_reg <- lm(punitiveness_agg ~ crime_rates_comp + fear_comp +
dueprocess_comp + uncertain_comp +
hatred_comp + anger_comp +
exclusion_comp + degradation_comp + suffering_comp +
prisonvi_comp + harsh_comp + revenge_comp +
rwa_comp + sdo_comp + venge_comp + vprone_comp +
raceresent_comp + bloodsports_comp,
data = df_clean)
cat("Full Model Summary:\n")## Full Model Summary:print(summary(construct_reg))##
## Call:
## lm(formula = punitiveness_agg ~ crime_rates_comp + fear_comp +
## dueprocess_comp + uncertain_comp + hatred_comp + anger_comp +
## exclusion_comp + degradation_comp + suffering_comp + prisonvi_comp +
## harsh_comp + revenge_comp + rwa_comp + sdo_comp + venge_comp +
## vprone_comp + raceresent_comp + bloodsports_comp, data = df_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5987 -0.4074 -0.0139 0.4384 1.5275
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.82108 0.14046 -12.96 < 0.0000000000000002 ***
## crime_rates_comp 0.07249 0.02536 2.86 0.00444 **
## fear_comp -0.00734 0.02408 -0.30 0.76069
## dueprocess_comp 0.06154 0.03323 1.85 0.06462 .
## uncertain_comp 0.00705 0.03274 0.22 0.82968
## hatred_comp 0.10116 0.03345 3.02 0.00263 **
## anger_comp -0.01882 0.03046 -0.62 0.53709
## exclusion_comp 0.08399 0.03152 2.66 0.00797 **
## degradation_comp -0.06156 0.03441 -1.79 0.07420 .
## suffering_comp 0.06728 0.02808 2.40 0.01697 *
## prisonvi_comp 0.04238 0.02623 1.62 0.10682
## harsh_comp 0.05388 0.03171 1.70 0.08995 .
## revenge_comp 0.01848 0.02994 0.62 0.53732
## rwa_comp 0.05236 0.02600 2.01 0.04460 *
## sdo_comp -0.06853 0.02657 -2.58 0.01020 *
## venge_comp -0.04976 0.02351 -2.12 0.03480 *
## vprone_comp 0.03496 0.02843 1.23 0.21945
## raceresent_comp 0.08190 0.02262 3.62 0.00032 ***
## bloodsports_comp -0.04353 0.02025 -2.15 0.03204 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.604 on 477 degrees of freedom
## Multiple R-squared: 0.472, Adjusted R-squared: 0.452
## F-statistic: 23.7 on 18 and 477 DF, p-value: <0.0000000000000002# Top predictors by standardized coefficient
cat("\nSignificant Predictors (p < .05):\n")##
## Significant Predictors (p < .05):construct_summary <- tidy(construct_reg) %>%
filter(p.value < .05, term != "(Intercept)") %>%
arrange(desc(abs(estimate)))
print(as.data.frame(construct_summary))## term estimate std.error statistic p.value
## 1 hatred_comp 0.10116 0.03345 3.024 0.0026294
## 2 exclusion_comp 0.08399 0.03152 2.665 0.0079668
## 3 raceresent_comp 0.08190 0.02262 3.622 0.0003243
## 4 crime_rates_comp 0.07249 0.02536 2.859 0.0044371
## 5 sdo_comp -0.06853 0.02657 -2.579 0.0102005
## 6 suffering_comp 0.06728 0.02808 2.396 0.0169695
## 7 rwa_comp 0.05236 0.02600 2.014 0.0446022
## 8 venge_comp -0.04976 0.02351 -2.117 0.0347963
## 9 bloodsports_comp -0.04353 0.02025 -2.150 0.0320437cat("\n")10.3 Hierarchical Regression (Incremental R²)
cat("=== HIERARCHICAL REGRESSION ===\n\n")## === HIERARCHICAL REGRESSION ===# Step 1: Crime Concerns only
step1 <- lm(punitiveness_agg ~ crime_concerns_agg, data = df_clean)
# Step 2: Add Process Violations
step2 <- lm(punitiveness_agg ~ crime_concerns_agg + process_agg, data = df_clean)
# Step 3: Add Harshness
step3 <- lm(punitiveness_agg ~ crime_concerns_agg + process_agg + hostile_agg, data = df_clean)
# Step 4: Add Emotions
step4 <- lm(punitiveness_agg ~ crime_concerns_agg + process_agg + hostile_agg + emotions_agg, data = df_clean)
# Step 5: Add Personality
step5 <- lm(punitiveness_agg ~ crime_concerns_agg + process_agg + hostile_agg + emotions_agg +
personality_agg, data = df_clean)
cat("Hierarchical Regression Results:\n")## Hierarchical Regression Results:cat(sprintf("Step 1 (Crime Concerns): R² = %.3f\n", summary(step1)$r.squared))## Step 1 (Crime Concerns): R² = 0.103cat(sprintf("Step 2 (+Process): R² = %.3f, ΔR² = %.3f\n",
summary(step2)$r.squared,
summary(step2)$r.squared - summary(step1)$r.squared))## Step 2 (+Process): R² = 0.221, ΔR² = 0.117cat(sprintf("Step 3 (+Harshness): R² = %.3f, ΔR² = %.3f\n",
summary(step3)$r.squared,
summary(step3)$r.squared - summary(step2)$r.squared))## Step 3 (+Harshness): R² = 0.363, ΔR² = 0.142cat(sprintf("Step 4 (+Emotions): R² = %.3f, ΔR² = %.3f\n",
summary(step4)$r.squared,
summary(step4)$r.squared - summary(step3)$r.squared))## Step 4 (+Emotions): R² = 0.386, ΔR² = 0.022cat(sprintf("Step 5 (+Personality): R² = %.3f, ΔR² = %.3f\n\n",
summary(step5)$r.squared,
summary(step5)$r.squared - summary(step4)$r.squared))## Step 5 (+Personality): R² = 0.390, ΔR² = 0.004# ANOVA comparing models
cat("Model Comparison (ANOVA):\n")## Model Comparison (ANOVA):print(anova(step1, step2, step3, step4, step5))## Analysis of Variance Table
##
## Model 1: punitiveness_agg ~ crime_concerns_agg
## Model 2: punitiveness_agg ~ crime_concerns_agg + process_agg
## Model 3: punitiveness_agg ~ crime_concerns_agg + process_agg + hostile_agg
## Model 4: punitiveness_agg ~ crime_concerns_agg + process_agg + hostile_agg +
## emotions_agg
## Model 5: punitiveness_agg ~ crime_concerns_agg + process_agg + hostile_agg +
## emotions_agg + personality_agg
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 494 296
## 2 493 257 1 38.7 94.29 < 0.0000000000000002 ***
## 3 492 210 1 47.0 114.27 < 0.0000000000000002 ***
## 4 491 203 1 7.4 17.97 0.000027 ***
## 5 490 201 1 1.3 3.24 0.072 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n")PHASE 11: SAVE OUTPUT
cat("=== SAVING OUTPUT FILES ===\n\n")## === SAVING OUTPUT FILES ===# Save cleaned dataset
write.csv(df_clean, "punishment_212_cleaned_data.csv", row.names = FALSE)
cat("Saved: punishment_212_cleaned_data.csv\n")## Saved: punishment_212_cleaned_data.csv# Save alpha table
write.csv(alpha_table, "reliability_alphas.csv", row.names = FALSE)
cat("Saved: reliability_alphas.csv\n")## Saved: reliability_alphas.csv# Save descriptives
write.csv(as.data.frame(descriptives_table), "descriptive_statistics.csv", row.names = FALSE)
cat("Saved: descriptive_statistics.csv\n")## Saved: descriptive_statistics.csv# Save H1 results
write.csv(as.data.frame(fdr_table), "H1_punitiveness_correlations.csv", row.names = FALSE)
cat("Saved: H1_punitiveness_correlations.csv\n")## Saved: H1_punitiveness_correlations.csv# Save vignette correlation stability
write.csv(as.data.frame(vignette_cors), "vignette_correlation_stability.csv", row.names = FALSE)
cat("Saved: vignette_correlation_stability.csv\n")## Saved: vignette_correlation_stability.csvcat("\n")SUMMARY
cat("\n", paste(rep("=", 70), collapse = ""), "\n")##
## ======================================================================cat("SUMMARY OF KEY FINDINGS\n")## SUMMARY OF KEY FINDINGScat(paste(rep("=", 70), collapse = ""), "\n\n")## ======================================================================cat("SAMPLE:\n")## SAMPLE:cat(sprintf(" Final N = %d (of %d collected)\n", nrow(df_clean), n_raw))## Final N = 496 (of 514 collected)cat(sprintf(" Exclusions: %d total\n\n", n_raw - nrow(df_clean)))## Exclusions: 18 totalcat("HYPOTHESIS TESTS:\n\n")## HYPOTHESIS TESTS:cat("H1: Punitiveness positively correlated with all correlates\n")## H1: Punitiveness positively correlated with all correlatescat(sprintf(" Construct-level: %d/%d (%.1f%%) positive correlations\n",
h1_n_positive, h1_n_total, h1_n_positive/h1_n_total*100))## Construct-level: 18/18 (100.0%) positive correlationscat(sprintf(" Significant after FDR: %d/%d\n", h1_n_sig_fdr, h1_n_total))## Significant after FDR: 16/18cat(sprintf(" SUPPORTED: %s\n\n",
ifelse(h1_n_positive/h1_n_total > 0.9, "YES", "PARTIAL")))## SUPPORTED: YEScat("H2: Harshness > Crime Concerns\n")## H2: Harshness > Crime Concernscat(sprintf(" r(Punitiveness, Harshness) = %.3f\n", r_hostile))## r(Punitiveness, Harshness) = 0.586cat(sprintf(" r(Punitiveness, Crime) = %.3f\n", r_crime))## r(Punitiveness, Crime) = 0.322cat(sprintf(" Difference = %.3f, p = %.4f\n", r_hostile - r_crime, h2_pvalue))## Difference = 0.264, p = 0.0000cat(sprintf(" SUPPORTED: %s\n\n",
ifelse(r_hostile > r_crime & h2_pvalue < .05, "YES (Harshness > Crime)", "NO")))## SUPPORTED: YES (Harshness > Crime)cat("H3: Correlates positively intercorrelated\n")## H3: Correlates positively intercorrelatedcat(sprintf(" Positive correlations: %d/%d (%.1f%%)\n",
h3_n_positive, h3_n_correlations, h3_n_positive/h3_n_correlations*100))## Positive correlations: 153/153 (100.0%)cat(sprintf(" SUPPORTED: %s\n\n",
ifelse(h3_n_positive/h3_n_correlations > 0.5, "YES", "NO")))## SUPPORTED: YEScat("KEY CORRELATES (Top 5 by |r|):\n")## KEY CORRELATES (Top 5 by |r|):print(as.data.frame(head(punitiveness_cors %>% select(Construct, r, sig), 5)), row.names = FALSE)## Construct r sig
## Harsh Conditions 0.56 ***
## Exclusion 0.54 ***
## Hatred 0.53 ***
## Suffering 0.52 ***
## RWA 0.50 ***cat("\n")cat("REGRESSION SUMMARY:\n")## REGRESSION SUMMARY:cat(sprintf(" Total R² (all clusters): %.3f\n", summary(cluster_reg)$r.squared))## Total R² (all clusters): 0.390cat(sprintf(" Harshness unique variance: %.3f\n", unique_var["hostile_agg"]))## Harshness unique variance: 0.037cat("\n")cat(paste(rep("=", 70), collapse = ""), "\n")## ======================================================================cat("ANALYSIS COMPLETE\n")## ANALYSIS COMPLETEcat(paste(rep("=", 70), collapse = ""), "\n")## ======================================================================