attempting quick Cox survival model in python and R

Author: Elye Bliss

notebooks/2025_09_16_child_mortality_investigation.py

@@ -3,12 +3,14 @@
 import pandas as pd
 import matplotlib.pyplot as plt
 import matplotlib.ticker as mticker
+from lifelines import CoxPHFitter  # for Cox survival models
 
 
 DATA_PATH = "/mnt/team/rapidresponse/pub/population/modeling/climate_malnutrition/child_mortality/training_data/2025_09_15.01/data.parquet"
 
 df = pd.read_parquet(DATA_PATH)
 
+## 1. Get basic info about data
 df["line_id"] = df["line_id"].astype(int)
 
 df["indv_id"] = df[["nid", "psu", "hh_id", "line_id"]].astype(str).agg("_".join, axis=1)
@@ -48,6 +50,7 @@
     f"who did not die were under 5 at survey ({under5_at_survey/(over5_at_survey + under5_at_survey):.1%})"
 )
 
+## 2. Make scatterplots and heatmaps based on raw data
 # Aggregate data
 agg_df = (
     df.groupby(["nid", "ihme_loc_id", "int_year"], as_index=False)
@@ -98,3 +101,28 @@
     fmt=".2f",
     cmap="YlOrBr",
 )
+
+## 3. Make simple model
+df.rename(columns={"ldipc_weighted_no_match": "consumption"}, inplace=True)
+event_col = "child_mortality"
+id_col = "indv_id"
+time_col = "age_month_at_year_end"
+covariate_cols = [
+    "consumption",
+    "ihme_loc_id",  # r.e. not yet supported in lifelines package
+    "mean_temperature",
+    # "total_precipitation",
+    # "relative_humidity",
+    # "mean_high_temperature",
+    # "mean_low_temperature",
+    # "precipitation_days",
+    "days_over_30C",
+    # "days_over_26C",
+]
+df_model_data = df[[event_col, time_col] + covariate_cols]
+
+cph = CoxPHFitter()
+cph.fit(df_model_data, duration_col=time_col, event_col=event_col)
+cph.print_summary()
+
+## 4. Plot data sources by country

notebooks/2025_09_16_child_mortality_investigation_cox_model.R

@@ -0,0 +1,75 @@
+################################################################################
+# DESCRIPTION: Test R's coxme package to perform a Cox proportional hazard model
+# with mixed effects (r.e. on location, not yet available in Python packages)
+# PROJECT: Climate nutrition
+# DATE: 2025-09-16
+################################################################################
+
+#==============================================================================
+# SECTION 0: PACKAGE LOADING AND ENVIRONMENT SETUP
+#==============================================================================
+# Clear workspace
+rm(list = ls())
+
+# Username is pulled automatically
+username <- Sys.info()[["user"]]
+if (Sys.info()["sysname"] == "Linux") {
+  j <- "/home/j/"
+  h <- paste0("/homes/", username, "/")
+  r <- "/mnt/"
+  l <-"/ihme/limited_use/"
+} else {
+  j <- "J:/"
+  h <- "H:/"
+  r <- "R:/"
+  l <- "L:/"
+}
+
+install.packages('coxme',lib = "/homes/elyeb/rlibs") # for survival analysis with mixed effects
+library(coxme,lib.loc = "/homes/elyeb/rlibs") 
+library(data.table)
+library(arrow) # to read parquet
+
+options(scipen = 999) # turn off scientific notation
+ 
+#==============================================================================
+# SECTION 1: DATA LOADING AND PREPROCESSING
+#==============================================================================
+
+df <- read_parquet("/mnt/team/rapidresponse/pub/population/modeling/climate_malnutrition/child_mortality/training_data/2025_09_15.01/data.parquet")
+
+df <- data.table(df)
+
+# flip child_alive so 1 = died, 0 = alive for easier interpretation
+df[,child_mortality := 1-child_alive]
+  
+setnames(df,old="ldipc_weighted_no_match",new="consumption")
+
+df_model <- df[,.(child_mortality,age_month_at_year_end,sex_id,ihme_loc_id,consumption,mean_temperature,days_over_30C)]
+
+#==============================================================================
+# SECTION 2: MAIN PROCESSING
+#==============================================================================
+
+fit <- coxme(Surv(age_month_at_year_end, child_mortality) ~ consumption + mean_temperature + days_over_30C + sex_id + (1|ihme_loc_id), data = df_model)
+summary(fit)
+# Mixed effects coxme model
+# Formula: Surv(age_month_at_year_end, child_mortality) ~ consumption +      mean_temperature + days_over_30C + sex_id + (1 | ihme_loc_id) 
+# Data: df_model 
+# 
+# events, n = 42130, 4893786
+# 
+# Random effects:
+#   group  variable       sd variance
+# 1 ihme_loc_id Intercept 1.135688 1.289787
+# Chisq    df p   AIC   BIC
+# Integrated loglik 16118  5.00 0 16108 16065
+# Penalized loglik 16241 18.78 0 16204 16041
+# 
+# Fixed effects:
+#   coef    exp(coef)     se(coef)      z                   p
+# consumption      -0.000229703  0.999770323  0.000004872 -47.15 <0.0000000000000002
+# mean_temperature -0.015280408  0.984835745  0.001395616 -10.95 <0.0000000000000002
+# days_over_30C     0.003450759  1.003456720  0.000152407  22.64 <0.0000000000000002
+# sex_id           -0.017863722  0.982294889  0.009749483  -1.83              0.0669
+