# ========================================================
# FROM A RESEARCH QUESTION TO AN ANALYSIS PLAN
# Canadian Social Connection Survey (CSCS), 2021 wave
# ========================================================
# Question: among adults in Canada, is living with a dog
# associated with lower loneliness?
# This script checks feasibility, draws the causal model
# with dagitty, and writes the analysis plan as code.
# It does not test the main association: that comes
# only after the plan is written down.
# ========================================================


# ========================================================
# SECTION 1: LOAD THE PACKAGE AND THE DATA
# ========================================================

# install.packages("dagitty")   # run once on each computer
library(dagitty)

github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData"
load(url(github_url))

# Keep the 2021 wave, which gives one row per person
data <- data[data$SURVEY_collection_year == 2021, ]
dim(data)


# ========================================================
# SECTION 2: FEASIBILITY CHECKS
# ========================================================

# --- 2.1 The exposure: how many dogs each person lives with
table(data$GEO_housing_live_with_dogs, useNA = "ifany")

# --- 2.2 The outcome: the UCLA loneliness score, 3 to 9
summary(data$LONELY_ucla_loneliness_scale_score)

# --- 2.3 Do implausible answers travel together?
many_dogs <- data$GEO_housing_live_with_dogs >= 5
big_house <- data$GEO_housing_household_size >= 15
table(many_dogs, big_house)


# ========================================================
# SECTION 3: PREPARE THE PLANNING VARIABLES
# ========================================================

plan <- data.frame(
  dogs       = data$GEO_housing_live_with_dogs,
  loneliness = data$LONELY_ucla_loneliness_scale_score,
  age        = data$DEMO_age,
  household  = data$GEO_housing_household_size,
  children   = data$GEO_housing_live_with_children,
  fish       = data$GEO_housing_live_with_fish)

# --- 3.1 Exposure: lives with at least one dog (1) or none (0)
plan$dog <- NA
plan$dog[plan$dogs == 0] <- 0
plan$dog[plan$dogs >= 1] <- 1

# --- 3.2 Confounders from the causal diagram
plan$alone <- NA                       # lives alone
plan$alone[plan$household == 0] <- 1
plan$alone[plan$household >= 1] <- 0

plan$kids <- NA                        # lives with children
plan$kids[plan$children == 0] <- 0
plan$kids[plan$children >= 1] <- 1

plan$gender <- NA                      # a non-answer stays missing
plan$gender[data$DEMO_gender == "Man"] <- "Man"
plan$gender[data$DEMO_gender == "Woman"] <- "Woman"
plan$gender[data$DEMO_gender == "Non-binary"] <- "Non-binary"

plan$disability <- NA                  # chronic health problem or disability
plan$disability[data$DEMO_identity_disability == "Not selected"] <- 0
plan$disability[data$DEMO_identity_disability ==
                "People with chronic health problems or disabilities"] <- 1

income_level <- as.numeric(data$DEMO_household_income)   # 1 to 26 in order
plan$income <- NA
plan$income[income_level <= 8] <- "1: under $40k"
plan$income[income_level >= 9  & income_level <= 13] <- "2: $40k to $79k"
plan$income[income_level >= 14 & income_level <= 17] <- "3: $80k to $119k"
plan$income[income_level >= 18 & income_level <= 26] <- "4: $120k or more"

# --- 3.3 How much is missing, and how many complete cases?
key_vars <- c("dog", "loneliness", "age", "gender", "income",
              "alone", "kids", "disability")
colSums(is.na(plan[, key_vars]))
sum(complete.cases(plan[, key_vars]))

# --- 3.4 Small cells: dog ownership by living alone
table(dog = plan$dog, alone = plan$alone)


# ========================================================
# SECTION 4: THE CAUSAL DIAGRAM (DAG) IN DAGITTY
# ========================================================

dag_text <- '
  dog        [exposure, pos="0,0"]
  loneliness [outcome,  pos="4,0"]
  age        [pos="1,-1.6"]
  gender     [pos="0.2,-1"]
  income     [pos="2,-1.9"]
  alone      [pos="3,-1.6"]
  kids       [pos="3.8,-1"]
  disability [pos="2,-0.8"]
  activity   [pos="2,1.1"]
  age -> dog          age -> loneliness    age -> income
  age -> alone        age -> kids          age -> disability
  gender -> dog       gender -> loneliness
  income -> dog       income -> loneliness
  alone -> dog        alone -> loneliness  kids -> alone
  kids -> dog         kids -> loneliness
  disability -> dog   disability -> loneliness
  disability -> activity
  dog -> activity     activity -> loneliness
  dog -> loneliness
'
dag <- dagitty(paste("dag {", dag_text, "}"))
plot(dag)

# --- 4.1 Classify the variables
common_causes <- intersect(ancestors(dag, "dog"), ancestors(dag, "loneliness"))
setdiff(common_causes, "dog")                       # confounders
on_path <- intersect(descendants(dag, "dog"), ancestors(dag, "loneliness"))
setdiff(on_path, c("dog", "loneliness"))            # mediators

# --- 4.2 What must be adjusted for?
adjustmentSets(dag)                       # the total effect
adjustmentSets(dag, effect = "direct")    # the effect not through activity

# --- 4.3 A tempting extra variable: depressive symptoms
dag_d <- dagitty(paste("dag {", dag_text,
  "depression  activity -> depression  loneliness -> depression }"))
descendants(dag_d, "loneliness")
adjustmentSets(dag_d)

# --- 4.4 An unmeasured confounder: personality
dag_u <- dagitty(paste("dag {", dag_text,
  "personality [latent]  personality -> dog  personality -> loneliness }"))
adjustmentSets(dag_u)

# --- 4.5 What the DAG implies we should see in the data
head(impliedConditionalIndependencies(dag), 4)


# ========================================================
# SECTION 5: THE ANALYSIS PLAN, WRITTEN AS CODE
# ========================================================

# Primary model: the adjustment set from section 4.2
primary <- loneliness ~ dog + age + gender + income +
  alone + kids + disability
primary

# Secondary: does the association differ for people who live alone?
by_alone <- loneliness ~ dog * alone + age + gender + income +
  kids + disability

# Sensitivity: a negative control exposure (living with fish)
plan$fish_any <- NA
plan$fish_any[plan$fish == 0] <- 0
plan$fish_any[plan$fish >= 1] <- 1
neg_control <- loneliness ~ fish_any + age + gender + income +
  alone + kids + disability

# --- 5.1 The decisions each sensitivity analysis will check
sensitivity <- data.frame(
  decision    = c("Blank dog counts", "Implausible answers",
                  "Outcome coding", "Adjustment set",
                  "Unmeasured confounding", "Model"),
  primary     = c("Missing", "Kept", "Score, 3 to 9",
                  "Minimal set", "Assumed absent", "Linear regression"),
  alternative = c("Counted as 0", "Excluded", "Lonely (6 or more)",
                  "Extended set", "E-value; fish", "Ordinal logistic"))
sensitivity