*********************************************************;
* Date: 01/22/18                                        *;
* Programmer: Yunshen Jiao                              *;
*                                                       *;
* Purpose: This program was created for the analysis in *;
* Yunshen Jiao's thesis in partial fulfillment of the   *;
* requirements for the degree of Master of Public Health*; 
* in Epidemiology 2019. The original dataset is from    *;
* project Closed Loop Vagal Stimulation in Patients with*; 
* Posttraumatic Stress Disorder (did not publish yet in *;
* April, 2019)											*;
*********************************************************;


/*Data Cleaning did not included in this file*/

/*Recode confounder variables: race, education level, biological sex for the following analysis*/
/*Create new variable: day of blood draw for the following analysis, based on the procedure 
described in Figure 1*/

data combine;
	set darpa.biomarker;
/*Race*/
	if race_id___1 = 1 then race = 1;
	if race_id___2 = 1 then race = 2;
	if race_id___4 = 1 then race = 3;
	if race_id___5 = 1 then race = 3;
	if race_id___6 = 1 then race = 3;
	if race_id___7 = 1 then race = 3;
	if race_id___8 = 1 then race = 3;
/*Education Level*/
	if education in (2,3,4) then educ = 0;
	if education in (5,6) then educ = 1;
/*Biological Sex*/
	if biol_sex = 1 then gender = 0;
	else if biol_sex = 2 then gender = 1;
/*New variable: Blood draw day*/
	if timepoints = 0 then timeday = 0;
	else if timepoints <= 14 then timeday = 1;
	else if timepoints <= 18 then timeday = 2;
	else if timepoints <= 22 then timeday = 3;
run;

/*Baseline Level test*/

/*Model 1-4: Association between PACAP blood levels and PTSD diagnosis*/

/*Model 1*/
proc reg data = combine;
	model logPACAP = PTSDorNOT;
run;
/*Model 2*/
proc glm data = combine;
	model logPACAP = PTSDorNOT gender PTSDorNOT*gender;
run;
/*Model 3*/
proc glm data = combine;
	model logPACAP = PTSDorNOT gender age BMI race educ;
run;
/*Model 4*/
proc glm data = combine;
	model logPACAP = PTSDorNOT gender age BMI race educ PTSDorNOT*gender;
run;

/*Model 5: Association between PACAP blood levels and psychological scores (linear)*/
/*Create series of variables for psychological scores*/
proc sql noprint;
   		select NAME
   		into :psyscore_1 - :psyscore_30
		from vars
		where varnum in (77,82,84,87,91,96,99,104,108,112,113,115,116,118,119,121,122,124,125,127,128,130,131,135,162,163,175,178,182,197)
		order by varnum; 
quit;
%macro psybio;
	data test;
		set _null_;
	run;
%do j = 1 %to 23;
	data &&ps_&j;
		set _null_;
	run;
	proc glm data = baseline;
		ods output ParameterEstimates = &&ps_&j ;
		class gender(ref = "0") educ (ref="0") race (ref = "1");
		model &&ps_&j = logPACAP age gender bmi educ race/ solution CLPARM;
	run;
	ods output close;
	data test;
		set test &&ps_&j;
	run;
	%end;
%mend;
%psybio;

/*Model 6: Association between PACAP blood levels and psychological scores (dichotomous)*/
/*Rank all the scores needed*/
proc rank data = rank out = rank groups = 2;
var PSS_score HAMA_score HAMD_score ETI_ALL ETI_ALL_2 BAECKE_work BAECKE_SportIndex ESSI_Score anger_state CADSS_total;
ranks PSS_2c HAMA_2c HAMD_2c ETI_rank ETI_rank_2 BAECKE_work_rank BAECKE_sport_rank ESSI_Score_rank anger_state_rank CADSS_rank;
run;
/*Models*/
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model PSS_2c(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model CADSS_rank(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model anger_state_rank(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model ESSI_Score_rank(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model BAECKE_sport_rank(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model BAECKE_work_rank(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model BDI_2c(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model ETI_2c(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model HAMA_2c(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model HAMD_2c(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model ETI_rank(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;
proc logistic data=baseline_rank1;
   class biol_sex race (ref = '1') educ(ref = '1');
   model ETI_rank_2(event = '1') = logPACAP age biol_sex race educ bmi/cl expb;
run;

/*Model 7: Mixed model to examine the assocation between treatment or PTSD diagnosis with PACAP concentration over time*/
proc glimmix data = combine method = rspl PLOTS=pearsonpanel;
	class caps_index timeday(ref = '0') race(ref = '1') educ (ref = '1');
	model logPACAP = PTSDorNOT treatment age gender BMI race educ timeday PTSDorNOT*timeday treatment*timeday/ s cl;
	random intercept / subject = caps_index type = un g v vcorr;
	covtest 'Chunk test for random effects' zerog;
estimate 'Day 1 change in PACAP concentration among healthy sham patients ' timeday 1 0 0 -1 / cl exp; 
estimate 'Day 2 change in PACAP concentration among healthy sham patients ' timeday 0 1 0 -1 / cl exp; 
estimate 'Day 3 change in PACAP concentration among healthy sham patients ' timeday 0 0 1 -1 / cl exp; 
estimate 'Day 1 change in PACAP concentration among PTSD sham patients ' timeday 1 0 0 -1 PTSDorNOT*timeday 1 0 0 -1/ cl exp; 
estimate 'Day 2 change in PACAP concentration among PTSD sham patients ' timeday 0 1 0 -1 PTSDorNOT*timeday 0 1 0 -1/ cl exp; 
estimate 'Day 3 change in PACAP concentration among PTSD sham patients ' timeday 0 0 1 -1 PTSDorNOT*timeday 0 0 1 -1/ cl exp; 
estimate 'Day 1 change in PACAP concentration among healthy VNS patients ' timeday 1 0 0 -1 treatment*timeday 1 0 0 -1/ cl exp; 
estimate 'Day 2 change in PACAP concentration among healthy VNS patients ' timeday 0 1 0 -1 treatment*timeday 0 1 0 -1/ cl exp; 
estimate 'Day 3 change in PACAP concentration among healthy VNS patients ' timeday 0 0 1 -1 treatment*timeday 0 0 1 -1/ cl exp; 
estimate 'Day 1 change in PACAP concentration among PTSD VNS patients ' timeday 1 0 0 -1 PTSDorNOT*timeday 1 0 0 -1 treatment*timeday 1 0 0 -1/ cl exp; 
estimate 'Day 2 change in PACAP concentration among PTSD VNS patients ' timeday 0 1 0 -1 PTSDorNOT*timeday 0 1 0 -1 treatment*timeday 0 1 0 -1/ cl exp; 
estimate 'Day 3 change in PACAP concentration among PTSD VNS patients ' timeday 0 0 1 -1 PTSDorNOT*timeday 0 0 1 -1 treatment*timeday 0 0 1 -1/ cl exp; 
run;

/*Model 8: Mixed model to examine the assocation between gender or PTSD diagnosis with PACAP concentration over time*/
proc glimmix data = combine method = rspl PLOTS=pearsonpanel;
	class caps_index timeday(ref = '0') race(ref = '1') educ (ref = '1');
	model logPACAP = PTSDorNOT gender age BMI race educ timeday PTSDorNOT*timeday gender*timeday/ s cl;
	random intercept / subject = caps_index type = un g v vcorr;
	covtest 'Chunk test for random effects' zerog;
estimate 'Day 1 change in PACAP concentration among healthy male patients ' timeday 1 0 0 -1 / cl exp; 
estimate 'Day 2 change in PACAP concentration among healthy male patients ' timeday 0 1 0 -1 / cl exp; 
estimate 'Day 3 change in PACAP concentration among healthy male patients ' timeday 0 0 1 -1 / cl exp; 
estimate 'Day 1 change in PACAP concentration among PTSD male patients ' timeday 1 0 0 -1 PTSDorNOT*timeday 1 0 0 -1/ cl exp; 
estimate 'Day 2 change in PACAP concentration among PTSD male patients ' timeday 0 1 0 -1 PTSDorNOT*timeday 0 1 0 -1/ cl exp; 
estimate 'Day 3 change in PACAP concentration among PTSD male patients ' timeday 0 0 1 -1 PTSDorNOT*timeday 0 0 1 -1/ cl exp; 
estimate 'Day 1 change in PACAP concentration among healthy female patients ' timeday 1 0 0 -1 gender*timeday 1 0 0 -1/ cl exp; 
estimate 'Day 2 change in PACAP concentration among healthy female patients ' timeday 0 1 0 -1 gender*timeday 0 1 0 -1/ cl exp; 
estimate 'Day 3 change in PACAP concentration among healthy female patients ' timeday 0 0 1 -1 gender*timeday 0 0 1 -1/ cl exp; 
estimate 'Day 1 change in PACAP concentration among PTSD female patients ' timeday 1 0 0 -1 PTSDorNOT*timeday 1 0 0 -1 gender*timeday 1 0 0 -1/ cl exp; 
estimate 'Day 2 change in PACAP concentration among PTSD female patients ' timeday 0 1 0 -1 PTSDorNOT*timeday 0 1 0 -1 gender*timeday 0 1 0 -1/ cl exp; 
estimate 'Day 3 change in PACAP concentration among PTSD female patients ' timeday 0 0 1 -1 PTSDorNOT*timeday 0 0 1 -1 gender*timeday 0 0 1 -1/ cl exp; 
run;

/*Model 9: Mixed model to examine the assocation between gender or treatment effect with PACAP concentration over time*/
proc glimmix data = combine method = rspl PLOTS=pearsonpanel;
	class caps_index timeday(ref = '0') race(ref = '1') educ (ref = '1');
	model logPACAP = treatment gender age BMI race educ timeday treatment*timeday gender*timeday/ s cl;
	random intercept / subject = caps_index type = un g v vcorr;
	covtest 'Chunk test for random effects' zerog;
estimate 'Day 1 change in PACAP concentration among male patients received Sham Stimulation' timeday 1 0 0 -1 / cl exp; 
estimate 'Day 2 change in PACAP concentration among male patients received Sham Stimulation ' timeday 0 1 0 -1 / cl exp; 
estimate 'Day 3 change in PACAP concentration among male patients received Sham Stimulation ' timeday 0 0 1 -1 / cl exp; 
estimate 'Day 1 change in PACAP concentration among male patients received VNS' timeday 1 0 0 -1 treatment*timeday 1 0 0 -1/ cl exp; 
estimate 'Day 2 change in PACAP concentration among male patients received VNS' timeday 0 1 0 -1 treatment*timeday 0 1 0 -1/ cl exp; 
estimate 'Day 3 change in PACAP concentration among male patients received VNS' timeday 0 0 1 -1 treatment*timeday 0 0 1 -1/ cl exp; 
estimate 'Day 1 change in PACAP concentration among female patients received Sham Stimulation' timeday 1 0 0 -1 gender*timeday 1 0 0 -1/ cl exp; 
estimate 'Day 2 change in PACAP concentration among female patients received Sham Stimulation' timeday 0 1 0 -1 gender*timeday 0 1 0 -1/ cl exp; 
estimate 'Day 3 change in PACAP concentration among female patients received Sham Stimulation' timeday 0 0 1 -1 gender*timeday 0 0 1 -1/ cl exp; 
estimate 'Day 1 change in PACAP concentration among female patients received VNS' timeday 1 0 0 -1 treatment*timeday 1 0 0 -1 gender*timeday 1 0 0 -1/ cl exp; 
estimate 'Day 2 change in PACAP concentration among female patients received VNS' timeday 0 1 0 -1 treatment*timeday 0 1 0 -1 gender*timeday 0 1 0 -1/ cl exp; 
estimate 'Day 3 change in PACAP concentration among female patients received VNS' timeday 0 0 1 -1 treatment*timeday 0 0 1 -1 gender*timeday 0 0 1 -1/ cl exp; 
run;

/*Model 10: Association between PACAP and the other cytokines' concentration*/
/*Create series of variables for log-transformed biomarker concentration*/
proc sql noprint;
   		select NAME
   		into :vars_1 - :vars_13
		from vars
		where 53 <= varnum <= 65
		order by varnum; 
quit;
/*note: vars_13 is the log transformed PACAP concentration*/
/*Models*/
%macro linear;
	data test;
		set _null_;
	run;
%do i = 1 %to 12;
	data &&vars_&i;
		set _null_;
	run;
	proc glm data = combine;
		ods output ParameterEstimates = &&vars_&i (where=(Parameter not eq"Intercept"));
		class biol_sex educ(ref = '1') race(ref="1");
		model logPACAP = &&vars_&i age biol_sex bmi educ race/ solution;
	run;
	ods output close;
	data test;
		set test &&vars_&i;
	run;
%end;
%mend;
%linear;