---
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
output: html_document
editor: visual
format:
  html:
    toc: true
    code-fold: true
    code-summary: "Show the code"
    comments:
      hypothesis: true
    fig-caption: true
    fig-numbering: true
bibliography: references.bib
---

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

## Review of Machine Learning Methods

Machine learning (ML) methods were applied to identify patterns associated with health information avoidance. Each model differs in complexity, interpretability, and suitability for different types of data structures. 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 |

In summary, Linear / Logistic Regression provides transparency and simplicity, Random Forest offers high predictive power for complex data, and MARS strikes a balance between flexibility and interpretability.

So in our study, 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 under IRB approval managed by Dr. Shane McCarty. All analyses were conducted in R, and reproducibility was ensured through a version-controlled Git repository and 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 — `Ethnicity`, `Political_Party`, `Gender4`, `Job_Classification`, `Education_Level`, `Age`, `Income`, `Race`, and `MacArthur_Numeric`.

2.  Media Use Model — `Social_Media_Usage`, `AI_Use`, `Video_Games_Hours`, `Listening_Podcasts`, `Facebook_Usage_num`, `TikTok_Use`, `X_Twitter_Usage`, `Social_Media_type`, and `Influencer_Following`.

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

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

5.  Other Factors Model — `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.

## IMPORT

```{r}
#| label: load-tidymodels
#| message: false
#| warning: false
library(tidymodels)
library(ranger)
```

```{r}
#| label: load-Health-Avoiders
#| message: false
#| warning: false
library(readr)

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

## TRANSFORM

### Predictors

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

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+")),
    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
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: distribution-of-age

#  (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}
#| 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 <= 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 = 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,
    ),
    Supplement_Consumption_Reason_num = 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
    ),
    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_
    )),
    Facebook_Usage_num = 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
    )
)
```

### Health Avoidance

```{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),
    
    # 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")

```

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

#### Matrix

```{r}
#| label: correlation-of-cancer-mean-with-main-predictors
#| echo: true
#| 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
full_cor_matrix_clean <- cor(numeric_vars_clean, use = "complete.obs")

# Extract correlations with Cancer_Avoidance_Mean
avoidance_corrs <- full_cor_matrix_clean["Cancer_Avoidance_Mean", , drop = FALSE]

# Drop self-correlation
avoidance_corrs <- avoidance_corrs[, colnames(avoidance_corrs) != "Cancer_Avoidance_Mean", drop = FALSE]

# Convert to table
avoidance_corrs_df <- data.frame(
  Variable = colnames(avoidance_corrs),
  Correlation = round(as.numeric(avoidance_corrs), 2)
)

# Print table
knitr::kable(
  avoidance_corrs_df,
  caption = "Correlations of Cancer_Avoidance_Mean with main predictors"
)

```

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

library(ggplot2)

ggplot(avoidance_corrs_df, aes(x = reorder(Variable, 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)

```

## REGRESSION MODELS

### 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
selectdata <- 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(selectdata)), size = 0.7 * nrow(selectdata))
train <- selectdata[train_idx, ]
test  <- selectdata[-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 + Age + 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)
metrics

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

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

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

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

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

-   Correlation = 0.128 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 Age is the most predict variable.

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

-   Age coefficient = 0.0015164 meaning for every one year increase in age, the predicted avoidance score goes up by 0.0015164 units.

-   Since p-value = 0.009154, so is conventionally significant at the 0.05 level.

-   But look at the scatter plot, the trend suggests a slight nonlinear pattern, avoidance rises slowly with age in younger adults, peaks, then decreases slightly in the oldest participants. Despite statistical significance (p = 0.009), the practical effect is very small, as most points cluster near the mean and the trend line is almost flat.

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

age_cancer_logistic <- glm(Cancer_Avoiders01 ~ Age, data = selectdata, family = binomial)
summary(age_cancer_logistic)
```

-   Age coefficient = 0.017472 meaning for every one year increase in age, the predicted avoidance score goes up by 0.017472 units.

-   Since p-value = 5.31e-05, so is conventionally significant at the 0.05 level. Which meaning older participants tend to have slightly higher avoidance scores.

-   But look at the scatter plot, the trend suggests a slight nonlinear pattern, avoidance rises slowly with age in younger adults, peaks, then decreases slightly in the oldest participants. Despite statistical significance (p = 5.31e-05), the practical effect is very small, as most points cluster near the mean and the trend line is almost flat.

```{r}
#| label: scatter-plot-Age-predicting-Cancer-Avoidance-Mean

library(ggplot2)

ggplot(data = selectdata, aes(x = Age, y = Cancer_Avoidance_Mean)) +
  geom_point(alpha = 0.6, position = position_jitter(width=20), size = 0.5) +  
  geom_smooth(method = "loess", se = TRUE, color = "red") +  
  labs(x = "Age", 
       y = "Cancer Avoidance Mean",
       title = "Relationship between Age and Cancer Avoidance") +
  scale_x_continuous(limits = c(min(selectdata$Age), 
                              max(selectdata$Age))) +
  theme_minimal()
```

#### 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
selectdata <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, Social_Media_Usage, AI_Use, Video_Games_Hours, Listening_Podcasts,
    Facebook_Usage_num, 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(selectdata)), size = 0.7 * nrow(selectdata))
train <- selectdata[train_idx, ]
test  <- selectdata[-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_num + 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)
metrics

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

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

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

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

```

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

Correlation = -0.0024 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 Facebook Usage is the most predictive variable.

```{r}
#| label: linear-regression-model-Facebook_Usage_num-predicting-Cancer_Avoidance_Mean
facebook_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Facebook_Usage_num, data = selectdata)
summary(facebook_cancer_linear)
```

-   Facebook Usage coefficient = 0.028 meaning for every one 5 hours increase in facebook usage, the predicted avoidance score goes up by 0.028 units.

-   Since p-value = 0.1211, so is not conventionally significant at the 0.05 level.

-   But look at the scatter plot,

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

facebook_cancer_logistic <- glm(Cancer_Avoiders01 ~ Facebook_Usage_num, data = selectdata, family = binomial)
summary(facebook_cancer_logistic)
```

-   For Facebook_Usage_num, there is not a conventionally significant at the 0.05 level.

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

library(ggplot2)

ggplot(data = selectdata, aes(x = Facebook_Usage_num, y = Cancer_Avoidance_Mean)) +
  geom_point() +
  labs(x = "Facebook Usage", 
       y = "Cancer Avoidance Mean",
       title = "Relationship between Facebook Usage and Cancer Avoidance") +
  theme_minimal()

ggplot(data = selectdata, aes(x = Facebook_Usage_num, y = Cancer_Avoidance_Mean)) +
  geom_point(alpha = 0.6, position = position_jitter(width=0.5), size = 0.5) +  
  geom_smooth(method = "loess", se = TRUE, color = "red") +  
  labs(x = "Facebook Usage", 
       y = "Cancer Avoidance Mean",
       title = "Relationship between Facebook Usage and Cancer Avoidance") +
  scale_x_continuous(limits = c(min(selectdata$Facebook_Usage_num), 
                              max(selectdata$Facebook_Usage_num))) +
  theme_minimal()
```

#### Health Condition Model

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


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


# need to drop NA to get accuracy
selectdata <- 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(selectdata)), size = 0.7 * nrow(selectdata))
train <- selectdata[train_idx, ]
test  <- selectdata[-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)
metrics

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

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

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

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

```

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

Correlation = 0.054 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 = selectdata)
summary(stress_cancer_linear)
```

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

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

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

library(ggplot2)

ggplot(data = selectdata, aes(x = Stress_TotalScore, y = Cancer_Avoidance_Mean)) +
  geom_point() +
  labs(x = "Stress_TotalScore", 
       y = "Cancer Avoidance Mean",
       title = "Relationship between Stress_TotalScore and Cancer Avoidance") +
  theme_minimal()

ggplot(data = selectdata, aes(x = Stress_TotalScore, y = Cancer_Avoidance_Mean)) +
  geom_point(alpha = 0.6, position = position_jitter(width=5), size = 0.5) +  
  geom_smooth(method = "loess", se = TRUE, color = "red") +  # Adds a smooth curve
  labs(x = "Stress_TotalScore", 
       y = "Cancer Avoidance Mean",
       title = "Relationship between Stress_TotalScore and Cancer Avoidance") +
  scale_x_continuous(limits = c(min(selectdata$Stress_TotalScore), 
                              max(selectdata$Stress_TotalScore))) +
  theme_minimal()
```

#### Health Behavior Model

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


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


# need to drop NA to get accuracy
selectdata <- 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(selectdata)), size = 0.7 * nrow(selectdata))
train <- selectdata[train_idx, ]
test  <- selectdata[-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)
metrics

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

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

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

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

```

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

-   Correlation = 0.078 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 = selectdata)
summary(smoking_cancer_linear)
```

-   Cigarette Smoking coefficient = 0.20695 meaning for every one unit increase in cigarette smoking, the predicted avoidance score goes up by 0.20695 units. So, individuals who smoke more tend to have slightly higher cancer avoidance scores.

-   Since p-value = 2.882e-08, so is conventionally significant at the 0.05 level.

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

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

-   Cigarette Smoking coefficient = 0.8318 meaning for every one unit increase in cigarette smoking, the predicted avoidance score goes up by 0.8318 units. So, individuals who smoke more tend to have slightly higher cancer avoidance scores.

-   Since p-value = 0.000345, so is conventionally significant at the 0.05 level.

-   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 or express greater concern about cancer in general.

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

library(ggplot2)

ggplot(selectdata, aes(x = factor(Cigarette_Smoking_num), y = Cancer_Avoidance_Mean)) +
  geom_boxplot() +
  geom_jitter(width = 0.1, alpha = 0.4, size = 0.7) +
  labs(
    x = "Cigarette Smoking (0 = No, 1 = Yes)",
    y = "Cancer Avoidance Mean",
    title = "Cancer Avoidance by Cigarette Smoking Status"
  ) +
  theme_minimal()

```

#### 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
selectdata <- selectdata %>%
  drop_na(
    Cancer_Avoidance_Mean, 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(selectdata)), size = 0.7 * nrow(selectdata))
train <- selectdata[train_idx, ]
test  <- selectdata[-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 ~ 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 <- predict(rf_other_fit, test) %>%
  bind_cols(test %>% select(Cancer_Avoidance_Mean))


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

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

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

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

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

```

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

Correlation = 0.107 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 Voter Registration (No, but I'm eligiable to vote) is the most predict variable.

```{r}
#| label: linear-regression-model-Voter_Registration-predicting-Cancer_Avoidance_Mean
voter_cancer_linear <- lm(Cancer_Avoidance_Mean ~ Voter_Registration, data = selectdata)
summary(voter_cancer_linear)

```

-   On average, those eligible but not registered have a cancer avoidance mean 0.18343 points higher than registered voters, and this difference is statistically significant (p = 0.0018).

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

voter_cancer_logistic <- glm(Cancer_Avoiders01 ~ Voter_Registration, data = selectdata, family = binomial)
summary(voter_cancer_logistic)
```

-   On average, those eligible but not registered have a cancer avoidance mean 0.8903 points higher than registered voters, and this difference is statistically significant (p = 0.0094).

-   Look at the box plot, There appears to be a modest association between voter registration status and cancer avoidance. Eligible non-registrants show the highest cancer avoidance scores, while registered voters and those ineligible to vote show similar, lower levels. Registered voters shows a heavily concentrated distribution in the lower ranges with many overlapping points.

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

library(ggplot2)

ggplot(data = selectdata, aes(x = Voter_Registration, y = Cancer_Avoidance_Mean)) +
  geom_boxplot() +  
  geom_jitter(width = 0.1, alpha = 0.4, size = 0.7) +
  labs(x = "Voter_Registration", 
       y = "Cancer Avoidance Mean",
       title = "Relationship between Voter_Registration and Cancer Avoidance") +
  theme_minimal()
```

```{r}
#| 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

# need to drop NA to get accuracy
selectdata <- 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(selectdata)), size = 0.7 * nrow(selectdata))
train <- selectdata[train_idx, ]
test  <- selectdata[-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)
```

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

`MacArthur_Numeric` is about 59% as important.

`Education_Level` is weak (\~12%).

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

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)


```

Interactions appear as products:

EthnicityPrefer not to say \* Political_PartyDemocrat = an interaction term between ethnicity and political party.

predictive power is weak: training R-square \~ 0.07, CV R-square: -0.076 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_ClassificationBlue Collar and MacArthur_Numeric \> 3.

#### Health condition

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

# need to drop NA to get accuracy
selectdata <- 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(selectdata)), size = 0.7 * nrow(selectdata))
train <- selectdata[train_idx, ]
test  <- selectdata[-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)
```

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

`Stress_TotalScore` is about 34% as important.

`Anxiety_Severity_num` is weak (\~12%).

```{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
)
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)


```

h(PTSD5_Score-2): When PTSD score \> 2, Cancer Avoidance decreases by 0.51.

h(Stress_TotalScore - 7): When stress score \> 7, avoidance decreases by 0.16.

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.

CV R-square: -0.023. This suggests health condition variables don’t explain Cancer_Avoidance_Mean ve

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

```{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.02), indicating that although weak nonlinear relationships exist, the overall predictive power of the model remains low.
