/******************************************************************************
 * Programme : Causal Effect Estimation by Doubly Robust Methods with PROC CAEFFECT
 * Reference : CAUSALE10B
 * Source    : https://www.wearecas.eu/de/sampleCode/CAUSALE10B
 ******************************************************************************/

/* --- BLOC 1 --- */
/* Création d'une session CAS si elle n'existe pas */
cas;

/* Assurez-vous que la caslib mylib est disponible et mappée */
caslib _all_ assign;

/* Création de la table de données SmokingWeight */
data mylib.SmokingWeight;
   infile datalines;
   input Sex $ Race $ Education $ Exercise $ Activity $ Age YearsSmoke PerDay Quit Change;
   datalines;
Female White HighSchool Light Active 30 10 15 0 2.5
Male Black College Moderate Sedentary 45 20 20 1 5.0
Female Asian GradSchool Heavy Active 55 30 25 1 7.2
Male White HighSchool Light Sedentary 35 12 10 0 1.8
Female Black College Moderate Active 40 18 18 1 4.5
Male Asian GradSchool Heavy Sedentary 60 25 22 0 3.1
Female White HighSchool Light Active 28 8 12 1 6.0
Male Black College Moderate Sedentary 50 22 19 0 2.0
Female Asian GradSchool Heavy Active 48 28 24 1 8.0
Male White HighSchool Light Sedentary 32 10 11 0 1.5
;
run;

/* 1. Modélisation de la variable de traitement 'Quit' avec PROC LOGSELECT */
proc logselect data=mylib.SmokingWeight;
   class Sex Race Education Exercise Activity;
   model Quit(Event='1') = Sex Race Education
         Exercise Activity Age YearsSmoke PerDay;
   output out=mylib.swDREstData pred=pTrt copyvars=(_ALL_);
run;

/* Calcul de la probabilité de la condition de contrôle (Quit=0) */
data mylib.swDREstData;
   set mylib.swDREstData;
   pCnt = 1 - pTrt;
run;

/* 2. Modélisation de la variable de résultat 'Change' avec PROC BART */
proc bart data=mylib.swDREstData nTree=100 nMC=200 seed=2156;
   class Sex Race Education Exercise Activity Quit;
   model Change = Sex Race Education Exercise Quit Activity
         Age YearsSmoke PerDay;
   store out=mylib.bartOutMod;
run;

/* 3. Estimation des effets causaux avec PROC CAEFFECT (AIPW) */
proc caeffect data=mylib.swDREstData method=aipw;
   treatvar Quit;
   outcomevar Change;
   outcomemodel restore=mylib.bartOutMod predName=P_Change;
   pom treatLev=1 treatProb=pTrt;
   pom treatLev=0 treatProb=pCnt;
   difference evtLev=1;
run;


/* --- BLOC 2 --- */
/* Création d'une session CAS si elle n'existe pas */
cas;

/* Assurez-vous que la caslib mylib est disponible et mappée */
caslib _all_ assign;

/* Création de la table de données SmokingWeight */
data mylib.SmokingWeight;
   infile datalines;
   input Sex $ Race $ Education $ Exercise $ Activity $ Age YearsSmoke PerDay Quit Change;
   datalines;
Female White HighSchool Light Active 30 10 15 0 2.5
Male Black College Moderate Sedentary 45 20 20 1 5.0
Female Asian GradSchool Heavy Active 55 30 25 1 7.2
Male White HighSchool Light Sedentary 35 12 10 0 1.8
Female Black College Moderate Active 40 18 18 1 4.5
Male Asian GradSchool Heavy Sedentary 60 25 22 0 3.1
Female White HighSchool Light Active 28 8 12 1 6.0
Male Black College Moderate Sedentary 50 22 19 0 2.0
Female Asian GradSchool Heavy Active 48 28 24 1 8.0
Male White HighSchool Light Sedentary 32 10 11 0 1.5
;
run;

/* 1. Modélisation de la variable de traitement 'Quit' avec PROC LOGSELECT */
proc logselect data=mylib.SmokingWeight;
   class Sex Race Education Exercise Activity;
   model Quit(Event='1') = Sex Race Education
         Exercise Activity Age YearsSmoke PerDay;
   output out=mylib.swDREstData pred=pTrt copyvars=(_ALL_);
run;

/* Calcul de la probabilité de la condition de contrôle (Quit=0) */
data mylib.swDREstData;
   set mylib.swDREstData;
   pCnt = 1 - pTrt;
run;

/* 2. Modélisation de la variable de résultat 'Change' avec PROC BART */
proc bart data=mylib.swDREstData nTree=100 nMC=200 seed=2156;
   class Sex Race Education Exercise Activity Quit;
   model Change = Sex Race Education Exercise Quit Activity
         Age YearsSmoke PerDay;
   store out=mylib.bartOutMod;
run;

/* 3. Estimation des effets causaux avec PROC CAEFFECT (TMLE) */
proc caeffect data=mylib.swDREstData method=tmle;
   treatvar Quit;
   outcomevar Change;
   outcomemodel restore=mylib.bartOutMod predName=P_Change;
   pom treatLev=1 treatProb=pTrt;
   pom treatLev=0 treatProb=pCnt;
   difference evtLev=1;
run;


/* --- BLOC 3 --- */
/* Création d'une session CAS si elle n'existe pas */
cas;

/* Assurez-vous que la caslib mylib est disponible et mappée */
caslib _all_ assign;

/* Création de la table de données SmokingWeight */
data mylib.SmokingWeight;
   infile datalines;
   input Sex $ Race $ Education $ Exercise $ Activity $ Age YearsSmoke PerDay Quit Change;
   datalines;
Female White HighSchool Light Active 30 10 15 0 2.5
Male Black College Moderate Sedentary 45 20 20 1 5.0
Female Asian GradSchool Heavy Active 55 30 25 1 7.2
Male White HighSchool Light Sedentary 35 12 10 0 1.8
Female Black College Moderate Active 40 18 18 1 4.5
Male Asian GradSchool Heavy Sedentary 60 25 22 0 3.1
Female White HighSchool Light Active 28 8 12 1 6.0
Male Black College Moderate Sedentary 50 22 19 0 2.0
Female Asian GradSchool Heavy Active 48 28 24 1 8.0
Male White HighSchool Light Sedentary 32 10 11 0 1.5
;
run;

/* 1. Modélisation de la variable de traitement 'Quit' avec PROC LOGSELECT */
proc logselect data=mylib.SmokingWeight;
   class Sex Race Education Exercise Activity;
   model Quit(Event='1') = Sex Race Education
         Exercise Activity Age YearsSmoke PerDay;
   output out=mylib.swDREstData pred=pTrt copyvars=(_ALL_);
run;

/* Calcul de la probabilité de la condition de contrôle (Quit=0) */
data mylib.swDREstData;
   set mylib.swDREstData;
   pCnt = 1 - pTrt;
run;

/* 2. Modélisation de la variable de résultat 'Change' avec PROC BART */
proc bart data=mylib.swDREstData nTree=100 nMC=200 seed=2156;
   class Sex Race Education Exercise Activity Quit;
   model Change = Sex Race Education Exercise Quit Activity
         Age YearsSmoke PerDay;
   store out=mylib.bartOutMod;
run;

/* 3. Estimation des effets causaux avec PROC CAEFFECT (AIPW) par sous-groupe */
proc caeffect data=mylib.swDREstData method=aipw;
   class Sex; /* Ajout d'une variable de classification pour l'analyse par groupe */
   treatvar Quit;
   outcomevar Change;
   outcomemodel restore=mylib.bartOutMod predName=P_Change;
   pom treatLev=1 treatProb=pTrt;
   pom treatLev=0 treatProb=pCnt;
   difference evtLev=1;
   by Sex; /* Analyse séparée par le sexe */
   ods output Estimates=mylib.aipwEstimates OutDiff=mylib.aipwDifferences; /* Sauvegarde des résultats */
run;


/* --- BLOC 4 --- */
/* Création d'une session CAS si elle n'existe pas */
cas;

/* Assurez-vous que la caslib mylib est disponible et mappée */
caslib _all_ assign;

/* Nettoyage des tables existantes pour s'assurer d'un redémarrage propre */
proc casutil;
   droptable casdata='SmokingWeight' incaslib='mylib' quiet;
   droptable casdata='swDREstData' incaslib='mylib' quiet;
   droptable casdata='bartOutMod' incaslib='mylib' quiet;
run;

/* Création d'un ensemble de données plus grand simulé */
data mylib.SmokingWeight(promote=yes);
   array char_vars Sex $ Race $ Education $ Exercise $ Activity $;
   do i = 1 to 10000; /* Simule 10000 observations */
      call streaminit(12345 + i);
      do _n_ = 1 to dim(char_vars);
         char_vars(_n_) = ranword(char_vars(_n_)); /* Générer des valeurs de caractères aléatoires */
      end;

      if rand('UNIFORM') < 0.5 then Sex = 'Female'; else Sex = 'Male';
      if rand('UNIFORM') < 0.7 then Race = 'White'; else Race = 'Black';
      if rand('UNIFORM') < 0.6 then Education = 'HighSchool'; else Education = 'College';
      if rand('UNIFORM') < 0.5 then Exercise = 'Light'; else Exercise = 'Moderate';
      if rand('UNIFORM') < 0.5 then Activity = 'Active'; else Activity = 'Sedentary';
      Age = round(rand('NORMAL', 40, 10));
      YearsSmoke = round(rand('NORMAL', 15, 5));
      PerDay = round(rand('NORMAL', 15, 5));
      if rand('UNIFORM') < 0.5 then Quit = 0; else Quit = 1;
      Change = rand('NORMAL', 3, 1.5) + (Quit * rand('NORMAL', 2, 0.5)) - (PerDay * 0.1);
      output;
   end;
   drop i;
   label Quit = "Arrêt du Tabac (0=Non, 1=Oui)";
   label Change = "Changement de Poids (kg)";
run;

/* 1. Modélisation de la variable de traitement 'Quit' avec PROC LOGSELECT */
/* Utilisation de l'option NTHREADS pour le parallélisme en CAS */
proc logselect data=mylib.SmokingWeight nthreads=4;
   class Sex Race Education Exercise Activity;
   model Quit(Event='1') = Sex Race Education
         Exercise Activity Age YearsSmoke PerDay;
   output out=mylib.swDREstData pred=pTrt copyvars=(_ALL_);
run;

/* Calcul de la probabilité de la condition de contrôle (Quit=0) */
data mylib.swDREstData;
   set mylib.swDREstData;
   pCnt = 1 - pTrt;
run;

/* 2. Modélisation de la variable de résultat 'Change' avec PROC BART */
/* Utilisation d'un Ntree et Nmc appropriés pour de grands jeux de données */
proc bart data=mylib.swDREstData nTree=200 nMC=500 seed=2156;
   class Sex Race Education Exercise Activity Quit;
   model Change = Sex Race Education Exercise Quit Activity
         Age YearsSmoke PerDay;
   store out=mylib.bartOutMod;
run;

/* 3. Estimation des effets causaux avec PROC CAEFFECT (TMLE) pour grands ensembles de données */
proc caeffect data=mylib.swDREstData method=tmle;
   treatvar Quit;
   outcomevar Change;
   outcomemodel restore=mylib.bartOutMod predName=P_Change;
   pom treatLev=1 treatProb=pTrt;
   pom treatLev=0 treatProb=pCnt;
   difference evtLev=1;
run;