---
title: "Predicting health information avoidance using machine learning models"
author: 
  - name: Zihan Hei
    corresponding: false
    orcid: 0009-0007-1617-2666 
    roles:
      - Software
      - Formal analysis
      - Visualization
      - Writing - original draft
    affiliations:
      - Binghamton University
      
editor: visual
execute:
  echo: false
  warning: false
  message: false

format:
  html:
    toc: true
    code-tools: false
  pdf:
    toc: true    
    number-sections: true
    geometry: margin=1in
    fontsize: 11pt  
    include-in-header:
      text: |
        \usepackage{microtype}
        \usepackage{longtable}
        \usepackage{array}
    
bibliography: references.bib
csl: apa.csl
---

# Introduction

Health information avoidance—defined as any behavior designed to prevent or delay access to available but potentially unwanted health information (@howell2020; @sweeny2010)—remains a significant barrier to realizing the public health benefits of personalized risk communication. Despite the rapid expansion of access to medical information in the Information Age, many people still choose not to learn about their personal health risks (@gigerenzer2017; @ho2021; @kelly2021). Health information avoidance is not a uniform behavior, but can take multiple forms, varying in its duration and degree of intentionality. Individuals may avoid health information by delaying the decision to learn about their screening results or they may avoid it completely, choosing to never know. Information avoidance may also manifest through both active and passive means. People may explicitly ask others not to disclose information, physically remove themselves from situations were information might be revealed, or passively refrain from questions that could reveal unwanted knowledge. These avoidance behaviors may manifest as refusing cancer screening, delaying medical care, or not requesting test results (@sweeny2010).

Research indicates that health information avoidance is common. For example, approximately 15% of US adults avoid personalized health risk information across various contexts ([Meese et al., 2022](https://www.sciencedirect.com/science/article/pii/S2667321524000702#bib34)), and nearly 39% express a reluctance to learn about their cancer risk (@emanuel2015). Similar avoidance rates (e.g., approximately 40%) have been observed for non-personalized health information (i.e., general health information not relevant to personal risk) (@chae2019; @orom2020).

Such widespread avoidance underscores the need to understand how and why individuals avoid health information. Some individuals may avoid health information to protect themselves from unpleasant emotions, prevent exposure to information that conflicts with their worldview or creates an obligation to act. Even when the information may be critical to health,avoidance eliminates the discomfort of decision-making and the emotional burden of confronting potential illness (@sweeny2010). Previous research has identified various psychological and cognitive factors underlying this avoidance phenomenon and explored potential explanations for this behavior. @obrien2024 found that self-perceptions of health, such as low perceived risk, engagement in healthy behaviors, and demographic characteristics, often guide people’s decisions to avoid learning about their health risks. Other studies have shown that information overload can increase anxiety and cognitive dissonance, leading to avoidance behavior (@dattilo2022; @song2021; @soroya2023). Furthermore, heightened risk perceptions can exacerbate anxiety and sadness, which may further hinder people from seeking health information (@sultana2023; @zhao2009).

However, few studies have attempted to use machine learning methods to identify predictors of health information avoidance. This study applies predictive modeling and machine learning methods to examine patterns of cancer screening avoidance (“health avoiders”) using sociodemographic and psychological data. By integrating behavioral and belief factors, this study aimed to better understand the complex dynamics that lead to health information avoidance. The key finding is that while many predictors showed statistical significance, none provided meaningful predictive power. This suggests that the key determinants leading people to avoid health information remain unidentified.

## Review of Machine Learning Methods

Machine learning (ML) methods were applied to identify patterns associated with health information avoidance. This study focuses on three main models: **Logistic Regression**, **Random Forest**, and **Multivariate Adaptive Regression Splines (MARS)**. Each approach offers different advantages and trade-offs, summarized in Table 1 below.

### Table 1. Comparison of Machine Learning Models

| **Model** | **Description** | **Advantages** | **Limitations** | **Interpretability** |
|:--------------|:--------------|:--------------|:--------------|:--------------|
| [**Linear Regression**](https://www.mathsjournal.com/pdf/2023/vol8issue6/PartB/8-6-22-868.pdf) | Models the relationship between predictors and a continuous outcome using a straight-line equation. | Simple, interpretable, and efficient to train; effective for linearly related data; allows understanding of variable relationships. | Assumes linearity and independence among variables; sensitive to outliers; limited for complex or nonlinear data. | High |
| [**Logistic Regression**](https://www.geeksforgeeks.org/data-science/advantages-and-disadvantages-of-logistic-regression/) | Estimates the probability of a binary outcome based on a linear combination of predictors. | Simple, interpretable, good baseline; easy to assess predictor importance. | Assumes linearity; struggles with nonlinear or high-dimensional data. | High |
| [**Random Forest**](https://en.wikipedia.org/wiki/Random_forest#Disadvantages) | Ensemble of decision trees using bootstrapped samples and random feature selection. | Handles nonlinear and complex data; reduces overfitting; provides feature importance. | Less interpretable; slower with large datasets; harder to explain model logic. | Moderate |
| [**MARS (Multivariate Adaptive Regression Splines)**](https://bradleyboehmke.github.io/HOML/mars.html#) | Builds flexible regression models using piecewise linear splines. | Captures nonlinear relationships; performs automatic feature selection; minimal preprocessing needed. | Computationally intensive; interpretation can be challenging with correlated predictors. | Moderate - High |

Linear / Logistic Regression served as a baseline interpretable model, Random Forest captured complex nonlinear interactions, and MARS modeled adaptive spline-based relationships between predictors and cancer avoidance. By comparing these models, this study aimed to identify the most effective approach to predict health information avoidance while maintaining interpretability.

# Method

## Data Source

This study used the de-identified Health Avoiders dataset provided through [Cloud Research](https://www.cloudresearch.com/pricing/) in collaboration with Dr. Heather Orom (University at Buffalo). The dataset includes sociodemographic, psychological, and behavioral variables collected by Cloud Research. All analyses were conducted in R, and reproducibility was ensured through a README file and a Quarto documentation workflow.

## Outcome Variable

The outcome variable, `Cancer_Avoidance_Mean`, represents the average score across 8 items (`Avoid_Cancer_`)measuring participants’ avoidance of cancer-related health information. Because this variable exhibited non-normality, both logarithmic and square-root transformations were tested; however, these transformations did not improve model performance or interpretability, so the untransformed variable was retained for analysis.

## Predictor Models

Heather Orom (University at Buffalo) developed a theoretical framework to group predictors as a health psychologist studying health information avoidance.

1.  Demographic Model (`demo_data`) — `Ethnicity`, `Political_Party`, `Gender4`, `Job_Classification`, `Education_Level`, `Age`, `Income`, `Race`, and `MacArthur_Numeric`.

2.  Media Use Model (`media_data`) — `Social_Media_Usage`, `AI_Use`, `Video_Games_Hours`, `Listening_Podcasts`, `Facebook_Usage_cat`, `TikTok_Use`, `X_Twitter_Usage`, `Social_Media_type`, and `Influencer_Following`.

3.  Health Condition Model (`health_condition_data`) — `Stressful_Events_Recent`, `Current_Depression`, `Anxiety_Severity_num`, `PTSD5_Score`, `Health_Depression_Severity_num`, and `Stress_TotalScore`.

4.  Health Behavior Model (`health_behavior_data`) — `Fast_Food_Consumption`, `Meditation_group`, `Physical_Activity_Guidelines`, `Cigarette_Smoking_num`, `Supplement_Consumption_Reason_num`, `Diet_Type`, and `Supplement_Consumption`.

5.  Other Factors Model (`other_data`) — `Home_Ownership`, `Voter_Registration`, `Climate_Change_Belief`, and `Mental_Health_of_Partner`.

## Variable Scaling and Scoring

Several psychological and health condition variables were standardized using validated scales:

-   [PC-PTSD-5](https://www.ptsd.va.gov/professional/assessment/screens/pc-ptsd.asp): a 5-item yes/no screen for post-traumatic stress.

-   [GAD-7](https://adaa.org/sites/default/files/GAD-7_Anxiety-updated_0.pdf): a 7-item measure of anxiety severity (0–3 scale).

-   [PHQ-9](https://www.socialworkportal.com/phq-9-questionnaire/): a 9-item measure of depression severity (0–3 scale).

-   Life Events Checklist: summed to represent cumulative stress exposure.

Composite averages were computed for each domain, resulting in 4 continuous variables: `PTSD5_Score`, `Anxiety_Severity_num`, `Health_Depression_Severity_num`, and `Stress_TotalScore`.

## Analysis Plan

Predictive modeling was conducted in R using the `tidymodels` package. Due to the non-normal outcome variable, two complimentary approaches were applied:

1.  Regression Models: predicted the continuous outcome `Cancer_Avoidance_Mean` were fitted using Linear Regression (baseline), Random Forest, and Multivariate Adaptive Regression Splines (MARS). Model performance was assessed using root mean squared error (RMSE), mean absolute error (MAE), and the correlation between predicted and observed values.

2.  Classification Models: predicted the binary outcome `Cancer_Avoiders01` (0 = non-avoider, 1 = avoider) were fitted using Logistic Regression (baseline) and Random Forest. Model performance was assessed using area under the curve (AUC) and accuracy with basic calibration applied when class imbalance occurred.

Some models used cross-validation, and Linear / Logistic Regression served as the baseline for comparison.

# Results

All analyses were conducted in Posit Cloud, a collaborative online platform for writing R scripts and developing Quarto markdown documents.

## Overview of Model Performance

Table 2 summarizes the predictive performance across all modeling approaches. Consistent with the study's central finding, all models demonstrated **statistically significant associations but no meaningful predictive power**.

### Table 2. Summary of Model Performance Across All Approaches

| **Domain** | **Model Type** | **Key Metric** | **Value** | **Interpretation** |
|:--------------|:--------------|:--------------|:--------------|:--------------|
| Demographics | Regression (RF) | Correlation | 0.109 | No predictive relationship |
| Demographics | Regression (MARS) | Cross-val R² | 0.017 | Poor generalization |
| Demographics | Classification (RF) | ROC AUC | 0.449 | Below chance performance |
| Media Usage | Regression (RF) | Correlation | 0.021 | No predictive relationship |
| Media Usage | Classification (RF) | ROC AUC | 0.574 | Marginally above chance |
| Health Condition | Regression (RF) | Correlation | 0.103 | No predictive relationship |
| Health Condition | Regression (MARS) | Cross-val R² | -0.003 | Negative (no pattern) |
| Health Condition | Classification (RF) | ROC AUC | 0.446 | Below chance performance |
| Health Behavior | Regression (RF) | Correlation | 0.060 | No predictive relationship |
| Health Behavior | Classification (RF) | ROC AUC | 0.472 | Below chance performance |
| Other Factors | Regression (RF) | Correlation | 0.145 | Weak relationship |
| Other Factors | Classification (RF) | ROC AUC | 0.420 | Below chance performance |
| **Full Model** | **Regression (RF)** | **Correlation** | **0.300** | **Weak relationship** |
| **Full Model** | **Classification (RF)** | **ROC AUC** | **0.230** | **Worse than individual models** |

*Note: RF = Random Forest; MARS = Multivariate Adaptive Regression Splines. Classification models showed high accuracy (94-96%) due to class imbalance but uniformly poor discrimination (AUC \< 0.60). All R² values were \< 0.03.*

The following sections detail findings for each predictor domain, focusing on the strongest associations identified despite their limited predictive utility.

## Overview of Findings

This analysis examined five predictor domains: demographic, health status, health behaviors, media usage, and other attitudinal variables. In all modeling approaches, predictors were statistically significant with cancer avoidance behavior, but their predictive power was limited. The following paragraphs detail the findings for each domain.

### Demographic Predictors

Linear regression analysis showed a significant negative association between socioeconomic status (MacArthur numerical score) and cancer avoidance behavior (β = -0.025, p = 1.89e-10), but the model explained only 0.5% of the variance (R² = 0.005). Age analysis identified a threshold of 35 years, with older individuals exhibiting significantly higher cancer avoidance scores (β = 0.063, p = 6.45e-06).

In both age groups (under 35 and over 35), Democrats had lower cancer avoidance scores than Republicans: under 35 (β = -0.165, p = 8.97e-08) and over 35 (β = -0.225, p = 2.67e-15). MARS analysis confirmed that party affiliation was the most important factor (importance = 100). Although the classification model showed high accuracy (95.5%), its discrimination was poor (ROC AUC = 0.449), indicating statistically significant but not practically predictive effects.

### Health Condition Predictors

Predictive relationships among health condition variables were weak. The total stress score showed no significant trend (β = 0.005, p = 0.073, R² = 0.0005). While the MARS model identified anxiety severity as the most important factor, cross-validation performance was negative (CVRSq = -0.003), suggesting no generalizable pattern. The classification model also failed to meaningfully differentiate cancer avoiders.

### Health Behavior Predictors

Cigarette smoking demonstrated a significant positive relationship with cancer avoidance. Smokers had higher cancer avoidance scores (β = 0.186, p < 2e-16, R² = 0.011) and 1.87 times higher odds of being cancer avoiders (β = 0.627, p = 1.05e-06). Approximately 8–9% of smokers were cancer avoiders compared to 4–5% of non-smokers. Despite statistical significance, the overall predictive power remained limited.

### Media Usage Predictors

Media usage variables showed minimal predictive value (overall correlation = 0.021). Influencer following was not significantly associated with cancer avoidance (β = -0.011, p = 0.613, R² = 0.00008). Facebook usage showed a statistically significant but weak effect (β = 0.236, p = 0.039), with nearly equal proportions of cancer avoiders among users (5%) and non-users (4%). These effects are statistically significant but not practically meaningful.

### Other Predictors: Climate Change Beliefs

Climate change beliefs were statistically associated with cancer avoidance (R² = 0.024). Individuals who strongly believe in human-caused climate change had lower cancer avoidance scores (β = -0.271, p = 3.32e-13) and 0.29 times the odds of being cancer avoiders compared to climate deniers (β = -1.236, p = 1.87e-10). Climate deniers showed 12–13% cancer avoiders compared to 3–4% among strong believers. Despite statistical significance, model differentiation remained poor (ROC AUC = 0.42).

### Comprehensive Model: Combining Major Predictors

Random forest models incorporating all main predictors showed modest improvement. The regression model achieved a correlation of 0.30 (RMSE = 0.62, MAE = 0.52). The classification model achieved 95% accuracy but failed to meaningfully differentiate cancer avoiders (ROC AUC = 0.23).

Across all analyses, predictors were statistically significant but not practically predictive (R² < 0.03, ROC AUC < 0.50), indicating that the measured variables have limited explanatory power for cancer avoidance behavior.

```{r}
#| label: setup
#| include: false
#| echo: false

# Set knitr options for PDF compatibility
knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE
)

# If rendering to PDF, use simple table format
if (knitr::is_latex_output()) {
  options(knitr.table.format = "latex")
} else {
  options(knitr.table.format = "pipe")
}
```


## Data Preprocessing and Descriptive Statistics

```{r}
#| label: run install scripts
#| echo: false
#| output: false
#| warning: false
#| error: false
#| results: false

# Run install.r to ensure packages are installed
source("install.r")
```

```{r}
#| label: load-Health-Avoiders
#| warning: false
#| include: true
#| echo: false

library(readr)

alldata <-
  read_csv("data/alldata.csv") %>%
  mutate(across(where(is.character), as.factor))
```

## Predictor Distributions

### Predictors

```{r selectdata}
#| label: selectdata-and-renaming-variable
#| echo: false

selectdata <- alldata %>%
  select("Ethnicity", "Political Party", "Gender", "Job Classification", "Education", 
         "Age", "Personal Income", "Race", "MacArthur Scale", "Social Media Usage", "AI Use", 
         'Video Games Hours', "Listening to Podcasts", "Facebook Usage", "Frequency of TikTok Use", 
         "X (Twitter) Usage", "Social Media", "Influencer Following",
         "Stressful Events - Reaction", "Stressful Events - Guilt", "Stressful Events - Detachment",
         "Stressful Events - Nightmares", "Stressful Events - Avoiding Situations",
         "Stressful Events - Recent Occurence", "Current Depression", 
         "Anxiety - Trouble Relaxing", "Anxiety - Irritable", "Anxiety - Restlessness",
         "Anxiety - Feeling Afraid", "Anxiety - Nervousness", "Anxiety - Worrying", 
         "Health - Feeling Failure", "Health - Hopelessness", "Health - Feeling Tired",
         "Health - Interest In Things", "Health - Pace", "Health - Poor Appetite",
         "Health - Thoughts Of Self Infliction", "Health - Concentration", "Health - Trouble Sleeping",
         "Anxiety - Worrying Different Things", "Stress Related Events - Most Stressful", "Stress Related Events",
         "Medical Diagnoses In Life", "Fast Food Consumption", "Physical Activity Guidelines", "Cigarette Smoking",
         "Supplement Consumption Reason", "Vaccinations", "Diet Type", "Supplement Consumption", "Meditation", 
         "Home Ownership", "Voter Registration", "Climate Change Belief", "Mental Health of Partner",
         "Information Avoidance - Cancer 1", "Information Avoidance - Cancer 2", 
         "Information Avoidance - Cancer 3 (R)", "Information Avoidance - Cancer 4", 
         "Information Avoidance - Cancer 5 (R)", "Information Avoidance - Cancer 6", 
         "Information Avoidance - Cancer 7 (R)", "Information Avoidance - Cancer 8 (R)") %>%
  rename(
    Political_Party = `Political Party`,
    Job_Classification = `Job Classification`,
    Personal_Income = `Personal Income`,
    MacArthur_Scale = `MacArthur Scale`,
    
    Social_Media_Usage = `Social Media Usage`,
    AI_Use = `AI Use`,
    Video_Games_Hours = `Video Games Hours`,
    Listening_Podcasts = `Listening to Podcasts`,
    Facebook_Usage = `Facebook Usage`,
    TikTok_Use = `Frequency of TikTok Use`,
    X_Twitter_Usage = `X (Twitter) Usage`,
    Social_Media_type =  `Social Media`,
    Influencer_Following = `Influencer Following`,
    
    Stressful_Events_Reaction = `Stressful Events - Reaction`,
    Stressful_Events_Guilt = `Stressful Events - Guilt`,
    Stressful_Events_Detachment = `Stressful Events - Detachment`,
    Stressful_Events_Nightmares = `Stressful Events - Nightmares`,
    Stressful_Events_Avoiding_Situations = `Stressful Events - Avoiding Situations`,
    Stressful_Events_Recent = `Stressful Events - Recent Occurence`,
    Current_Depression = `Current Depression`,
    Stressful_Events_Detachment = `Stressful Events - Detachment`,
    
    Anxiety_Trouble_Relaxing = `Anxiety - Trouble Relaxing`,
    Anxiety_Irritable = `Anxiety - Irritable`,
    Anxiety_Restlessness = `Anxiety - Restlessness`,
    Anxiety_Feeling_Afraid = `Anxiety - Feeling Afraid`,
    Anxiety_Nervousness = `Anxiety - Nervousness`,
    Anxiety_Worrying = `Anxiety - Worrying`,
    
    Health_Feeling_Failure = `Health - Feeling Failure`,
    Health_Hopelessness = `Health - Hopelessness`,
    Health_Feeling_Tired = `Health - Feeling Tired`,
    Health_Interest_In_Things = `Health - Interest In Things`,
    Health_Pace = `Health - Pace`,
    Health_Poor_Appetite = `Health - Poor Appetite`,
    Health_Thoughts_Of_Self_Infliction = `Health - Thoughts Of Self Infliction`,
    Health_Concentration = `Health - Concentration`,
    Health_Trouble_Sleeping = `Health - Trouble Sleeping`,
    Anxiety_Worrying_Different_Things = `Anxiety - Worrying Different Things`,
    
    Stressful_Events_Most = `Stress Related Events - Most Stressful`,
    Stress_Related_Events = `Stress Related Events`,
    Medical_Diagnoses_In_Life = `Medical Diagnoses In Life`,
    
    Fast_Food_Consumption = `Fast Food Consumption`,
    Physical_Activity_Guidelines = `Physical Activity Guidelines`,
    Cigarette_Smoking = `Cigarette Smoking`,
    Supplement_Consumption_Reason = `Supplement Consumption Reason`,
    Diet_Type = `Diet Type`,
    Supplement_Consumption = `Supplement Consumption`,
    
    Home_Ownership = `Home Ownership`,
    Voter_Registration = `Voter Registration`,
    Climate_Change_Belief = `Climate Change Belief`,
    Mental_Health_of_Partner = `Mental Health of Partner`
  )
```

```{r}
#| label: factoring-variables

selectdata <- selectdata %>%
  mutate(
    across(where(is.character), as.factor),
    Ethnicity = factor(Ethnicity, levels = c(
      "Yes, Cuban",
      "No, not of Hispanic, Latino, or Spanish origin",
      "Yes, Mexican, Mexican Am., Chicano",
      "Yes, another Hispanic, Latino, or Spanish origin – (for example, Salvadoran, Dominican, Colombian, Guatemalan, Spaniard, Ecuadorian, etc.)",
      "Yes, Puerto Rican",
      "Prefer not to say")),
    Political_Party = factor(Political_Party, levels = c(
      "Republican", "Democrat", "Independent", "Something else", "Prefer not to say")),
    Gender = factor(Gender, level = c(
        "Woman",
        "Man",
        "Non-binary",
        "Agender",
        "Two-spirit",
        "Additional gender category/identity not listed",
        "Prefer not to say")),
    Job_Classification = factor(Job_Classification, levels = c(
      "White Collar (e.g., accountant, software developer, human resources manager, marketing analyst, public safety)", "Creative or Cultural (e.g., writer, artist, musician, actor)", "Public Service and Government (e.g., soldier, teacher, police officer, public health worker, education, childcare)", "Information Technology (e.g., IT specialist, web developer, data scientist, cybersecurity analyst)", "Service Industry (e.g., retail worker, server, hotel staff, flight attendant, food services, personal care, funeral services, animal/veterinary care, leisure and hospitality)", "Freelance and Gig Economy (e.g., freelance writer, graphic designer, rideshare driver, delivery person)", "Blue Collar (e.g., electrician, plumber, mechanic, welder, manufacturing, oil and gas extraction, transportation, utilities, mining, waste collection/treatment/disposal, automotive services)", "Professional (e.g., doctor, lawyer, professor, engineer, nurse, healthcare)", "I am unemployed/a student/stay at home parent", " Manual Labor (e.g., farmer, construction worker, factory worker, miner, landscaping, agriculture, cleaning and custodial services)")),
    Education = factor(Education, levels = c(
      "No formal education",
      "Less than a high school diploma",
      "High school graduate - high school diploma or the equivalent (for example: GED)",
      "Some college, but no degree",
      "Associate degree (for example: AA, AS)",
      "Bachelor's degree (for example: BA, AB, BS)",
      "Master's degree (for example: MA, MS, MEng, MEd, MSW, MBA)",
      "Professional degree (for example: MD, DDS, DVM, LLB, JD)",
      "Doctorate degree (for example: PhD, EdD)",
      "Prefer not to say")),
     AgeGroup = case_when(
      Age < 25 ~ "18–24",
      Age < 35 ~ "25–34",
      Age < 45 ~ "35–44",
      Age < 55 ~ "45–54",
      Age < 65 ~ "55–64",
      TRUE     ~ "65+"
    ),
    AgeGroup = factor(AgeGroup, levels = c("18–24", "25–34", "35–44", "45–54", "55–64", "65+")),
    AgeBand = if_else(Age < 35, "Under 35", "35+"),
    AgeBand = factor(AgeBand, levels = c("Under 35", "35+")),
    Personal_Income = factor(Personal_Income, levels = c(
      "Less than $10,000",
        "$10,000-$19,999",
        "$20,000-$29,999",
        "$30,000-$39,999",
        "$40,000-$49,999",
        "$50,000-$59,999",
        "$60,000-$69,999",
        "$70,000-$79,999",
        "$80,000-$89,999",
        "$90,000-$99,999",
        "$100,000-$124,999",
        "$125,000-$149,999",
        "$150,000-$174,999",
        "$175,000-$199,999",
        "$200,000-$224,999",
        "$225,000-$249,999",
        "$250,000 or more",
        "Prefer not to say")),
    Race = factor(Race, levels = c(
    "American Indian or Alaska Native",
    "Asian Indian", "Chinese", "Filipino", "Japanese", "Korean", "Vietnamese", "Samoan", "Guamanian", "Hawaiian",
    "Black or African American",
    "White",
    "An ethnicity not listed here",
    "Other",
    "Prefer not to say")),
    MacArthur_Scale = factor(MacArthur_Scale, levels = paste("Rung", 1:10)),
    Social_Media_Usage = factor(Social_Media_Usage, levels = c(
             "Never",
             "A couple of times a year",
             "A couple of times a month",
             "Every week",
             "Every day")),
    Listening_Podcasts = factor(Listening_Podcasts, levels = c(
    "I have never listened to a podcast",
    "I do not regularly listen to podcasts",
    "At least every 4 weeks (22-30 days)",
    "At least every 3 weeks (15-21 days)",
    "At least every 2 weeks (8-14 days)",
    "At least every week (7 days)")),
    AI_Use = factor(AI_Use, levels = c(
      "Yes",
      "No")),
    Video_Games_Hours = factor(Video_Games_Hours, levels = c(
    "I do not play video games",
    "1-5 hours",
    "6-10 hours",
    "11-15 hours",
    "16-20 hours",
    "20+ hours")),
    Facebook_Usage = factor(Facebook_Usage, levels = c(
    "Never",
    "A few times a year",
    "A few times a month",
    "At least once a week",
    "Daily")),
    TikTok_Use = factor(TikTok_Use, levels = c(
    "Once a week or less",
    "Once a day",
    "A few times a week",
    "Several times a day")),
    X_Twitter_Usage = factor(X_Twitter_Usage, levels = c(
    "Never",
    "A few times a year",
    "A few times a month",
    "At least once a week",
    "Daily")),
    Influencer_Following = factor(Influencer_Following, levels = c(
    "No",
    "Unsure",
    "Yes")),
    Stressful_Events_Recent = factor(Stressful_Events_Recent, levels = c(
      "No, it was within the last 30 days", 
      "I did not experience a stressful or distressing event within the last 30 days", 
      "Yes, it occurred more than 30 days ago")),
    Current_Depression = factor(Current_Depression, levels = c(
      "Yes",
      "No")),
    Stressful_Events_Detachment = factor(Stressful_Events_Detachment, level = c(
      "Yes",
      "No")),
    Anxiety_Trouble_Relaxing = factor(Anxiety_Trouble_Relaxing, levels = c(
             "Not at all",
             "Several days",
             "More than half the days",
             "Nearly every day")),
    Anxiety_Feeling_Afraid = factor(Anxiety_Feeling_Afraid, levels = c(
             "Not at all",
             "Several days",
             "More than half the days",
             "Nearly every day")),
    Anxiety_Nervousness = factor(Anxiety_Nervousness, levels = c(
             "Not at all",
             "Several days",
             "More than half the days",
             "Nearly every day")),
    Anxiety_Worrying = factor(Anxiety_Worrying, levels = c(
             "Not at all",
             "Several days",
             "More than half the days",
             "Nearly every day")),
    Stressful_Events_Most = factor(Stressful_Events_Most, levels = c(
    "Life-threatening illness or injury",
    "Exposure to toxic substance (for example, dangerous chemicals, radiation)",
    "Sudden accidental death",
    "Transportation accident (for example, car accident, boat accident, train wreck, plane crash)",
    "Serious accident at work, home, or during recreational activity",
    "Captivity (for example, being kidnapped, abducted, held hostage, prisoner of war)",
    "Severe human suffering",
    "Combat or exposure to a war-zone (in the military or as a civilian)",
    "Assault with a weapon (for example, being shot, stabbed, threatened with a knife, gun, bomb)",
    "Sudden violent death (for example, homicide, suicide)",
    "Sexual assault (rape, attempted rape, made to perform any type of sexual act through force or threat of harm)",
    "Serious injury, harm, or death you caused to someone else",
    "Physical assault (for example, being attacked, hit, slapped, kicked, beaten up)",
    "Any other very stressful event or experience",
    "Other unwanted or uncomfortable sexual experience",
    "Fire or explosion",
    "Natural disaster (for example, flood, hurricane, tornado, earthquake)")),
    Fast_Food_Consumption = factor(Fast_Food_Consumption, levels = c(
    "I never eat fast food",
    "Less than once a month",
    "Once a month",
    "2-3 times a month",
    "Once a week or more")),
    Vaccinations = factor(Vaccinations, levels = c(
      "No",
      "It depends",
      "Yes")),
    Meditation = factor(Meditation, levels = c(
    "I've never done this before",
    "I've tried it in the past, but it wasn't for me",
    "I've done this in the past, but not regularly",
    "I did this regularly in the past, but not currently",
    "I do this at least once a month")),
    Physical_Activity_Guidelines = factor(Physical_Activity_Guidelines, levels = c(
        "No, and I do not intend to in the next 6 months.",
        "No, but I intend to in the next 6 months.",
        "No, but I intend to in the next 30 days.",
        "Yes, I have been for LESS than 6 months.",
        "Yes, I have been for MORE than 6 months.")),
    Cigarette_Smoking = factor(Cigarette_Smoking, levels = c(
        "I don't smoke cigarettes",
        "Less than one a day",
        "1 to 3 cigarettes",
        "4 to 6 cigarettes",
        "7 to 10 cigarettes",
        "More than 10 cigarettes")),
    Supplement_Consumption_Reason = factor(Supplement_Consumption_Reason, levels = c(
        "I do not take any supplements",
        "For maintaining good health in general",
        "For a specific issue")),
    Diet_Type = factor(Diet_Type, levels = c(
        "Vegetarian (I do not eat meat or fish)",
        "Vegan (I do not eat any animal products)",
        "Pollotarian (I do not eat red meat and fish, but eat poultry and fowl)",
        "None of the above",
        "Flexitarian (I eat vegetarian, but also occasionally meat or fish)",
        "Pescetarian (I do not eat meat, but I do eat fish)")),
    Supplement_Consumption = factor(Supplement_Consumption, levels = c(
      "I do not take any supplments", 
      "1", 
      "2", 
      "3", 
      "4 or more")),
    Home_Ownership = factor(Home_Ownership, levels = c(
    "Rent",
    "Lease",
    "Mortgage",
    "Full Ownership",
    "Other")),
    Voter_Registration = factor(Voter_Registration, levels = c(
        "Yes",
        "No, but I am eligible to vote.",
        "No, and I am not eligible to vote.")),
    Climate_Change_Belief = factor(Climate_Change_Belief, levels = c(
        "Strongly skeptical of claims about climate change and its link to human activities.",
        "Somewhat skeptical about the impact of human activities on climate change, believing that climate change is a natural cycle.",
        "Uncertain about the causes and extent of climate change.",
        "No opinion on the matter.",
        "Somewhat believe climate change is occurring and is influenced by human activities, but natural factors also play a significant role.",
        "Strongly believe climate change is occurring and is primarily caused by human activities.")),
    Mental_Health_of_Partner = factor(
      Mental_Health_of_Partner,
      levels = c(
        "My partner has not been diagnosed with a mental illness",
        "Less than a year",
        "1-2 years",
        "2-3 years",
        "3 or more years"
      )
    )
    
)
```

```{r}
#| label: simplify-job-class
#| echo: false


selectdata <- selectdata %>%
  mutate(
    Job_Classification = recode_factor(Job_Classification,
      `White Collar (e.g., accountant, software developer, human resources manager, marketing analyst, public safety)` = "White Collar",
      `Creative or Cultural (e.g., writer, artist, musician, actor)` = "Creative/Cultural",
      `Public Service and Government (e.g., soldier, teacher, police officer, public health worker, education, childcare)` = "Public Service/Government",
      `Information Technology (e.g., IT specialist, web developer, data scientist, cybersecurity analyst)` = "IT",
      `Service Industry (e.g., retail worker, server, hotel staff, flight attendant, food services, personal care, funeral services, animal/veterinary care, leisure and hospitality)` = "Service Industry",
      `Freelance and Gig Economy (e.g., freelance writer, graphic designer, rideshare driver, delivery person)` = "Freelance/Gig",
      `Blue Collar (e.g., electrician, plumber, mechanic, welder, manufacturing, oil and gas extraction, transportation, utilities, mining, waste collection/treatment/disposal, automotive services)` = "Blue Collar",
      `Professional (e.g., doctor, lawyer, professor, engineer, nurse, healthcare)` = "Professional",
      `I am unemployed/a student/stay at home parent` = "Unemployed/Student/Parent",
      ` Manual Labor (e.g., farmer, construction worker, factory worker, miner, landscaping, agriculture, cleaning and custodial services)` = "Manual Labor"
    )
  )
```

```{r}
#| label: dummy politial party
#| echo: false

library(dplyr)

selectdata <- selectdata %>%
  mutate(
    Republican = if_else(Political_Party == "Republican", 1, 0),
    Democrat = if_else(Political_Party == "Democrat", 1, 0),
    Independent = if_else(Political_Party == "Independent", 1, 0),
    Something_else = if_else(Political_Party == "Something else", 1, 0),
    Prefer_not_to_say = if_else(Political_Party == "Prefer not to say", 1, 0)
  )

```

```{r}
#| label: distribution-of-age
#| include: false
#| fig-height: 6
#| echo: false

#  (some of their age greater than 100)

library(ggplot2)

ggplot(selectdata, aes(x = Age)) +
  geom_histogram(binwidth = 5, color = "black") +
  labs(title = "Distribution of Age", x = "Age", y = "Count")

```

```{r}
library(ggplot2)

selectdata %>%
  count(AgeBand) %>%
  ggplot(aes(x = AgeBand, y = n)) +
  geom_col() +
  labs(title = "Age Band Distribution",
       x = "Age Band",
       y = "Count")

```

```{r}
#| label: distribution-of-personal-income 

# need renumercing the level

library(ggplot2)

ggplot(selectdata, aes(x = Personal_Income)) +
  geom_bar(fill = "steelblue", color = "black") +
  labs(title = "Distribution of Personal Income", 
       x = "Personal Income", 
       y = "Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

```

```{r}
#| label: numercing-gender-education-age-income-and-MacArthur-scale-variables

selectdata <- selectdata %>%
  mutate(
    Gender4 = case_when(
      Gender == "Woman" ~ "Woman",
      Gender == "Man" ~ "Man",
      Gender == "Non-binary" ~ "Non-binary",
      TRUE ~ "Other"   # all remaining categories go here
    ),
    Gender4 = factor(Gender4, levels = c("Woman", "Man", "Non-binary", "Other")),
    
  Education_Level = case_when(
      Education == "No formal education" ~ 0,
      Education == "Less than a high school diploma" ~ 1,
      Education == "High school graduate - high school diploma or the equivalent (for example: GED)" ~ 2,
      Education == "Some college, but no degree" ~ 3,
      Education == "Associate degree (for example: AA, AS)" ~ 4,
      Education == "Bachelor's degree (for example: BA, AB, BS)" ~ 5,
      Education == "Master's degree (for example: MA, MS, MEng, MEd, MSW, MBA)" ~ 6,
      Education == "Professional degree (for example: MD, DDS, DVM, LLB, JD)" ~ 7,
      Education == "Doctorate degree (for example: PhD, EdD)" ~ 7,      
      Education == "Perfer not to say" ~ 8,
      TRUE ~ NA
    ),
  Age = as.numeric(Age),
  Income = case_when(
    Personal_Income == "Less than $10,000" ~ 1,
    Personal_Income %in% c("$10,000-$19,999", "$20,000-$29,999") ~ 2,
    Personal_Income %in% c("$30,000-$39,999", "$40,000-$49,999") ~ 3,
    Personal_Income %in% c("$50,000-$59,999", "$60,000-$69,999") ~ 4,
    Personal_Income %in% c("$70,000-$79,999", "$80,000-$89,999", "$90,000-$99,999") ~ 5,
    Personal_Income %in% c("$100,000-$124,999", "$125,000-$149,999") ~ 6,
    Personal_Income %in% c("$150,000-$174,999", "$175,000-$199,999", "$200,000-$224,999", 
                           "$225,000-$249,999", "$250,000 or more") ~ 6,
    Personal_Income == "Prefer not to say" ~ 7,
    TRUE ~ NA
  ),
  
    MacArthur_Numeric = as.numeric(gsub("Rung ", "", MacArthur_Scale)),
  
    # Political Party: 1-4
    Political_Party_Group = case_when(
      Political_Party == "Republican" ~ 1,
      Political_Party == "Democrat" ~ 2,
      Political_Party == "Independent" ~ 3,
      Political_Party == "Something else" ~ 4,
      TRUE ~ NA_real_
    ),
    
    # Gender: 1-4
    Gender4 = case_when(
      Gender4 == "Woman" ~ 1,
      Gender4 == "Man" ~ 2,
      Gender4 == "Non-binary" ~ 3,
      Gender4 == "Other" ~ 4,
      TRUE ~ NA_real_
    ),
    
    # Race: white = 1, people of color = 0
    Race2 = case_when(
      Race %in% c("White") ~ 1,
      Race %in% c("American Indian or Alaska Native", "Asian Indian", "Chinese", "Filipino", 
                  "Japanese", "Korean", "Vietnamese", "Samoan", "Guamanian", "Hawaiian", 
                  "Black or African American", "An ethnicity not listed here", "Other") ~ 0,
      TRUE ~ NA_real_
    )
  ) %>%
  filter(Age > 18 & Age <= 100)
```

```{r}
#| label: numercing-summing-and-scaling-anxiety-score

selectdata <- selectdata %>%
  mutate(
    across(c(Anxiety_Trouble_Relaxing,
             Anxiety_Irritable,
             Anxiety_Restlessness,
             Anxiety_Feeling_Afraid,
             Anxiety_Nervousness,
             Anxiety_Worrying,
             Anxiety_Worrying_Different_Things),
           ~ case_when(
             . == "Not at all" ~ 0,
             . == "Several days" ~ 1,
             . == "More than half the days" ~ 2,
             . == "Nearly every day" ~ 3,
             TRUE ~ NA_real_
           ),
           .names = "{.col}_num"),
    
    # total score (0–21)
    Anxiety_Total_Score = Anxiety_Trouble_Relaxing_num +
                          Anxiety_Irritable_num +
                          Anxiety_Restlessness_num +
                          Anxiety_Feeling_Afraid_num +
                          Anxiety_Nervousness_num +
                          Anxiety_Worrying_num +
                          Anxiety_Worrying_Different_Things_num,
    
    # categorical version
    Anxiety_Severity_cat = case_when(
      Anxiety_Total_Score >= 0  & Anxiety_Total_Score <= 4  ~ "Minimal",
      Anxiety_Total_Score >= 5  & Anxiety_Total_Score <= 9  ~ "Mild",
      Anxiety_Total_Score >= 10 & Anxiety_Total_Score <= 14 ~ "Moderate",
      Anxiety_Total_Score >= 15 & Anxiety_Total_Score <= 21 ~ "Severe",
      TRUE ~ NA_character_
    ),
    
    # numerical version
     Anxiety_Severity_num = case_when(
      Anxiety_Total_Score >= 0  & Anxiety_Total_Score <= 4  ~ 1,
      Anxiety_Total_Score >= 5  & Anxiety_Total_Score <= 9  ~ 2,
      Anxiety_Total_Score >= 10 & Anxiety_Total_Score <= 14 ~ 3,
      Anxiety_Total_Score >= 15 & Anxiety_Total_Score <= 21 ~ 4,
      TRUE ~ NA_real_
    )
  )

```

```{r}
#| label: numercing-summing-and-scaling-stressful(ptsd)-score

selectdata <- selectdata %>%
  mutate(
    across(c(Stressful_Events_Reaction,
             Stressful_Events_Guilt,
             Stressful_Events_Detachment,
             Stressful_Events_Nightmares,
             Stressful_Events_Avoiding_Situations),
           ~ case_when(
             . == "Yes" ~ 1,
             . == "No"  ~ 0,
             TRUE ~ NA_real_
           ),
           .names = "{.col}_num"),
    
    # PC-PTSD-5 score: sum of the 5 recoded items
    PTSD5_Score = Stressful_Events_Reaction_num +
                     Stressful_Events_Guilt_num +
                     Stressful_Events_Detachment_num +
                     Stressful_Events_Nightmares_num +
                     Stressful_Events_Avoiding_Situations_num,
    
    # PTSD5 categorical
    PTSD5_cat = case_when(
      PTSD5_Score == 0 ~ "None",
      PTSD5_Score == 1 ~ "Minimal",
      PTSD5_Score == 2 ~ "Mild",
      PTSD5_Score == 3 ~ "Moderate",
      PTSD5_Score >= 4 ~ "Severe",
      TRUE ~ NA_character_
    )
  )

```

```{r}
#| label: numercing-summing-and-scaling-health-score

selectdata <- selectdata %>%
  mutate(
    across(c(Health_Feeling_Failure,
             Health_Hopelessness,
             Health_Feeling_Tired,
             Health_Interest_In_Things,
             Health_Pace,
             Health_Poor_Appetite,
             Health_Thoughts_Of_Self_Infliction,
             Health_Concentration,
             Health_Trouble_Sleeping),
           ~ case_when(
             . == "Not at all" ~ 0,
             . == "Several days" ~ 1,
             . == "More than half the days" ~ 2,
             . == "Nearly every day" ~ 3,
             TRUE ~ NA_real_
           ),
           .names = "{.col}_num"),
    
    # total score (range: 0–27 if 9 items)
    Health_Total_Score = Health_Feeling_Failure_num +
                         Health_Hopelessness_num +
                         Health_Feeling_Tired_num +
                         Health_Interest_In_Things_num +
                         Health_Pace_num +
                         Health_Poor_Appetite_num +
                         Health_Thoughts_Of_Self_Infliction_num +
                         Health_Concentration_num +
                         Health_Trouble_Sleeping_num,
    
    # categorical version
    Health_Depression_Severity = case_when(
      Health_Total_Score >= 0  & Health_Total_Score <= 4   ~ "None to Minimal Depression",
      Health_Total_Score >= 5  & Health_Total_Score <= 9   ~ "Mild Depression",
      Health_Total_Score >= 10 & Health_Total_Score <= 14  ~ "Moderate Depression",
      Health_Total_Score >= 15 & Health_Total_Score <= 19  ~ "Moderately Severe Depression",
      Health_Total_Score >= 20 & Health_Total_Score <= 27  ~ "Severe Depression",
      TRUE ~ NA_character_
    ),
    
    # numeric version of severity
    Health_Depression_Severity_num = case_when(
      Health_Total_Score >= 0  & Health_Total_Score <= 4   ~ 1,
      Health_Total_Score >= 5  & Health_Total_Score <= 9   ~ 2,
      Health_Total_Score >= 10 & Health_Total_Score <= 14  ~ 3,
      Health_Total_Score >= 15 & Health_Total_Score <= 19  ~ 4,
      Health_Total_Score >= 20 & Health_Total_Score <= 27  ~ 5,
      TRUE ~ NA_real_
    )
  )

```

```{r}
#| label: summing-life-events-anxiety score

library(dplyr)
library(stringr)
library(tidyr)

selectdata <- selectdata %>%
  mutate(
    StressEvents_Count = str_count(Stress_Related_Events, ";") + 1,
    MostStressful_Count = str_count(Stressful_Events_Most, ";") + 1,
    Stress_TotalScore = replace_na(StressEvents_Count, 0) + replace_na(MostStressful_Count, 0)
  )


```

```{r}
#| label: recode-cirgarette-smoking-meditation-group-and-facebook-usage-in-numeric

selectdata <- selectdata %>%
  mutate(
    Cigarette_Smoking_num = factor(case_when(
      Cigarette_Smoking %in% c("Less than one a day", "1 to 3 cigarettes", 
                               "4 to 6 cigarettes", "7 to 10 cigarettes", 
                               "More than 10 cigarettes") ~ 1,
      Cigarette_Smoking == "I don't smoke cigarettes" ~ 0,
      TRUE ~ NA_real_
    ), levels = c(0, 1)),
    
    Supplement_Consumption_Reason_num = factor(case_when(
      Supplement_Consumption_Reason %in% c("For a specific issue", 
                                           "For maintaining good health in general") ~ 1,
      Supplement_Consumption_Reason == "I do not take any supplements" ~ 0,
      TRUE ~ NA_real_
    ), levels = c(0, 1)),
    
    Meditation_group = factor(case_when(
      Meditation == "I've never done this before" ~ "No to past",
      Meditation %in% c("I've tried it in the past, but it wasn't for me",
                        "I've done this in the past, but not regularly",
                        "I did this regularly in the past, but not currently") ~ "Yes to past",
      Meditation == "I do this at least once a month" ~ "Yes to current",
      TRUE ~ NA_character_
    ), levels = c("No to past", "Yes to past", "Yes to current")),
    
    Facebook_Usage_cat = factor(case_when(
      Facebook_Usage == "Daily" ~ 1,
      Facebook_Usage %in% c("Never", "A few times a year", "A few times a month", 
                            "At least once a week") ~ 0,
      TRUE ~ NA_real_
    ), levels = c(0, 1))
  )

```

## Outcome Distribution and Item Correlations

```{r}
#| label: new-cancer-columns-with-new-names
selectdata <- selectdata %>%
  mutate(
    Avoid_Cancer1 = `Information Avoidance - Cancer 1`,
    Avoid_Cancer2 = `Information Avoidance - Cancer 2`,
    Avoid_Cancer3 = `Information Avoidance - Cancer 3 (R)`,
    Avoid_Cancer4 = `Information Avoidance - Cancer 4`,
    Avoid_Cancer5 = `Information Avoidance - Cancer 5 (R)`,
    Avoid_Cancer6 = `Information Avoidance - Cancer 6`,
    Avoid_Cancer7 = `Information Avoidance - Cancer 7 (R)`,
    Avoid_Cancer8 = `Information Avoidance - Cancer 8 (R)`
  )
## the R represent in the Name means reverse code
```

```{r}
#| label: convert-8-cancer-item-to-numeric
library(forcats)

# Reorder ALL 8 cancer items at once (1–4, higher = stronger agreement)
selectdata <- selectdata %>%
  mutate(across(starts_with("Avoid_Cancer"),
                ~ as.factor(.) %>%
                  fct_relevel("Strongly disagree",
                              "Somewhat disagree",
                              "Somewhat agree",
                              "Strongly agree") %>%
                  as.numeric()))

```

-   The outcome variable, `Cancer_Avoidance_Mean`, was computed as the average of 8 items measuring participants’ avoidance of cancer-related health information. Four of these items `Avoid_Cancer3`, `Avoid_Cancer5`, `Avoid_Cancer7`, `Avoid_Cancer8` were reverse-coded to ensure higher scores consistently reflected greater avoidance. Reverse scoring was calculated using the formula:

    -   Reversed score=(Max+Min)−Original score

```{r}
#| label: Reverse-the-R-coded-cancer-items
selectdata <- selectdata %>%
  mutate(
    Avoid_Cancer3 = 5 - Avoid_Cancer3,
    Avoid_Cancer5 = 5 - Avoid_Cancer5,
    Avoid_Cancer7 = 5 - Avoid_Cancer7,
    Avoid_Cancer8 = 5 - Avoid_Cancer8
  )
```

#### Correlation

```{r}
#| label: cancer-items-correlations
cancer_items <- select(selectdata, starts_with("Avoid_Cancer"))

# Correlation matrix
cor_matrix <- cor(cancer_items, use = "complete.obs")

print(cor_matrix)

```

```{r}
#| label: cancer-items-correlations-plot
library(corrplot)

corrplot(cor_matrix, 
         method = "color",
         type = "upper",
         order = "original",
         tl.cex = 0.8,
         tl.col = "black",
         tl.srt = 45,
         addCoef.col = "black",
         number.cex = 0.7)

# Alternative visualization with different style
corrplot(cor_matrix, 
         method = "circle",
         type = "full",
         order = "original",
         tl.cex = 0.8,
         tl.col = "black",
         tl.srt = 45,
         col = colorRampPalette(c("blue", "white", "red"))(100))
```

-   All the 8 of the cancer items are moderate to high correlate with each other, so we can add them up and get the average score.


#### Alpha

```{r}
#| label: cancer-items-alpha
library(psych)

alpha_identity <- psych::alpha(selectdata[, c("Avoid_Cancer1", "Avoid_Cancer2", "Avoid_Cancer3", "Avoid_Cancer4","Avoid_Cancer5", "Avoid_Cancer6", "Avoid_Cancer7", "Avoid_Cancer8")])
print(alpha_identity)
```

-   Based on the analysis, the alpha does not improve for the measure if any of the items are dropped. All final alphas will be equal or lower than the 0.89 raw alpha and 0.9 standardized alpha.

-   The internal consistency of the `Information Avoidance – Cancer` items was high (std.alpha = 0.895), suggesting that the items measured a coherent construct. The average inter-item correlation was (average_r = 0.55), indicating moderately strong associations among the items without redundancy. Reliability estimates (G6 = 0.894) further supported this consistency.

```{r}
#| label: average-the-8-cancer-items

# If one scale:
selectdata <- selectdata %>%
  mutate(Cancer_Avoidance_Mean = rowMeans(cancer_items, na.rm = TRUE))
```

```{r}
#| label: distribution-of-the-average-cancer-avoidance-score

library(ggplot2)

ggplot(selectdata, aes(x = Cancer_Avoidance_Mean)) +
  geom_histogram(aes(fill = after_stat(x)), binwidth = 0.5, color = "black") +
  scale_x_continuous(
    breaks = seq(1, 4, 0.5),
    limits = c(1, 4)
  ) +
  theme_minimal(base_size = 13) +
  scale_fill_gradientn(
    colours = c("white", "mistyrose", "lightcoral", "red", "darkred"),
    values = scales::rescale(c(1, 2, 3, 3.5, 4)),
    limits = c(1, 4),
    name = "Avoidance Level"
  ) +
  labs(
    x = "Average Cancer Avoidance Score",
    y = "Count",
    title = "Distribution of the Average Cancer Avoidance Score"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = "bottom",
    legend.key.width = unit(2, "cm")
  )
```

Percentage of participants with Cancer_Avoidance_Mean \>= 2.5

```{r}
#| label: percentage-of-participants-with-higher-cancer-avoidance-score
# Percentage of participants with Cancer_Avoidance_Mean >= 2.5
mean_25_or_higher <- mean(selectdata$Cancer_Avoidance_Mean >= 2.5, na.rm = TRUE) * 100
mean_25_or_higher

```

```{r}
#| label: binary-and-transformation-cancer-avoidance-mean-score

selectdata <- selectdata %>%
  mutate(
    # Binary variable: 1 if mean >= 3 (higher avoidance), 0 if 1–2
    Cancer_Avoiders01 = ifelse(Cancer_Avoidance_Mean >= 3, 1, 0),
    Cancer_Avoiders01 = factor(Cancer_Avoiders01, levels = c(0, 1)),
    
    # Log transformation (since variable is positive)
    Cancer_Avoidance_log = log10(Cancer_Avoidance_Mean + 1),
    
    # Square root transformation
    Cancer_Avoidance_sqrt = sqrt(Cancer_Avoidance_Mean)
  )
```

```{r}
#| label: plots-to-compare-original-log-sqrt

plot_df <- selectdata %>%
  select(Cancer_Avoidance_Mean, Cancer_Avoidance_log, Cancer_Avoidance_sqrt) %>%
  pivot_longer(everything(), names_to = "version", values_to = "value")

# histograms
ggplot(plot_df, aes(x = value)) +
  geom_histogram(bins = 30) +
  facet_wrap(~version, scales = "free") +
  labs(title = "Original vs log vs sqrt", x = NULL, y = "Count")

# densities
ggplot(plot_df, aes(x = value)) +
  geom_density() +
  facet_wrap(~version, scales = "free") +
  labs(title = "Density: original vs log vs sqrt", x = NULL, y = "Density")

```

#### Matrix

```{r}
#| label: correlation-cancer-mean-main-predictors
#| echo: false
#| message: false
#| warning: false

library(dplyr)
library(corrplot)
library(knitr)

exclude_patterns <- c("_num$", 
                      paste0("Avoid_Cancer", 1:8),
                      "Cancer_Avoiders01",
                      "Cancer_Avoidance_log",
                      "Cancer_Avoidance_sqrt",
                      "MostStressful_Count")

# Select only numeric vars
numeric_vars_clean <- selectdata %>%
  select(where(is.numeric)) %>%
  select(-matches(paste(exclude_patterns, collapse="|")))

# Recalculate correlation matrix
cor_matrix <- cor(numeric_vars_clean, use = "complete.obs")

# Drop self-correlation
avoidance_corrs <- cor_matrix[
  "Cancer_Avoidance_Mean",
  setdiff(colnames(cor_matrix), "Cancer_Avoidance_Mean"),
  drop = TRUE
]

# Convert to table

avoidance_corrs_df <- data.frame(
  Predictor = names(avoidance_corrs),
  Correlation = round(as.numeric(avoidance_corrs), 2)
)

kable(
  avoidance_corrs_df,
  caption = "Correlations between Cancer Avoidance Mean and main predictors",
  escape = TRUE,
  booktabs = TRUE
)
```

```{r}
#| label: plot-correlation-cancer-mean-with-main-predictors
#| fig-width: 8
#| fig-height: 9

library(ggplot2)

ggplot(avoidance_corrs_df, aes(x = reorder(Predictor, Correlation), y = Correlation)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Correlations of Cancer Avoidance Mean with Main Predictors",
       x = "Main Predictors",
       y = "Correlation") +
  theme_minimal(base_size = 12)
```

```{r}
#| label: save-selectdata
write.csv(selectdata, "data/selectdata.csv", row.names = FALSE)
```







## Descriptive tables for categorical variables

```{r}
#| label: categorical-table

library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
library(stringr)

# select all categorical variables
categorical_vars <- selectdata %>%
  select(where(~ is.factor(.) || is.character(.)))  

# counts and percentages, missing values
categorical_table <- categorical_vars %>%
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value") %>%
  mutate(Value = ifelse(is.na(Value), "Missing", Value)) %>%  
  group_by(Variable, Value) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Variable) %>%
  mutate(percent = paste0(round(100 * n / sum(n), 1), "%")) %>%
  arrange(Variable, Value) %>%
  ungroup()

categorical_table_head <- head(categorical_table, 15)

kable(
  categorical_table_head,
  caption = "Sample of Categorical Variable Counts and Percentages",
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```

## REGRESSION MODELS

(continuous outcome `Cancer_Avoidance_Mean`)

### Random Forest Analysis (tidymodel)

#### Demographic Model

```{r}
#| label: random-forest-demographic-model



library(tidymodels)
library(tidyr)
library(vip)

# need to drop NA to get accuracy
demo_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Ethnicity, Political_Party, Gender4, Job_Classification,
    Education_Level, AgeBand, Income, Race, MacArthur_Numeric
  )


# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(demo_data)), size = 0.7 * nrow(demo_data))
train <- demo_data[train_idx, ]
test  <- demo_data[-train_idx, ]

# Fit random forest model (only for demo)
rf_demo <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("regression")

# Set recipe
rf_demo_recipe <- recipe(
  Cancer_Avoidance_Mean ~  Ethnicity + Political_Party + Gender4 + Job_Classification + 
    Education_Level + AgeBand + Income + Race + MacArthur_Numeric, 
  data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 

# Create workflow
rf_demo_workflow <- workflow() %>%
  add_model(rf_demo) %>%
  add_recipe(rf_demo_recipe)

# Fit model
rf_demo_fit <- rf_demo_workflow %>%
  fit(data = train)


# Predict on test set
pred <- predict(rf_demo_fit, test) %>%
  bind_cols(test %>% select(Cancer_Avoidance_Mean))


# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoidance_Mean, estimate = .pred)


# Root Mean Squared Error
rmse <- pred %>%
  rmse(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)


# Mean Absolute Error
mae <- pred %>%
  mae(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)


# Correlation between predicted and actual
cor <- cor(pred$.pred, pred$Cancer_Avoidance_Mean)


# Variable importance
rf_demo_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)
```
```{r}
#| label: random-forest-demographic-model-kable

# Create performance table
rsq_val <- pred %>%
  rsq(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)

performance_df <- data.frame(
  Metric = c("RMSE", "MAE", "R-squared", "Correlation"),
  Value = c(rmse, mae, rsq_val, cor)
)

kable(
  performance_df,
  caption = "Random Forest Model Performance (Demographic Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```


-   RMSE (\~0.63) and MAE (\~0.53) are pretty close, meaning errors aren’t heavily dominated by outliers.

-   Correlation = 0.109 is extremely low. It basically means that predicted values do not track the actual values at all. So the model is not capturing the pattern in the test data.

-   If predictors have very weak relationships with the outcome, the model will predicts values near the mean of the training set.

-   The random forest model shows that `MacArthur_Numeric` is the most predict variable.



```{r}
#| label: linear-regression-model-MacArthur_Numeric-predicting-Cancer-Avoidance-Mean
#| warning: false


MacArthur_cancer_linear <- lm(
  Cancer_Avoidance_Mean ~ MacArthur_Numeric,
  data = demo_data
)

library(knitr)
library(kableExtra)

fit <- lm(Cancer_Avoidance_Mean ~ MacArthur_Numeric, data = demo_data)

# Use broom to get clean output
tidy_results <- tidy(fit)

# Clean names
tidy_results$term <- c("Intercept", "MacArthur Numeric")

kable(
  tidy_results,
  col.names = c("Term", "Estimate", "Std Error", "Statistic", "p value"),
  digits = c(3,3,3,3,10),
  booktabs = TRUE,
  escape = TRUE
) %>%
  kable_styling(latex_options = "hold_position")

```

-   The linear regression analysis shows a statistically significant negative relationship between MacArthur Scale scores and cancer avoidance mean ($\beta$ = -0.025, p = 1.89e-10). For each one-unit increase in the MacArthur score (indicating higher subjective socioeconomic status), cancer avoidance mean scores decreased by 0.025 units.

```{r}
library(ggplot2)

ggplot(demo_data, aes(x = MacArthur_Numeric, y = Cancer_Avoidance_Mean)) +
  geom_point(alpha = 0.6, color = "black") +
  geom_smooth(method = "lm", se = TRUE, color = "red", size = 1.2) +
  geom_jitter(height = 0, width = 0.2, alpha = 0.2, color = "black") +
  scale_x_continuous(breaks = 1:10, limits = c(1, 10)) +
  labs(
    x = "MacArthur Numeric Score",
    y = "Cancer Avoidance Mean",
    title = "Linear Regression: Cancer Avoidance vs MacArthur Score"
  ) +
  theme_minimal()

```

-   The scatter plot displays a slight negative trend (indicated by the red regression line), but the data points are widely dispersed across all MacArthur score levels.

-   This visual pattern confirms the weak correlation, showing that while a statistical relationship exists, MacArthur score (subjective socioeconomic status) is not a strong predictor of cancer avoidance behaviors on its own.

```{r}
#| label: linear-regression-model-AgeGroup-predicting-Cancer-Avoidance-Mean
age_cancer_linear <- lm(Cancer_Avoidance_Mean ~ AgeGroup, data = demo_data)

coef_df <- as.data.frame(summary(age_cancer_linear)$coefficients)

knitr::kable(
  coef_df,
  caption = "Linear Regression: Predicting Cancer Avoidance Mean from Age Group",
  digits = 3
)
```

-   The initial analysis examined cancer avoidance scores across six age categories, revealing that middle-aged adults (35-64) showed significantly higher cancer avoidance behaviors compared to the youngest group (18-24). However, the model explains only 0.33% of the variance ($R^2$= 0.003), indicating weak predictive power despite statistical significance.

-   So we break the `AgeGroup` into `AgeBand` with only 2 categories "above 35" and "under 35" to see if the result will be better.

```{r}
#| label: linear-regression-model-AgeBand-predicting-Cancer-Avoidance-Mean
age_cancer_linear <- lm(Cancer_Avoidance_Mean ~ AgeBand, data = demo_data)

coef_df <- as.data.frame(summary(age_cancer_linear)$coefficients)

knitr::kable(
  coef_df,
  caption = "Linear Regression: Predicting Cancer Avoidance Mean from Age Band",
  digits = 5
)
```

-   Individuals aged 35 and older had cancer avoidance scores 0.063 units higher than those under 35 (p = 6.45e-06), a highly significant difference. The model remains weak in explanatory power ($R^2$ = 0.0025), but the clear age threshold at 35 years provides a meaningful distinction for further investigation.

-   Now, we wanted to see if other predictors would influence the `AgeBand` for under 35 and over 35, so we created two multivariable linear regression models (Age_under35 and Age_over35) predicting Cancer Avoidance Mean.

```{r}
#| label: multivariable-linear-regression-model-Age_under35-predicting-Cancer-Avoidance-Mean

model_under35 <- lm(Cancer_Avoidance_Mean ~ Education_Level + Income + 
                    MacArthur_Numeric + Political_Party,
                    data = dplyr::filter(demo_data, AgeBand == "Under 35"))

coef_df <- as.data.frame(summary(model_under35)$coefficients)

knitr::kable(
  coef_df,
  caption = "Linear Regression: Under 35 Predicting Cancer Avoidance Mean",
  digits = 3
)
```

-   For individuals under 35, political party affiliation emerged as the strongest predictor of cancer avoidance. Democrats showed 0.165 lower scores than Republicans (p = 8.97e-08), Independents showed 0.120 lower scores (p = 0.0003), and those selecting "Something else" showed 0.155 lower scores (p = 0.0001).

-   But the overall model explaining less than 1% of variance ($R^2$ = 0.0097).

```{r}
library(ggeffects)

pred_party <- ggpredict(model_under35, terms = "Political_Party")
ggplot(pred_party, aes(x = x, y = predicted, fill = x)) +
  geom_col() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  labs(
    x = "Political Party",
    y = "Predicted Cancer Avoidance Mean",
    title = "Effect of Political Party on Cancer Avoidance (Age < 35)"
  ) +
  theme_minimal()

```

-   The bar plots shows predicted cancer avoidance means across political affiliations, stratified by age. For those under 35, Republicans show the highest cancer avoidance behaviors (\~1.80), while all other political groups cluster around 1.65, with error bars indicating moderate variability.

```{r}
#| label: multivariable-linear-regression-model-Age-over35-predicting-Cancer-Avoidance-Mean

model_over35 <- lm(Cancer_Avoidance_Mean ~ Education_Level + Income + 
                   MacArthur_Numeric + Political_Party,
                   data = dplyr::filter(demo_data, AgeBand == "35+"))

coef_df <- as.data.frame(summary(model_over35)$coefficients)

knitr::kable(
  coef_df,
  caption = "Linear Regression: Over 35 Predicting Cancer Avoidance Mean",
  digits = 3
)
```

-   For individuals aged 35 and older, political party affiliation showed even stronger effects than in the younger cohort. Democrats had 0.225 lower cancer avoidance scores than Republicans (p = 2.67e-15), Independents had 0.131 lower scores (p = 1.28e-05), and "Something else" respondents had 0.138 lower scores (p = 0.0007).

-   Additionally, both income ($\beta$ = -0.016, p = 0.022) and MacArthur score ($\beta$ = -0.029, p = 7.00e-06) showed significant negative associations, with this model explaining 2.8% of the variance ($R^2$ = 0.028), the highest among all models tested.

```{r}
pred_party <- ggpredict(model_over35, terms = "Political_Party")
ggplot(pred_party, aes(x = x, y = predicted, fill = x)) +
  geom_col() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  labs(
    x = "Political Party",
    y = "Predicted Cancer Avoidance Mean",
    title = "Effect of Political Party on Cancer Avoidance (Age > 35)"
  ) +
  theme_minimal()

```

-   The bar plots shows predicted cancer avoidance means across political affiliations, stratified by age. For those 35 and older, Republicans maintain the highest scores (\~1.93), while Democrats show notably lower scores (\~1.68), and other groups fall in between.

-   There has been a consistent gap between Republicans and Democrats in the two age groups. Republicans exhibited higher levels of cancer avoidance scores, suggesting that political party is strongly associated with these avoidance behaviors regardless of age, although this effect is more pronounced in older adults.

#### Media Use Model

```{r}
#| label: random-forest-meida-use-model
#| fig-width: 8
#| fig-height: 4


library(tidymodels)
library(tidyr)
library(vip)


# need to drop NA to get accuracy
media_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Social_Media_Usage, AI_Use, Video_Games_Hours, Listening_Podcasts,
    Facebook_Usage_cat, TikTok_Use, X_Twitter_Usage, Social_Media_type, Influencer_Following
  )

# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(media_data)), size = 0.7 * nrow(media_data))
train <- media_data[train_idx, ]
test  <- media_data[-train_idx, ]

# Fit random forest model (only for demo)
rf_media <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("regression")

# Set recipe
rf_media_recipe <- recipe(
  Cancer_Avoidance_Mean ~  Social_Media_Usage + AI_Use + Video_Games_Hours + Listening_Podcasts +
    Facebook_Usage_cat + TikTok_Use + X_Twitter_Usage + Influencer_Following,
  data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 

# Create workflow
rf_media_workflow <- workflow() %>%
  add_model(rf_media) %>%
  add_recipe(rf_media_recipe)

# Fit model
rf_media_fit <- rf_media_workflow %>%
  fit(data = train)

# Predict on test set
pred <- predict(rf_media_fit, test) %>%
  bind_cols(test %>% select(Cancer_Avoidance_Mean))


# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoidance_Mean, estimate = .pred)


# Root Mean Squared Error
rmse <- pred %>%
  rmse(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)


# Mean Absolute Error
mae <- pred %>%
  mae(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)


# Correlation between predicted and actual
cor <- cor(pred$.pred, pred$Cancer_Avoidance_Mean)


# Variable importance
rf_media_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)

```

```{r}
# Create performance table
rsq_val <- pred %>%
  rsq(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)

performance_df <- data.frame(
  Metric = c("RMSE", "MAE", "R-squared", "Correlation"),
  Value = c(rmse, mae, rsq_val, cor)
)

kable(
  performance_df,
  caption = "Random Forest Model Performance (Media Usage Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```

-   RMSE (\~0.66) and MAE (\~0.55) are pretty close, meaning errors aren’t heavily dominated by outliers.

-   Correlation = 0.021 is extremely low. It basically means that predicted values do not track the actual values at all. So the model is not capturing the pattern in the test data.

-   If predictors have very weak relationships with the outcome, the model will predicts values near the mean of the training set.

-   The random forest model shows that `Influencer_Following` is the most predictive variable.

```{r}
#| label: linear-regression-model-Influencer_Following-predicting-Cancer_Avoidance_Mean
influencer_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Influencer_Following, data = media_data)

coef_df <- as.data.frame(summary(influencer_cancer_linear)$coefficients)

knitr::kable(
  coef_df,
  caption = "Linear Regression: Influencer Following Predicting Cancer Avoidance Mean",
  digits = 3
)
```

-   The linear regression analysis shows a none meaningful relationship between Influencer_Following scores and cancer avoidance mean (p \> 0.05).

```{r}
#| label: scatter-plot-Influencer_Following-predicting-Cancer_Avoidance_Mean

library(ggplot2)

ggplot(media_data, aes(x = factor(Influencer_Following), y = Cancer_Avoidance_Mean)) +
  geom_boxplot(fill = "white", alpha = 0.6) +
  geom_jitter(aes(color = Cancer_Avoidance_Mean), width = 0.1, alpha = 0.6) +
  scale_color_gradientn(
    colours = c("white", "mistyrose", "lightcoral", "red", "darkred"),
    values = scales::rescale(c(1, 2, 3, 3.5, 4)),
    limits = c(1, 4),
    name = "Avoidance Level"
  ) +
  scale_x_discrete(labels = c("0" = "No", "1" = "Yes")) +
  labs(
    x = "Influencer Following",
    y = "Cancer Avoidance Mean",
    title = "Cancer Avoidance by Influencer Following"
  ) +
  theme_minimal()

```

-   The boxplot examining cancer avoidance scores across influencer following categories (No, Unsure, Yes) reveals remarkably similar distributions across all three groups. This visual pattern strongly supports the conclusion that influencer following status has no substantial relationship with cancer avoidance score.

#### Health Condition Model

```{r}
#| label: random-forest-Health-Condition-model
#| fig-width: 12
#| fig-height: 5


library(tidymodels)
library(tidyr)
library(vip)


# need to drop NA to get accuracy
health_condition_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Stressful_Events_Recent, Current_Depression, Anxiety_Severity_num, PTSD5_Score,
    Health_Depression_Severity_num, Stress_TotalScore
  )


# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(health_condition_data)), size = 0.7 * nrow(health_condition_data))
train <- health_condition_data[train_idx, ]
test  <- health_condition_data[-train_idx, ]


# Fit random forest model (only for health condition)
rf_health_condition <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("regression")


# Set recipe
rf_health_condition_recipe <- recipe(
  Cancer_Avoidance_Mean ~ Stressful_Events_Recent + Current_Depression + Anxiety_Severity_num +
  PTSD5_Score +  Health_Depression_Severity_num + Stress_TotalScore,
  data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 


# Create workflow
rf_health_condition_workflow <- workflow() %>%
  add_model(rf_health_condition) %>%
  add_recipe(rf_health_condition_recipe)

# Fit model
rf_health_condition_fit <- rf_health_condition_workflow %>%
  fit(data = train)

# Predict on test set
pred <- predict(rf_health_condition_fit, test) %>%
  bind_cols(test %>% select(Cancer_Avoidance_Mean))


# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoidance_Mean, estimate = .pred)


# Root Mean Squared Error
rmse <- pred %>%
  rmse(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)


# Mean Absolute Error
mae <- pred %>%
  mae(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)


# Correlation between predicted and actual
cor <- cor(pred$.pred, pred$Cancer_Avoidance_Mean)


# Variable importance
rf_health_condition_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)

```

```{r}
# Create performance table
rsq_val <- pred %>%
  rsq(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)

performance_df <- data.frame(
  Metric = c("RMSE", "MAE", "R-squared", "Correlation"),
  Value = c(rmse, mae, rsq_val, cor)
)

kable(
  performance_df,
  caption = "Random Forest Model Performance (Health Condition Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```

-   RMSE (\~0.65) and MAE (\~0.55) are pretty close, meaning errors aren’t heavily dominated by outliers.

-   Correlation = 0.103 is extremely low. It basically means that predicted values do not track the actual values at all. So the model is not capturing the pattern in the test data.

-   If predictors have very weak relationships with the outcome, the model will predicts values near the mean of the training set.

-   The random forest model shows that `Stress_TotalScore` is the most predict variable.

```{r}
#| label: linear-regression-model-Stress_TotalScore-predicting-Cancer_Avoidance_Mean
stress_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Stress_TotalScore, data = health_condition_data)

coef_df <- as.data.frame(summary(stress_cancer_linear)$coefficients)

knitr::kable(
  coef_df,
  caption = "Linear Regression: Stress Score Predicting Cancer Avoidance Mean",
  digits = 3
)
```

-   The linear regression analysis shows a none meaningful relationship between stress total scores and cancer avoidance mean (p \> 0.05).

```{r}
#| label: scatter-plot-Stress_TotalScore-predicting-Cancer_Avoidance_Mean

library(ggeffects)

pred <- ggpredict(stress_cancer_linear, terms = "Stress_TotalScore")

ggplot(pred, aes(x = x, y = predicted)) +
  geom_line(color = "red", size = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  labs(
    x = "Stress Total Score",
    y = "Predicted Cancer Avoidance Mean",
    title = "Predicted Cancer Avoidance by Stress Score"
  ) +
  theme_minimal()

```

-   The plot displays predicted cancer avoidance scores across the range of stress total scores (0-18), showing a slight positive trend indicated by the red regression line. The broad gray confidence band around the regression line reflects the considerable uncertainty in these predictions. While higher stress levels may be weakly associated with slightly higher cancer avoidance behaviors, but this relationship is not statistically significant.

#### Health Behavior Model

```{r}
#| label: random-forest-Health-Behavior-model
#| fig-width: 12
#| fig-height: 4


library(tidymodels)
library(tidyr)
library(vip)


# need to drop NA to get accuracy
health_behavior_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Fast_Food_Consumption, Meditation_group, Physical_Activity_Guidelines,
    Cigarette_Smoking_num, Supplement_Consumption_Reason_num, Diet_Type, Supplement_Consumption
  )


# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(health_behavior_data)), size = 0.7 * nrow(health_behavior_data))
train <- health_behavior_data[train_idx, ]
test  <- health_behavior_data[-train_idx, ]

# Fit random forest model (only for health_condition)
rf_health_behavior <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("regression")


# Set recipe
rf_health_behavior_recipe <- recipe(
  Cancer_Avoidance_Mean ~ Fast_Food_Consumption + Meditation_group + Physical_Activity_Guidelines +
    Cigarette_Smoking_num + Supplement_Consumption_Reason_num + Diet_Type + Supplement_Consumption,
data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 


# Create workflow
rf_health_behavior_workflow <- workflow() %>%
  add_model(rf_health_behavior) %>%
  add_recipe(rf_health_behavior_recipe)

# Fit model
rf_health_behavior_fit <- rf_health_behavior_workflow %>%
  fit(data = train)

# Predict on test set
pred <- predict(rf_health_behavior_fit, test) %>%
  bind_cols(test %>% select(Cancer_Avoidance_Mean))


# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoidance_Mean, estimate = .pred)


# Root Mean Squared Error
rmse <- pred %>%
  rmse(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)


# Mean Absolute Error
mae <- pred %>%
  mae(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)


# Correlation between predicted and actual
cor <- cor(pred$.pred, pred$Cancer_Avoidance_Mean)


# Variable importance
rf_health_behavior_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)

```

```{r}
# Create performance table
rsq_val <- pred %>%
  rsq(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)

performance_df <- data.frame(
  Metric = c("RMSE", "MAE", "R-squared", "Correlation"),
  Value = c(rmse, mae, rsq_val, cor)
)

kable(
  performance_df,
  caption = "Random Forest Model Performance (Health Behavior Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```


-   RMSE (\~0.67) and MAE (\~0.56) are pretty close, meaning errors aren’t heavily dominated by outliers.

-   Correlation = 0.06 is extremely low. It basically means that predicted values do not track the actual values at all. So the model is not capturing the pattern in the test data.

-   If predictors have very weak relationships with the outcome, the model will predicts values near the mean of the training set.

-   The random forest model shows that number of `Cigarette_Smoking` is the most predict variable.

```{r}
#| label: linear-regression-model-Cigarette_Smoking_num-predicting-Cancer_Avoidance_Mean
smoking_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Cigarette_Smoking_num, data = health_behavior_data)

coef_df <- as.data.frame(summary(smoking_cancer_linear)$coefficients)

knitr::kable(
  coef_df,
  caption = "Linear Regression: Smoking Predicting Cancer Avoidance Mean",
  digits = 3
)
```

-   The linear regression analysis shows a statistically significant positive relationship between cigarette smoking and cancer avoidance behaviors ($\beta$ = 0.186, p \< 2e-16). For every one-unit increase in cigarette smoking status (from non-smoker to smoker), the predicted cancer avoidance score increases by 0.186 units, meaning individuals who smoke tend to have slightly higher cancer avoidance scores.

-   While this relationship is conventionally significant at the 0.05 level, the model explains only 1.1% of the variance ($R^2$ = 0.011), indicating that smoking status alone provides minimal predictive power for cancer avoidance scores.

```{r}
#| label: scatter-plot-Cigarette_Smoking_num-predicting-Cancer_Avoidance_Mean

library(ggplot2)

ggplot(health_behavior_data, aes(x = factor(Cigarette_Smoking_num), y = Cancer_Avoidance_Mean)) +
  geom_boxplot(fill = "white", alpha = 0.6) +
  geom_jitter(aes(color = Cancer_Avoidance_Mean), width = 0.1, alpha = 0.6) +
  scale_color_gradientn(
    colours = c("white", "mistyrose", "lightcoral", "red", "darkred"),
    values = scales::rescale(c(1, 2, 3, 3.5, 4)),
    limits = c(1, 4),
    name = "Avoidance Level"
  ) +
  labs(
    x = "Cigarette Smoking (0 = No, 1 = Yes)",
    y = "Cancer Avoidance Mean",
    title = "Cancer Avoidance by Cigarette Smoking Status"
  ) +
  theme_minimal()
```

-   Look at the box plot, smokers show higher cancer avoidance scores on average compared to non-smokers. This could mean that smokers are more aware of or concerned about cancer risks, possibly leading them to engage in other cancer-preventive behaviors. However, the low variance reminds us that while this pattern is statistically significant, smoking status explains only a small fraction of the variation in cancer avoidance scores.

#### Other Model

```{r}
#| label: random-forest-other-model
#| fig-width: 15
#| fig-height: 4


library(tidymodels)
library(tidyr)
library(vip)

# need to drop NA to get accuracy
other_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Home_Ownership, Voter_Registration, Climate_Change_Belief,
    Mental_Health_of_Partner, Education_Level, Political_Party
  )


# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(other_data)), size = 0.7 * nrow(other_data))
train <- other_data[train_idx, ]
test  <- other_data[-train_idx, ]

# Fit random forest model (only for other variable)
rf_other <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("regression")


# Set recipe
rf_other_recipe <- recipe(
  Cancer_Avoidance_Mean ~  Climate_Change_Belief + Education_Level + Political_Party,
  data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 


# Create workflow
rf_other_workflow <- workflow() %>%
  add_model(rf_other) %>%
  add_recipe(rf_other_recipe)

# Fit model
rf_other_fit <- rf_other_workflow %>%
  fit(data = train)

# Predict on test set
pred <- predict(rf_other_fit, test) %>%
  bind_cols(test %>% select(Cancer_Avoidance_Mean))


# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoidance_Mean, estimate = .pred)

# Root Mean Squared Error
rmse <- pred %>%
  rmse(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)

# Mean Absolute Error
mae <- pred %>%
  mae(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)

# Correlation between predicted and actual
cor <- cor(pred$.pred, pred$Cancer_Avoidance_Mean)

# Variable importance
rf_other_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)

```

```{r}
# Create performance table
rsq_val <- pred %>%
  rsq(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)

performance_df <- data.frame(
  Metric = c("RMSE", "MAE", "R-squared", "Correlation"),
  Value = c(rmse, mae, rsq_val, cor)
)

kable(
  performance_df,
  caption = "Random Forest Model Performance (Other Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```


-   RMSE (\~0.63) and MAE (\~0.52) are pretty close, meaning errors aren’t heavily dominated by outliers.

-   Correlation = 0.145 is extremely low. It basically means that predicted values do not track the actual values at all. So the model is not capturing the pattern in the test data.

-   If predictors have very weak relationships with the outcome, the model will predicts values near the mean of the training set.

-   The random forest model shows that `Climate_Change_Belief` (Strongly believe climate change is occurring and is primarily caused by human activities) is the most predict variable.

```{r}
#| label: linear-regression-model-Climate_Change_Belief-predicting-Cancer_Avoidance_Mean
#| wrap: true
#| tidy: true

voter_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Climate_Change_Belief, data = other_data)

coef_df <- as.data.frame(summary(voter_cancer_linear)$coefficients)

knitr::kable(
  coef_df,
  caption = "Linear Regression: Cliamte Change Belief Predicting Cancer Avoidance Mean",
  digits = 3
)
```

-   The linear regression shows significant associations between climate change beliefs and cancer avoidance scores. Individuals who strongly believe in human-caused climate change showed significantly lower cancer avoidance scores ($\beta$ = -0.271, p = 3.32e-13), meaning their scores were 0.271 units lower on average. The model explains 2.4% of the variance ($R^2$ = 0.024), indicating a weak but statistically significant relationship.

```{r}
#| label: scatter-plot-Climate_Change_Belief-predicting-Cancer_Avoidance_Mean
#| fig-width: 9

library(ggplot2)

ggplot(data = other_data, aes(x = Climate_Change_Belief, y = Cancer_Avoidance_Mean)) +
  geom_boxplot(fill = "white", alpha = 0.6) +
  geom_jitter(aes(color = Cancer_Avoidance_Mean), width = 0.1, alpha = 0.6) +
  scale_color_gradientn(
    colours = c("white", "mistyrose", "lightcoral", "red", "darkred"),
    values = scales::rescale(c(1, 2, 3, 3.5, 4)),
    limits = c(1, 4),
    name = "Avoidance Level"
  ) +
  labs(x = "Climate_Change_Belief", 
       y = "Cancer Avoidance Mean",
       title = "Relationship between Climate_Change_Belief and Cancer Avoidance") +
  scale_x_discrete(labels = c("Strongly skeptical of claims about climate change and its link to human activities." = "Strongly disbelieve", 
                              "Somewhat skeptical about the impact of human activities on climate change, believing that climate change is a natural cycle." = "Somewhat disbelieve", 
                              "Uncertain about the causes and extent of climate change." = "Uncertain",
                              "No opinion on the matter." = "No opinion",
                              "Somewhat believe climate change is occurring and is influenced by human activities, but natural factors also play a significant role." = "Somewhat believe",
                              "Strongly believe climate change is occurring and is primarily caused by human activities." = "Strongly believe")) +
  theme_minimal()
```

-   The boxplot displays cancer avoidance behaviors across six climate change belief categories, revealing a clear gradient pattern. Climate change deniers and no opinion show higher median cancer avoidance scores around 1.9 - 2, with relatively compact distributions. Which means Individuals who strongly believe in human-caused climate change showed significantly lower cancer avoidance scores. However, the low variance reminds us that while this pattern is statistically significant, undefined explains only a small fraction of the variation in cancer avoidance scores.

```{r}
#| include: FALSE
#| label: check-how-many-categorical-predictors-can-take 

sapply(selectdata, function(x) if(is.factor(x)) length(levels(x)) else NA)

## Can not handle categorical predictors with more than 53 categories. ASk Dr.Shane which variables can be remove, Medical_Diagnoses_In_Life need to be delete (that's too much!)
```

### Multivariate Adaptive Regression Splines ([MARS](https://bradleyboehmke.github.io/HOML/mars.html))

#### Demographic

```{r}
#| label: package-for-MARS
# Modeling packages
library(earth)     # for fitting MARS models
library(caret)     # for automating the tuning process

# Model interpretability packages
library(vip)       # for variable importance
library(pdp)       # for variable relationships
```

```{r}
#| label: mars-demographic-model
#| results: hide

# need to drop NA to get accuracy
demo_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Ethnicity, Political_Party, Gender4, Job_Classification,
    Education_Level, Age, Income, Race, MacArthur_Numeric
  )

# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(demo_data)), size = 0.7 * nrow(demo_data))
train <- demo_data[train_idx, ]
test  <- demo_data[-train_idx, ]

# Fit MARS model
mars_model <- earth(
  Cancer_Avoidance_Mean ~ Ethnicity + Political_Party + Gender4 + Job_Classification +
    Education_Level + Age + Income + Race + MacArthur_Numeric,
  data = train
)

summary(mars_model)
# Predict on test set
pred <- predict(mars_model, newdata = test)

# Variable importance
evimp(mars_model)
```

```{r}
#| label: mars-demographic-model-kable

library(knitr)
library(kableExtra)

# need to drop NA to get accuracy
demo_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Ethnicity, Political_Party, Gender4, Job_Classification,
    Education_Level, Age, Income, Race, MacArthur_Numeric
  )

# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(demo_data)), size = 0.7 * nrow(demo_data))
train <- demo_data[train_idx, ]
test  <- demo_data[-train_idx, ]

# Fit MARS model
mars_model <- earth(
  Cancer_Avoidance_Mean ~ Ethnicity + Political_Party + Gender4 + Job_Classification +
    Education_Level + Age + Income + Race + MacArthur_Numeric,
  data = train
)

# Extract and display coefficients
coef_df <- data.frame(
  Term = names(coef(mars_model)),
  Coefficient = as.numeric(coef(mars_model))
)

kable(
  head(coef_df, 15),
  caption = "MARS Model Coefficients (Demographic Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Variable importance table
vip_result <- evimp(mars_model)

# Create clean data frame
vip_df <- data.frame(
  Variable = rownames(vip_result),
  nsubsets = vip_result[, 1],
  gcv = vip_result[, 2],
  rss = vip_result[, 3]
)

kable(
  vip_df,
  caption = "MARS Variable Importance (Demographic Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE,
  col.names = c("Variable", "N Subsets", "GCV", "RSS")
) %>%
  kable_styling(latex_options = "hold_position")

# Model summary statistics
summary_stats <- data.frame(
  Metric = c("Selected Terms", "R-squared", "GRSq"),
  Value = c(
    length(coef(mars_model)),
    summary(mars_model)$rsq,
    summary(mars_model)$grsq
  )
)

kable(
  summary_stats,
  caption = "MARS Model Summary Statistics",
  digits = 3,
  booktabs = TRUE
) %>%
  kable_styling(latex_options = "hold_position")
```

-   `Political_PartyDemocrat` is the most influential predictor (baseline = 100).

-   `MacArthur_Numeric` is about 75% as important.

-   `RaceWhite` is about 63%.

```{r}
#| label: mars-cross-validation-demographic
#| wrap: true
#| tidy: true
#| results: hide

library(earth)

mars_model <- earth(
  Cancer_Avoidance_Mean ~ Ethnicity + Political_Party + Gender4 + Job_Classification +
    Education_Level + Age + Income + Race + MacArthur_Numeric,
  data = train,
  degree = 2,       # allow up to 2-way interactions
  nfold = 10,       # 10-fold CV
  keepxy = TRUE
)
summary(mars_model) %>% .$coefficients %>% head(10)
```

```{r}
#| label: mars-cross-validation-demographic-kable

library(earth)
library(knitr)
library(kableExtra)

mars_model <- earth(
  Cancer_Avoidance_Mean ~ Ethnicity + Political_Party + Gender4 + Job_Classification +
    Education_Level + Age + Income + Race + MacArthur_Numeric,
  data = train,
  degree = 2,       # allow up to 2-way interactions
  nfold = 10,       # 10-fold CV
  keepxy = TRUE
)

# Extract coefficients
coef_matrix <- summary(mars_model)$coefficients
coef_df <- data.frame(
  Term = rownames(coef_matrix),
  Coefficient = coef_matrix[, 1]
)

kable(
  head(coef_df, 10),
  caption = "MARS Model Coefficients with Cross-Validation (Top 10 Terms)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Model summary statistics
summary_stats <- data.frame(
  Metric = c("Selected Terms", "R-squared", "GRSq"),
  Value = c(
    length(coef(mars_model)),
    summary(mars_model)$rsq,
    summary(mars_model)$grsq
  )
)

kable(
  summary_stats,
  caption = "MARS Model Summary Statistics",
  digits = 4,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Cross-validation table (if available)
if (!is.null(mars_model$cv.rsq.tab)) {
  cv_table <- as.data.frame(mars_model$cv.rsq.tab)
  cv_table$Terms <- rownames(cv_table)
  cv_table <- cv_table[, c("Terms", setdiff(names(cv_table), "Terms"))]
  
  kable(
    tail(cv_table, 5),  # Show last 5 rows
    caption = "Cross-Validation Results by Number of Terms",
    digits = 4,
    booktabs = TRUE,
    row.names = FALSE
  ) %>%
    kable_styling(latex_options = "hold_position")
}
```



Interactions appear as products (such as):

-   Ethnicity Prefer not to say \* Political_Party Democrat = an interaction term between ethnicity and political party.

-   Selected model at 17 terms

-   GRSq/RSq converging around 0.03-0.04

-   Mean out-of-fold RSq near zero (poor predictive performance)

-   High cross-validation error (CVRSq = 0.017, sd = 0.019)

-   This suggests demographic variables don’t explain Cancer_Avoidance_Mean very well.

```{r}
#| label: mars-importance-plot-demographic

library(ggplot2)

vip <- evimp(mars_model)

# Convert to data.frame for ggplot
vip_df <- data.frame(
  Variable = rownames(vip),
  GCV = vip[, "gcv"]
)

# Wrap labels at ~20 characters
vip_df$Variable_wrapped <- str_wrap(vip_df$Variable, width = 60)

ggplot(vip_df, aes(x = reorder(Variable_wrapped, GCV), y = GCV)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "MARS Variable Importance (GCV)",
       x = "Predictor",
       y = "Importance (scaled GCV)")
```

[Visualizing MARS models](https://rpubs.com/erblast/mars)

```{r}
#| label: mars-plot-demographic
plot(mars_model, which = 1)   # Predicted vs observed
```

-   The model selection plot shows that both RSq (training fit) and GRSq (generalized R-square) values remain low. This indicates that demographic characteristics provide limited explanatory power for cancer avoidance mean, and adding nonlinear terms does not meaningfully improve predictive accuracy.

-   Demographics status explain very little of the variation in cancer avoidance. Where effects appear, they are narrow interactions, such as between political affiliation and ethnicity, or between Job_Classification Blue Collar and MacArthur_Numeric \> 3.

#### Health condition

```{r}
#| label: mars-health-condition-model
#| results: hide

# need to drop NA to get accuracy
health_condition_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Stressful_Events_Recent, Current_Depression, Anxiety_Severity_num, PTSD5_Score,
    Health_Depression_Severity_num, Stress_TotalScore
  )

# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(health_condition_data)), size = 0.7 * nrow(health_condition_data))
train <- health_condition_data[train_idx, ]
test  <- health_condition_data[-train_idx, ]

# Fit MARS model
mars_model_health_condition <- earth(
    Cancer_Avoidance_Mean ~ Stressful_Events_Recent + Current_Depression + Anxiety_Severity_num +
    PTSD5_Score +  Health_Depression_Severity_num + Stress_TotalScore,
  data = train
)

summary(mars_model)
# Predict on test set
pred <- predict(mars_model_health_condition, newdata = test)

# Variable importance
evimp(mars_model_health_condition)
```

```{r}
#| label: mars-health-condition-model-kable

library(knitr)
library(kableExtra)

# need to drop NA to get accuracy
health_condition_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Stressful_Events_Recent, Current_Depression, Anxiety_Severity_num, PTSD5_Score,
    Health_Depression_Severity_num, Stress_TotalScore
  )

# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(health_condition_data)), size = 0.7 * nrow(health_condition_data))
train <- health_condition_data[train_idx, ]
test  <- health_condition_data[-train_idx, ]

# Fit MARS model
mars_model_health_condition <- earth(
    Cancer_Avoidance_Mean ~ Stressful_Events_Recent + Current_Depression + Anxiety_Severity_num +
    PTSD5_Score +  Health_Depression_Severity_num + Stress_TotalScore,
  data = train
)

# Extract coefficients
coef_df <- data.frame(
  Term = names(coef(mars_model_health_condition)),
  Coefficient = as.numeric(coef(mars_model_health_condition))
)

kable(
  coef_df,
  caption = "MARS Model Coefficients (Health Condition Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Variable importance table
vip_result <- evimp(mars_model_health_condition)

vip_df <- data.frame(
  Variable = rownames(vip_result),
  nsubsets = vip_result[, 1],
  gcv = vip_result[, 2],
  rss = vip_result[, 3]
)

kable(
  vip_df,
  caption = "MARS Variable Importance (Health Condition Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE,
  col.names = c("Variable", "N Subsets", "GCV", "RSS")
) %>%
  kable_styling(latex_options = "hold_position")

# Model summary statistics
summary_stats <- data.frame(
  Metric = c("Selected Terms", "R-squared", "GRSq"),
  Value = c(
    length(coef(mars_model_health_condition)),
    summary(mars_model_health_condition)$rsq,
    summary(mars_model_health_condition)$grsq
  )
)

kable(
  summary_stats,
  caption = "MARS Model Summary Statistics (Health Condition)",
  digits = 3,
  booktabs = TRUE
) %>%
  kable_styling(latex_options = "hold_position")

# Predict on test set
pred <- predict(mars_model_health_condition, newdata = test)
```


-   `Anxiety_Severity_num` is the most influential predictor (baseline = 100).

-   `PTSD5_Score` is about 61%

```{r}
#| label: mars-cross-validation-health-condition

library(earth)

mars_model_health_condition <- earth(
    Cancer_Avoidance_Mean ~ Stressful_Events_Recent + Current_Depression + Anxiety_Severity_num +
    PTSD5_Score +  Health_Depression_Severity_num + Stress_TotalScore,
  data = train,
  degree = 2,       # allow up to 2-way interactions
  nk = 100,
  nfold = 10,       # 10-fold CV
  keepxy = TRUE
)
pred <- predict(mars_model_health_condition, newdata = test)
rmse <- sqrt(mean((pred - test$Cancer_Avoidance_Mean)^2))
cor(pred, test$Cancer_Avoidance_Mean)


```


Interactions appear as products (such as):

-   h(Anxiety_Severity_num-3) = 0.191: When anxiety severity score \> 3, cancer avoidance increases by 0.191 units.

-   h(1-PTSD5_Score) = -0.088: When PTSD score \< 1, cancer avoidance decreases by 0.088 units, meaning higher PTSD scores are associated with slightly higher avoidance.

-   Selected model at 5 terms

-   GRSq/RSq barely above zero (GRSq = 0.0073, RSq = 0.012)

-   Mean out-of-fold RSq is negative - confirms no predictive power

-   CVRSq = -0.003, MaxErr = 2.35, This suggests health condition variables don’t explain Cancer_Avoidance_Mean very well.

-   Predictive power is weak: training R-square \~ 0.065, The model explains about 6.5% of the variance in Cancer_Avoidance_Mean on the training data and is a weak fit.

-   The correlation between Cancer_Avoidance_Mean and health condition is also very week (0.1452411)

```{r}
#| label: mars-importance-plot-health-condition

library(ggplot2)

vip <- evimp(mars_model_health_condition)

# Convert to data.frame for ggplot
vip_df <- data.frame(
  Variable = rownames(vip),
  GCV = vip[, "gcv"]
)

# Wrap labels at ~20 characters
vip_df$Variable_wrapped <- str_wrap(vip_df$Variable, width = 30)

ggplot(vip_df, aes(x = reorder(Variable_wrapped, GCV), y = GCV)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "MARS Variable Importance (GCV)",
       x = "Predictor",
       y = "Importance (scaled GCV)")

```

```{r}
#| label: mars-plot-health-condition

plot(mars_model_health_condition, which = 1)   # Predicted vs observed
```

-   The model selection plot shows that both RSq (training fit) and GRSq (generalized R-square) initially increase, indicating that adding the first few hinge functions improves model performance. However, after approximately 4 terms, GRSq reaches its maximum and then slightly declines, suggesting that additional terms provide minimal benefit and may cause overfitting.

-   The selected model at 4 terms uses 4 predictors and achieves a modest GRSq (0.01), indicating that although weak nonlinear relationships exist, the overall predictive power of the model remains low.

## CLASSIFICATION MODEL

(binary outcome `Cancer_Avoiders01`)

### Random Forest

#### Demographic Model

```{r}
#| label: random-forest-demographic-model_class

library(tidymodels)
library(tidyr)
library(vip)

# need to drop NA to get accuracy
demo_data <- selectdata %>%
  drop_na(
    Cancer_Avoiders01, Ethnicity, Political_Party, Gender4, Job_Classification,
    Education_Level, AgeBand, Income, Race, MacArthur_Numeric
  )


# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(demo_data)), size = 0.7 * nrow(demo_data))
train <- demo_data[train_idx, ]
test  <- demo_data[-train_idx, ]

# Fit random forest model (only for demo)
rf_demo <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("classification")

# Set recipe
rf_demo_recipe <- recipe(
  Cancer_Avoiders01 ~  Ethnicity + Political_Party + Gender4 + Job_Classification + 
    Education_Level + AgeBand + Income + Race + MacArthur_Numeric, 
  data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 

# Create workflow
rf_demo_workflow <- workflow() %>%
  add_model(rf_demo) %>%
  add_recipe(rf_demo_recipe)

# Fit model
rf_demo_fit <- rf_demo_workflow %>%
  fit(data = train)


# Predict on test set
pred_probs <- predict(rf_demo_fit, test, type = "prob")
pred_class <- predict(rf_demo_fit, test, type = "class")

# Combine everything
pred <- bind_cols(
  test %>% select(Cancer_Avoiders01),
  pred_class,
  pred_probs
)

# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoiders01, estimate = .pred_class)

# Confusion matrix
conf_mat_result <- conf_mat(pred, truth = Cancer_Avoiders01, estimate = .pred_class)

# ROC AUC
roc_auc_val <- roc_auc(pred, truth = Cancer_Avoiders01, .pred_1)

# Variable importance
rf_demo_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)
```
```{r}
# Create performance metrics table
performance_df <- data.frame(
  Metric = metrics$.metric,
  Value = metrics$.estimate
)

kable(
  performance_df,
  caption = "Random Forest Classification Model Performance (Demographic Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Create confusion matrix table
conf_mat_tbl <- as.data.frame(conf_mat_result$table)
colnames(conf_mat_tbl) <- c("Truth", "Prediction", "Count")

conf_mat_wide <- tidyr::pivot_wider(
  conf_mat_tbl,
  names_from = Prediction,
  values_from = Count
)

kable(
  conf_mat_wide,
  caption = "Confusion Matrix",
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Add ROC AUC to a separate table
roc_df <- data.frame(
  Metric = "ROC AUC",
  Value = roc_auc_val$.estimate
)

kable(
  roc_df,
  caption = "ROC AUC Score",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```

-   High accuracy (0.955) simply reflects predicting majority class.

-   ROC AUC = 0.449 , which means no better than random

-   The random forest model shows that `MacArthur_Numeric` is the most predict variable.

```{r}
#| label: logistic-regression-model-MacArthur_Numeric-predicting-Cancer-Avoiders01

MacArthur_cancer_logistic <- glm(Cancer_Avoiders01 ~ MacArthur_Numeric, data = demo_data, family = binomial)

coef_df <- as.data.frame(summary(MacArthur_cancer_logistic)$coefficients)

knitr::kable(
  coef_df,
  caption = "Logistic Regression: MacArthur Numeric Predicting Cancer Avoiders",
  digits = 3
)
```

-   The logistic regression shows a statistically significant negative relationship between MacArthur scores and the probability of being a cancer avoider ($\beta$ = -0.108, p = 0.0006). For each one-unit increase in the MacArthur score, the log-odds of being a cancer avoider decrease by 0.108, meaning individuals with higher subjective socioeconomic status are less likely to be classified as cancer avoiders.

-   However, the model shows minimal improvement with only 11.8 point reduction in deviance over the null model (null deviance = 2785.6, residual deviance = 2773.8).

```{r}
#| label: auc-macarthur
#| warning: false

library(pROC)
MacArthur_pred <- predict(MacArthur_cancer_logistic, type = "response")

roc_obj_MacArthur <- roc(demo_data$Cancer_Avoiders01, MacArthur_pred)

auc(roc_obj_MacArthur)
```

-   ROC AUC = 0.55, which means slightly better than random

```{r}
#| label: scatter-plot-MacArthur_Numeric-predicting-Cancer-Avoiders01

library(ggplot2)
library(ggeffects)

pred <- ggpredict(MacArthur_cancer_logistic, terms = "MacArthur_Numeric")

ggplot(pred, aes(x = x, y = predicted)) +
  geom_line(color = "red", size = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  scale_x_continuous(breaks = 1:10) +
  labs(
    x = "MacArthur Score",
    y = "Predicted Probability",
    title = "Predicted Probability of Cancer Avoiders by MacArthur Score"
  ) +
  theme_minimal()

```

-   The plot displays the predicted probability of being a cancer avoider across MacArthur scores ranging from 1 to 10. The red line shows a clear downward trend. The gray confidence band widens slightly at the extremes but remains relatively narrow, indicating consistent uncertainty across the score range. This visualization confirms the negative relationship identified in the logistic regression, higher socioeconomic status is associated with lower probability of being a cancer avoider.

-   However, the poor ROC AUC value (0.449) indicates that although the relationship is statistically significant, the model has little practical use in identifying cancer avoiders.

#### Media Usage Model

```{r}
#| label: random-forest-media-use-model-class

library(tidymodels)
library(tidyr)
library(vip)

# need to drop NA to get accuracy
media_data <- selectdata %>%
  drop_na(
    Cancer_Avoiders01, Social_Media_Usage, AI_Use, Video_Games_Hours, Listening_Podcasts,
    Facebook_Usage_cat, TikTok_Use, X_Twitter_Usage, Social_Media_type, Influencer_Following
  )


# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(media_data)), size = 0.7 * nrow(media_data))
train <- media_data[train_idx, ]
test  <- media_data[-train_idx, ]

# Fit random forest model (only for demo)
rf_media <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("classification")

# Set recipe
rf_media_recipe <- recipe(Cancer_Avoiders01 ~  Social_Media_Usage + AI_Use + Video_Games_Hours + Listening_Podcasts + Facebook_Usage_cat + TikTok_Use + X_Twitter_Usage + Influencer_Following, 
  data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 

# Create workflow
rf_media_workflow <- workflow() %>%
  add_model(rf_media) %>%
  add_recipe(rf_media_recipe)

# Fit model
rf_media_fit <- rf_media_workflow %>%
  fit(data = train)


# Predict on test set
pred_probs <- predict(rf_media_fit, test, type = "prob")
pred_class <- predict(rf_media_fit, test, type = "class")

# Combine everything
pred <- bind_cols(
  test %>% select(Cancer_Avoiders01),
  pred_class,
  pred_probs
)

# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoiders01, estimate = .pred_class)

# Confusion matrix
conf_mat_result <- conf_mat(pred, truth = Cancer_Avoiders01, estimate = .pred_class)

# ROC AUC
roc_auc_val  <- roc_auc(pred, truth = Cancer_Avoiders01, .pred_1)

# Variable importance
rf_media_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)
```
```{r}
# Create performance metrics table
performance_df <- data.frame(
  Metric = metrics$.metric,
  Value = metrics$.estimate
)

kable(
  performance_df,
  caption = "Random Forest Classification Model Performance (Media Usage Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Create confusion matrix table
conf_mat_tbl <- as.data.frame(conf_mat_result$table)
colnames(conf_mat_tbl) <- c("Truth", "Prediction", "Count")

conf_mat_wide <- tidyr::pivot_wider(
  conf_mat_tbl,
  names_from = Prediction,
  values_from = Count
)

kable(
  conf_mat_wide,
  caption = "Confusion Matrix",
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Add ROC AUC to a separate table
roc_df <- data.frame(
  Metric = "ROC AUC",
  Value = roc_auc_val$.estimate
)

kable(
  roc_df,
  caption = "ROC AUC Score",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```

-   High accuracy (0.947) simply reflects predicting majority class.

-   ROC AUC = 0.574 , which means slightly better than random

-   The random forest model shows that `Facebook_Usage` is the most predict variable.

```{r}
#| label: logistic-regression-model-Facebook_Usage_cat-predicting-Cancer_Avoiders01

facebook_cancer_logistic <- glm(Cancer_Avoiders01 ~ Facebook_Usage_cat, data = media_data, family = binomial)

coef_df <- as.data.frame(summary(facebook_cancer_logistic)$coefficients)

knitr::kable(
  coef_df,
  caption = "Logistic Regression: Facebook Usage Predicting Cancer Avoiders",
  digits = 3
)
```

-   The logistic regression shows a statistically significant positive relationship between Facebook usage and the probability of being a cancer avoider ($\beta$ = 0.2358, p = 0.0388). Facebook users have odds of being a cancer avoider that are 1.27 times higher ($e^0.2358$ = 1.266) than non-users, representing a 26.6% increase in odds.

-   However, the model shows minimal improvement over the null model, with only a 4.3 point reduction in deviance (null deviance = 2654.9, residual deviance = 2650.6), indicating that Facebook usage explains none of the variation in cancer avoider.

```{r}
#| label: auc-facebook
#| warning: false
facebook_pred <- predict(facebook_cancer_logistic, type = "response")

roc_obj_facebook <- roc(media_data$Cancer_Avoiders01, facebook_pred)

auc(roc_obj_facebook)
```

-   ROC AUC = 0.53, which means slightly better than random

```{r}
ggplot(media_data, aes(x = Facebook_Usage_cat, 
                       fill = factor(Cancer_Avoiders01))) +
  geom_bar(position = "fill") +  
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(
    values = c("0" = "lightgray", "1" = "red"),
    labels = c("No", "Yes")
  ) +
  labs(
    x = "Facebook Usage",
    y = "Percentage",
    fill = "Cancer Avoider",
    title = "Percentage of Cancer Avoiders by Facebook Usage"
  ) +
  theme_minimal()

```

-   The stacked bar chart displays the percentage distribution of cancer avoiders and non-avoiders across Facebook usage groups (0 = No, 1 = Yes). Both groups show similar distributions, with cancer avoiders (red) representing a very small percentage at the bottom of each bar. While the logistic regression identified a statistically significant difference (p = 0.039), the visual similarity between the two bars reinforces that this difference has minimal value for predicting cancer avoider status.

#### Health Condition Model

```{r}
#| label: random-forest-health-condition-model-class
#| fig-width: 12
#| fig-height: 5

library(tidymodels)
library(tidyr)
library(vip)

# need to drop NA to get accuracy
health_condition_data <- selectdata %>%
  drop_na(
    Cancer_Avoiders01, Stressful_Events_Recent, Current_Depression, Anxiety_Severity_num, PTSD5_Score,
    Health_Depression_Severity_num, Stress_TotalScore
  )


# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(health_condition_data)), size = 0.7 * nrow(health_condition_data))
train <- health_condition_data[train_idx, ]
test  <- health_condition_data[-train_idx, ]

# Fit random forest model (only for demo)
rf_health_condition <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("classification")

# Set recipe
rf_health_condition_recipe <- recipe(Cancer_Avoiders01 ~  Stressful_Events_Recent + Current_Depression + Anxiety_Severity_num + PTSD5_Score +  Health_Depression_Severity_num + Stress_TotalScore,
  data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 

# Create workflow
rf_health_condition_workflow <- workflow() %>%
  add_model(rf_health_condition) %>%
  add_recipe(rf_health_condition_recipe)

# Fit model
rf_health_condition_fit <- rf_health_condition_workflow %>%
  fit(data = train)


# Predict on test set
pred_probs <- predict(rf_health_condition_fit, test, type = "prob")
pred_class <- predict(rf_health_condition_fit, test, type = "class")

# Combine everything
pred <- bind_cols(
  test %>% select(Cancer_Avoiders01),
  pred_class,
  pred_probs
)

# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoiders01, estimate = .pred_class)


# Confusion matrix
conf_mat_result <- conf_mat(pred, truth = Cancer_Avoiders01, estimate = .pred_class)

# ROC AUC
roc_auc_val <- roc_auc(pred, truth = Cancer_Avoiders01, .pred_1)

# Variable importance
rf_health_condition_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)
```

```{r}
# Create performance metrics table
performance_df <- data.frame(
  Metric = metrics$.metric,
  Value = metrics$.estimate
)

kable(
  performance_df,
  caption = "Random Forest Classification Model Performance (Health COndition Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Create confusion matrix table
conf_mat_tbl <- as.data.frame(conf_mat_result$table)
colnames(conf_mat_tbl) <- c("Truth", "Prediction", "Count")

conf_mat_wide <- tidyr::pivot_wider(
  conf_mat_tbl,
  names_from = Prediction,
  values_from = Count
)

kable(
  conf_mat_wide,
  caption = "Confusion Matrix",
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Add ROC AUC to a separate table
roc_df <- data.frame(
  Metric = "ROC AUC",
  Value = roc_auc_val$.estimate
)

kable(
  roc_df,
  caption = "ROC AUC Score",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```

-   High accuracy (0.948) simply reflects predicting majority class.

-   ROC AUC = 0.446, which means no better than random

-   The random forest model shows that `Stress_TotalScore` is the most predict variable.

```{r}
#| label: logistic-regression-model-Stress_TotalScore-predicting-Cancer_Avoiders01

stress_cancer_logistic <- glm(Cancer_Avoiders01 ~ Stress_TotalScore, data = health_condition_data, family = binomial)

coef_df <- as.data.frame(summary(stress_cancer_logistic)$coefficients)

knitr::kable(
  coef_df,
  caption = "Logistic Regression: Stress Total Score Predicting Cancer Avoiders",
  digits = 3
)

```

-   The logistic regression shows a none significant relationship between stress total score and the probability of being a cancer avoider ($\beta$ = 0.01475, p = 0.447).

```{r}
#| label: auc-stress
#| warning: false

stress_pred <- predict(stress_cancer_logistic, type = "response")

roc_obj_stress <- roc(health_condition_data$Cancer_Avoiders01, stress_pred)

auc(roc_obj_stress)
```

-   ROC AUC = 0.51, which means slightly better than random

```{r}
pred <- ggpredict(stress_cancer_logistic, terms = "Stress_TotalScore")

# Plot
ggplot(pred, aes(x = x, y = predicted)) +
  geom_line(color = "red", size = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  scale_x_continuous(breaks = 1:15) +   # set x-axis from 0 to 15
  labs(
    x = "Stress Total Score",
    y = "Predicted Probability",
    title = "Predicted Probability of Cancer Avoiders by Stress Total Score"
  ) +
  theme_minimal()

```

-   The plot displays the predicted probability of being a cancer avoider across stress total scores ranging from 1 to 18. The red line shows a very slight upward trend. The gray confidence band widens at higher stress scores, indicating increasing uncertainty in predictions for individuals with very high stress levels. Which demonstrates that stress total score provides essentially no predictive information for identifying cancer avoiders, consistent with the lack of statistical significance in the model.

#### Health Behavior Model

```{r}
#| label: random-forest-health-behavior-model-class
#| fig-width: 12
#| fig-height: 5

library(tidymodels)
library(tidyr)
library(vip)

# need to drop NA to get accuracy
health_behavior_data <- selectdata %>%
  drop_na(
    Cancer_Avoiders01, Fast_Food_Consumption, Meditation_group, Physical_Activity_Guidelines,
    Cigarette_Smoking_num, Supplement_Consumption_Reason_num, Diet_Type, Supplement_Consumption
  )


# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(health_behavior_data)), size = 0.7 * nrow(health_behavior_data))
train <- health_behavior_data[train_idx, ]
test  <- health_behavior_data[-train_idx, ]

# Fit random forest model (only for demo)
rf_health_behavior <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("classification")

# Set recipe
rf_health_behavior_recipe <- recipe(Cancer_Avoiders01 ~  Fast_Food_Consumption + Meditation_group + Physical_Activity_Guidelines + Cigarette_Smoking_num + Supplement_Consumption_Reason_num + Diet_Type + Supplement_Consumption,
  data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 

# Create workflow
rf_health_behavior_workflow <- workflow() %>%
  add_model(rf_health_behavior) %>%
  add_recipe(rf_health_behavior_recipe)

# Fit model
rf_health_behavior_fit <- rf_health_behavior_workflow %>%
  fit(data = train)


# Predict on test set
pred_probs <- predict(rf_health_behavior_fit, test, type = "prob")
pred_class <- predict(rf_health_behavior_fit, test, type = "class")

# Combine everything
pred <- bind_cols(
  test %>% select(Cancer_Avoiders01),
  pred_class,
  pred_probs
)

# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoiders01, estimate = .pred_class)

# Confusion matrix
conf_mat_result <- conf_mat(pred, truth = Cancer_Avoiders01, estimate = .pred_class)

# ROC AUC
roc_auc_val <- roc_auc(pred, truth = Cancer_Avoiders01, .pred_1)

# Variable importance
rf_health_behavior_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)
```

```{r}
# Create performance metrics table
performance_df <- data.frame(
  Metric = metrics$.metric,
  Value = metrics$.estimate
)

kable(
  performance_df,
  caption = "Random Forest Classification Model Performance (Health Behavior Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Create confusion matrix table
conf_mat_tbl <- as.data.frame(conf_mat_result$table)
colnames(conf_mat_tbl) <- c("Truth", "Prediction", "Count")

conf_mat_wide <- tidyr::pivot_wider(
  conf_mat_tbl,
  names_from = Prediction,
  values_from = Count
)

kable(
  conf_mat_wide,
  caption = "Confusion Matrix",
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Add ROC AUC to a separate table
roc_df <- data.frame(
  Metric = "ROC AUC",
  Value = roc_auc_val$.estimate
)

kable(
  roc_df,
  caption = "ROC AUC Score",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```

-   High accuracy (0.946) simply reflects predicting majority class.

-   ROC AUC = 0.472, which means no better than random

-   The random forest model shows that `Cigarette_Smoking` is the most predict variable.

```{r}
#| label: logistic-regression-model-Cigarette_Smoking_num-predicting-Cancer_Avoiders01

smoking_cancer_logistic <- glm(Cancer_Avoiders01 ~ Cigarette_Smoking_num, data = health_behavior_data, family = binomial)

coef_df <- as.data.frame(summary(smoking_cancer_logistic)$coefficients)

knitr::kable(
  coef_df,
  caption = "Logistic Regression: Cigarette Smoking Predicting Cancer Avoiders",
  digits = 8
)
```

-   The logistic regression shows a statistically significant positive relationship between smokers and the probability of being a cancer avoider ($\beta$ = 0.62694, p = 1.05e-06). Smokers have odds of being a cancer avoider that are 1.87 times higher ($e^0.62694$ = 1.872) compared to non-smokers, representing an 87.2% increase in odds

-   However, the model shows modest improvement over the null model, with a 21.7 point reduction in deviance (null deviance = 2709.8, residual deviance = 2688.1), indicating that Facebook usage explains modest of the variation in cancer avoider.

```{r}
#| label: auc-cigarette
#| warning: false

library(pROC)

smoking_pred <- predict(smoking_cancer_logistic, type = "response")

roc_obj <- roc(health_behavior_data$Cancer_Avoiders01, smoking_pred)
auc(roc_obj)
```

-   ROC AUC = 0.55, which means slightly better than random

```{r}
ggplot(health_behavior_data, aes(x = factor(Cigarette_Smoking_num), fill = factor(Cancer_Avoiders01))) +
  geom_bar(position = "fill") +  # stacked bar normalized to proportion
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("0" = "lightgray", "1" = "red"),
                    labels = c("No", "Yes")) +
  labs(
    x = "Cigarette Smoking",
    y = "Percentage",
    fill = "Cancer Avoiders",
    title = "Distribution of Cancer Avoiders by Smoking Status"
  ) +
  theme_minimal()

```

-   The stacked bar chart displays the percentage distribution of cancer avoiders and non-avoiders across smoking status (0 = Non-smoker, 1 = Smoker). Smokers display a larger red segment at approximately 8-9% cancer avoiders. This visual pattern confirms the logistic regression finding, smokers are roughly twice as likely to be cancer avoiders compared to non-smokers.

-   However, it's important to note that despite this significant association, cancer avoiders remain a minority in both groups, with over 90% of both smokers and non-smokers classified as non-avoiders.

#### Other Model

```{r}
#| label: random-forest-other-model-class
#| fig-width: 15
#| fig-height: 4


library(tidymodels)
library(tidyr)
library(vip)

# need to drop NA to get accuracy
other_data <- selectdata %>%
  drop_na(
    Cancer_Avoiders01, Home_Ownership, Voter_Registration, Climate_Change_Belief,
    Mental_Health_of_Partner
  )


# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(other_data)), size = 0.7 * nrow(other_data))
train <- other_data[train_idx, ]
test  <- other_data[-train_idx, ]

# Fit random forest model (only for demo)
rf_other <- rand_forest(
  trees = 500
) %>%
  set_engine("ranger", importance = "impurity") %>% 
  set_mode("classification")

# Set recipe
rf_other_recipe <- recipe(Cancer_Avoiders01 ~  Home_Ownership + Voter_Registration + Climate_Change_Belief + Mental_Health_of_Partner,
  data = train) %>%
  step_other(all_nominal_predictors(), threshold = 0.05) %>%  # collapse rare levels <5%
  step_dummy(all_nominal_predictors()) 

# Create workflow
rf_other_workflow <- workflow() %>%
  add_model(rf_other) %>%
  add_recipe(rf_other_recipe)

# Fit model
rf_other_fit <- rf_other_workflow %>%
  fit(data = train)


# Predict on test set
pred_probs <- predict(rf_other_fit, test, type = "prob")
pred_class <- predict(rf_other_fit, test, type = "class")

# Combine everything
pred <- bind_cols(
  test %>% select(Cancer_Avoiders01),
  pred_class,
  pred_probs
)

# Model performance
metrics <- pred %>%
  metrics(truth = Cancer_Avoiders01, estimate = .pred_class)

# Confusion matrix
conf_mat_result <- conf_mat(pred, truth = Cancer_Avoiders01, estimate = .pred_class)

# ROC AUC
roc_auc_val <- roc_auc(pred, truth = Cancer_Avoiders01, .pred_1)

# Variable importance
rf_other_fit %>%
  extract_fit_parsnip() %>%
  vip(num_features = 10)
```

```{r}
# Create performance metrics table
performance_df <- data.frame(
  Metric = metrics$.metric,
  Value = metrics$.estimate
)

kable(
  performance_df,
  caption = "Random Forest Classification Model Performance (Other Predictors)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Create confusion matrix table
conf_mat_tbl <- as.data.frame(conf_mat_result$table)
colnames(conf_mat_tbl) <- c("Truth", "Prediction", "Count")

conf_mat_wide <- tidyr::pivot_wider(
  conf_mat_tbl,
  names_from = Prediction,
  values_from = Count
)

kable(
  conf_mat_wide,
  caption = "Confusion Matrix",
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")

# Add ROC AUC to a separate table
roc_df <- data.frame(
  Metric = "ROC AUC",
  Value = roc_auc_val$.estimate
)

kable(
  roc_df,
  caption = "ROC AUC Score",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```


-   High accuracy (0.96) simply reflects predicting majority class.

-   ROC AUC = 0.42, which means no better than random

-   The random forest model shows that `Climate_Change_Belief`(other) is the most predict variable.

```{r}
#| label: logistic-regression-model-Climate_Change_Belief-predicting-Cancer_Avoiders01
#| wrap: true
#| tidy: true

climate_cancer_logistic <- glm(Cancer_Avoiders01 ~ Climate_Change_Belief, data = other_data, family = binomial)

coef_df <- as.data.frame(summary(climate_cancer_logistic)$coefficients)

knitr::kable(
  coef_df,
  caption = "Logistic Regression: Climate Change Belief Predicting Cancer Avoiders",
  digits = 5
)

```

-   The logistic regression shows significant negative associations between accepting human caused climate change and the probability of being a cancer avoider. Individuals who strongly believe climate change is primarily caused by human activities have 71% lower odds (1 - $e^-1.236$ = 0.71) of being cancer avoiders ($\beta$ = -1.236, p = 1.87e-10) compared to climate dennier.

-   The model shows modest improvement over the null model with a 50.1 point deviance reduction (null deviance = 2781.4, residual deviance = 2731.3)

```{r}
#| label: auc-climate
#| warning: false

climate_pred <- predict(climate_cancer_logistic, type = "response")

roc_obj_voter <- roc(other_data$Cancer_Avoiders01, climate_pred)

auc(roc_obj_voter)
```

-   ROC AUC = 0.58, which means slightly better than random

```{r}
#| fig.width: 12
library(ggplot2)

ggplot(other_data, aes(x = Climate_Change_Belief, fill = factor(Cancer_Avoiders01))) +
  geom_bar(position = "fill") +  # stacked bars scaled to proportion
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("0" = "lightgray", "1" = "red"),
                    labels = c("No", "Yes")) +
  labs(
    x = "Climate Change Belief",
    y = "Percentage of Cancer Avoiders",
    fill = "Cancer Avoiders",
    title = "Cancer Avoiders by Climate Change Belief"
  ) +
  scale_x_discrete(labels = c("Strongly skeptical of claims about climate change and its link to human activities." = "Strongly disbelieve", 
                              "Somewhat skeptical about the impact of human activities on climate change, believing that climate change is a natural cycle." = "Somewhat disbelieve", 
                              "Uncertain about the causes and extent of climate change." = "Uncertain",
                              "No opinion on the matter." = "No opinion",
                              "Somewhat believe climate change is occurring and is influenced by human activities, but natural factors also play a significant role." = "Somewhat believe",
                              "Strongly believe climate change is occurring and is primarily caused by human activities." = "Strongly believe")) +
  theme_minimal(base_size = 14)

```

-   The stacked bar chart displays the percentage distribution of cancer avoiders across six climate change belief categories. This visual pattern demonstrates that individuals who deny about human-caused climate change are 3-4 times more likely to be cancer avoiders compared to those who strongly accept about human-caused climate change.

## Comprehensive Model with All Predictors

```{r}
library(tidymodels)
# Predictor variables (based on previous feature importance)
predictors <- c("MacArthur_Numeric", "Facebook_Usage_cat", 
                "Stress_TotalScore", "Cigarette_Smoking_num", "Climate_Change_Belief")

# Outcomes
outcome_numeric <- "Cancer_Avoidance_Mean"
outcome_binary <- "Cancer_Avoiders01"

# Select relevant columns
rf_data <- selectdata %>%
  dplyr::select(all_of(c(outcome_numeric, outcome_binary, predictors)))


```

```{r}
library(recipes)

recipe_numeric <- recipe(as.formula(
  paste(outcome_numeric, "~", paste(predictors, collapse = "+"))
), data = rf_data) %>%
  step_impute_median(all_numeric_predictors()) %>%
  step_impute_mode(all_nominal_predictors()) %>%
  step_dummy(all_nominal_predictors())

recipe_binary <- recipe(as.formula(
  paste(outcome_binary, "~", paste(predictors, collapse = "+"))
), data = rf_data) %>%
  step_impute_median(all_numeric_predictors()) %>%
  step_impute_mode(all_nominal_predictors()) %>%
  step_dummy(all_nominal_predictors())

# Regression RF
rf_numeric_spec <- rand_forest(
  mode = "regression",
  trees = 500,
  mtry = 3,
  min_n = 5
) %>%
  set_engine("ranger", importance = "impurity")  # calculate variable importance


# Classification RF
rf_binary_spec <- rand_forest(
  mode = "classification",
  trees = 500,
  mtry = 3,
  min_n = 5
) %>%
  set_engine("ranger", importance = "impurity")  # calculate variable importance

# Workflow for numeric outcome
workflow_numeric <- workflow() %>%
  add_model(rf_numeric_spec) %>%
  add_recipe(recipe_numeric)

# Workflow for binary outcome
workflow_binary <- workflow() %>%
  add_model(rf_binary_spec) %>%
  add_recipe(recipe_binary)

rf_numeric_fit <- workflow_numeric %>% fit(data = rf_data)
rf_binary_fit <- workflow_binary %>% fit(data = rf_data)


```

```{r}
#| fig-width: 15
#| fig-height: 4

library(ranger)

# Predictions for regression model
pred_numeric <- predict(rf_numeric_fit, rf_data) %>%
  bind_cols(rf_data %>% dplyr::select(Cancer_Avoidance_Mean))

# Model performance (MAE, RMSE, R^2)
metrics_numeric <- pred_numeric %>%
  metrics(truth = Cancer_Avoidance_Mean, estimate = .pred)

rmse_value <- pred_numeric %>%
  rmse(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)


mae_value <- pred_numeric %>%
  mae(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)

rsq_value <- pred_numeric %>%
  rsq(truth = Cancer_Avoidance_Mean, estimate = .pred) %>%
  pull(.estimate)

cor_numeric <- cor(pred_numeric$.pred, pred_numeric$Cancer_Avoidance_Mean)

# Numeric outcome Cancer_Avoidance_Mean
vip::vip(rf_numeric_fit$fit$fit)

# Create performance table
performance_df <- data.frame(
  Metric = c("RMSE", "MAE", "R-squared", "Correlation"),
  Value = c(rmse_value, mae_value, rsq_value, cor_numeric)
)

kable(
  performance_df,
  caption = "Comprehensive Random Forest Model Performance (Regression)",
  digits = 3,
  booktabs = TRUE,
  row.names = FALSE
) %>%
  kable_styling(latex_options = "hold_position")
```

-   RMSE (\~0.62) and MAE (\~0.52) are pretty close, meaning errors aren’t heavily dominated by outliers.

-   Correlation = 0.3 is low. But it imporve a little bit compare to all the regression models. It basically means that predicted values do not track the actual values at all. So the model is not capturing the pattern in the test data.

-   If predictors have very weak relationships with the outcome, the model will predicts values near the mean of the training set.

-   The random forest model shows that `Stress_TotalScore` is the most predict variable.

```{r}
#| fig-width: 15
#| fig-height: 4

# Predictions for classification model
pred_binary <- predict(rf_binary_fit, rf_data, type = "class") %>%
  bind_cols(predict(rf_binary_fit, rf_data, type = "prob")) %>%
  bind_cols(rf_data %>% dplyr::select(Cancer_Avoiders01))

metrics_binary <- pred_binary %>%
  metrics(truth = Cancer_Avoiders01, estimate = .pred_class)

conf_mat_val <- conf_mat(pred_binary, truth = Cancer_Avoiders01, estimate = .pred_class)

roc_auc_value <- roc_auc(pred_binary, truth = Cancer_Avoiders01, .pred_1)

# Binary outcome Cancer_Avoiders01
vip::vip(rf_binary_fit$fit$fit)
```

```{r}
metrics_table <- metrics_binary %>%
  select(.metric, .estimate) %>%
  rename(
    Metric = .metric,
    Estimate = .estimate
)

kable(
  metrics_table,
  caption = "Performance Metrics for Binary Outcome Model",
  digits = 3,
  escape = TRUE,
  booktabs = TRUE,
  row.names = FALSE
) %>%
kable_styling(latex_options = "hold_position")
```

-   High accuracy (0.95) simply reflects predicting majority class.

-   ROC AUC = 0.23, which means no better than random, and is worse than all the classification models.

-   The random forest model shows that `Stress_TotalScore` is the most predict variable.

# Discussion

This analysis examined predictors of cancer avoidance scores across demographic, media-related, health, behavioral, and other domains using regression, classification, and MARS models. Across all approaches, the findings point to the same conclusion is that many predictors reached significance, but none meaningfully predicted cancer avoidance.

Regression models explained less than 3% of variance ($R^2$ \< 0.03), and the full model reached only a modest correlation of 0.30. Classification models showed high accuracy due to class imbalance, but poor distinction, with ROC AUC values between 0.23 and 0.49, at below chance. The full classification model performed worst (AUC = 0.23), reinforcing that adding strong predictors does not improve performance.

MARS models found some nonlinear patterns, but the negative cross-validated $R^2$ shows they likely don’t hold up and reflect overfitting. Overall, the small predictive value suggests important factors are missing. The low number of cancer avoiders (4–5%) also makes prediction harder. The unexpected links with smoking and climate beliefs may come from unmeasured factors or differences in how people understand cancer avoidance.

## Limitations

-   Cross-sectional design: Because the data were collected at one point in time, we cannot tell what causes what. Other unmeasured factors may also affect the results.

-   Self-reported outcome: People may interpret “cancer avoidance” differently, which means the score may not fully reflect their real behaviors.

-   Small effect sizes: Many predictors were statistically significant but had tiny effects, likely because the sample size was large.

-   Inconsistent samples across models: Missing data led to different subsets being used, limiting direct comparisons.

-   Random forest issues: The variable importance results were not consistent with simpler models and may exaggerate weak patterns.

-   Limited generalizability: Results may not generalize to other groups or contexts.

## Future Directions

Future studies should: (1) use longitudinal data to better understand cause-and-effect; (2) explore factors that might explain or change the relationships found; (3) oversample cancer avoiders to address class imbalance; (4) include qualitative work to understand unexpected patterns; and (5) validate these findings in independent samples.

## Acknowledgments {.appendix}

We thank the following people and organizations for their guidance, support, and resources in this project:

-   **Dr. Shane McCarty** (Binghamton University) – Principal Investigator and mentor\
-   **Dr. Heather Orom** (University at Buffalo) – Principal Investigator and mentor\
-   **Dr. Kargin Vladislav** (Binghamton University) – Principal Investigator and mentor\
-   **Cloud Research** – Owner of the Health Avoiders dataset and provider of access to de-identified survey data


## Appendix A: Detailed Statistical Output {.appendix}
### Linear Regression

#### MacArthur Scale vs cancer-avoidance score

```{r}
#| label: MacArthur_Numeric-predicting-Cancer-Avoidance-Mean
#| echo: true

MacArthur_cancer_linear <- lm(Cancer_Avoidance_Mean ~ MacArthur_Numeric, data = demo_data)
summary(MacArthur_cancer_linear)
```

#### Age vs cancer-avoidance score

```{r}
#| label: AgeGroup-predicting-Cancer-Avoidance-Mean
#| echo: true
agegroup_cancer_linear <- lm(Cancer_Avoidance_Mean ~ AgeGroup, data = demo_data)
summary(age_cancer_linear)
```

#### Age vs cancer-avoidance score
```{r}
#| label: AgeBand-predicting-Cancer-Avoidance-Mean
#| echo: true
ageband_cancer_linear <- lm(Cancer_Avoidance_Mean ~ AgeBand, data = demo_data)
summary(age_cancer_linear)
```

#### Age over 35 vs cancer-avoidance score
```{r}
#| label: Age-over35-predicting-Cancer-Avoidance-Mean
#| echo: true

model_over35 <- lm(Cancer_Avoidance_Mean ~ Education_Level + Income + 
                   MacArthur_Numeric + Political_Party,
                   data = dplyr::filter(demo_data, AgeBand == "35+"))

summary(model_over35)
```

#### Influencer Following vs cancer-avoidance score

```{r}
#| label: Influencer_Following-predicting-Cancer_Avoidance_Mean
#| echo: true
influencer_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Influencer_Following, data = media_data)
summary(influencer_cancer_linear)
```

#### Stress Score vs cancer-avoidance score

```{r}
#| label: Stress_TotalScore-predicting-Cancer_Avoidance_Mean
#| echo: true
stress_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Stress_TotalScore, data = health_condition_data)
summary(stress_cancer_linear)
```

#### Smoking vs cancer-avoidance score

```{r}
#| label: Cigarette_Smoking_num-predicting-Cancer_Avoidance_Mean
#| echo: true
smoking_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Cigarette_Smoking_num, data = health_behavior_data)
summary(smoking_cancer_linear)
```

#### Climate Change Belief vs cancer-avoidance score

```{r}
#| label: Climate_Change_Belief-predicting-Cancer_Avoidance_Mean
#| echo: true
#| wrap: true
#| tidy: true

voter_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Climate_Change_Belief, data = other_data)
summary(voter_cancer_linear)

```

### Mars
#### Demographic model
```{r}
#| label: demographic-model

# need to drop NA to get accuracy
demo_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Ethnicity, Political_Party, Gender4, Job_Classification,
    Education_Level, Age, Income, Race, MacArthur_Numeric
  )

# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(demo_data)), size = 0.7 * nrow(demo_data))
train <- demo_data[train_idx, ]
test  <- demo_data[-train_idx, ]

# Fit MARS model
mars_model <- earth(
  Cancer_Avoidance_Mean ~ Ethnicity + Political_Party + Gender4 + Job_Classification +
    Education_Level + Age + Income + Race + MacArthur_Numeric,
  data = train
)

# Predict on test set
pred <- predict(mars_model, newdata = test)

summary(mars_model)

# Variable importance
evimp(mars_model)
```

#### Cross validation demographic
```{r}
#| label: cross-validation-demographic
#| echo: true
#| wrap: true
#| tidy: true

library(earth)

mars_model <- earth(
  Cancer_Avoidance_Mean ~ Ethnicity + Political_Party + Gender4 + Job_Classification +
    Education_Level + Age + Income + Race + MacArthur_Numeric,
  data = train,
  degree = 2,       # allow up to 2-way interactions
  nfold = 10,       # 10-fold CV
  keepxy = TRUE
)
summary(mars_model) 
summary(mars_model) %>% .$coefficients %>% head(10)


```

#### Health condition Model

```{r}
#| label: health-condition-model
#| echo: true

# need to drop NA to get accuracy
health_condition_data <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Stressful_Events_Recent, Current_Depression, Anxiety_Severity_num, PTSD5_Score,
    Health_Depression_Severity_num, Stress_TotalScore
  )

# split into training and testing sets
set.seed(123)
train_idx <- sample(seq_len(nrow(health_condition_data)), size = 0.7 * nrow(health_condition_data))
train <- health_condition_data[train_idx, ]
test  <- health_condition_data[-train_idx, ]

# Fit MARS model
mars_model_health_condition <- earth(
    Cancer_Avoidance_Mean ~ Stressful_Events_Recent + Current_Depression + Anxiety_Severity_num +
    PTSD5_Score +  Health_Depression_Severity_num + Stress_TotalScore,
  data = train
)

# Predict on test set
pred <- predict(mars_model_health_condition, newdata = test)

summary(mars_model_health_condition)

# Variable importance
evimp(mars_model_health_condition)
```

#### Cross validation on Health Condition

```{r}
#| label: health-condition
#| echo: true

library(earth)

mars_model_health_condition <- earth(
    Cancer_Avoidance_Mean ~ Stressful_Events_Recent + Current_Depression + Anxiety_Severity_num +
    PTSD5_Score +  Health_Depression_Severity_num + Stress_TotalScore,
  data = train,
  degree = 2,       # allow up to 2-way interactions
  nk = 100,
  nfold = 10,       # 10-fold CV
  keepxy = TRUE
)
summary(mars_model_health_condition)

pred <- predict(mars_model_health_condition, newdata = test)
rmse <- sqrt(mean((pred - test$Cancer_Avoidance_Mean)^2))
cor(pred, test$Cancer_Avoidance_Mean)


```

### Logistic Regression
#### MacArthur Scale vs cancer-avoidance score

```{r}
#| label: MacArthur_Numeric-predicting-Cancer-Avoiders01
#| echo: true

MacArthur_cancer_logistic <- glm(Cancer_Avoiders01 ~ MacArthur_Numeric, data = demo_data, family = binomial)
summary(MacArthur_cancer_logistic)
```


#### Facebook usage vs cancer-avoidance score
```{r}
#| label: Facebook_Usage_cat-predicting-Cancer_Avoiders01
#| echo: true

facebook_cancer_logistic <- glm(Cancer_Avoiders01 ~ Facebook_Usage_cat, data = media_data, family = binomial)
summary(facebook_cancer_logistic)
```

#### Stress score vs cancer-avoidance score
```{r}
#| label: Stress_TotalScore-predicting-Cancer_Avoiders01
#| echo: true

stress_cancer_logistic <- glm(Cancer_Avoiders01 ~ Stress_TotalScore, data = health_condition_data, family = binomial)
summary(stress_cancer_logistic)
```

#### Smoking vs cancer-avoidance score
```{r}
#| label: Cigarette_Smoking_num-predicting-Cancer_Avoiders01
#| echo: true

smoking_cancer_logistic <- glm(Cancer_Avoiders01 ~ Cigarette_Smoking_num, data = health_behavior_data, family = binomial)
summary(smoking_cancer_logistic)
```

#### Climate Change vs cancer-avoidance score

```{r}
#| label: Climate_Change_Belief-predicting-Cancer_Avoiders01
#| echo: true
#| wrap: true
#| tidy: true

climate_cancer_logistic <- glm(Cancer_Avoiders01 ~ Climate_Change_Belief, data = other_data, family = binomial)
summary(climate_cancer_logistic)
```



