*using file "zambia pilot data for stats_resubmission";

*data prep;
proc format; 
	value groupf 	1='Orphans'
					0='Control'; 
	value sexf  1='F'
				0='M';
	value seasonf 0='Early Dry'
	              1='Late Dry';
run;

data pilot;
	set pilot;
	if not missing(ID); *Drop the blank row of data;

	format group groupf. sex sexf. season seasonf.;
run;

Data pilot;
	Set pilot;
	logCort = log10(cort);
	logThyroid =log10(thyroid);

	age2 = age*age; 
	label group = 'Group';
run;
*Descriptive tables;
*Overall;
Proc means data=pilot n mean std min max median maxdec=2;
var age thyroid cort BCS;
Run;
*by group & season;
 proc sort data=pilot;
by group season; 
run;
*By Group;
proc sort data=pilot;
by group; 
run;
Proc means data=pilot n mean std min max median maxdec=2;
By group;
var age thyroid cort BCS;
Run;
*By Group & Season;
proc sort data=pilot;
by group season; 
run;
Proc means data=pilot n mean std min max median maxdec=2;
By group season;
var age thyroid cort BCS;
Run;

* Get values of low and high age to use in slicing interaction;
proc means data=pilot min max median mean stddev maxdec=2;
	var age;
run;
*Models;

*** fGCM - Cortisol ***;
proc mixed data=pilot;
	class season id group;
	model logcort = group age season age*group group*season/ solution cl;
	Repeated season/subject= ID;
	lsmeans group / at age=4 pdiff=all cl; *min;
	lsmeans group / at age=46.2 pdiff=all cl; *1SD below mean;
	lsmeans group / pdiff=all cl; *at mean;
	lsmeans group / at age=137.4 pdiff=all cl; *1SD above mean;
	lsmeans group / at age=186 pdiff=all cl; *max;
	lsmeans group|season / pdiff = all cl;
run;
** PLOT **;
*For the sake of making a PLOT don't need season or group*season or random effects;
ods exclude all;  *Suppress output;
proc mixed data=pilot;
	class id group;
	model logcort = group age age*group / 
		solution outpred=predCort cl;
run;
ods select all; * Turn output on again;
proc sgplot data=predCort;
	styleattrs
   	datacontrastcolors= (teal orange);
    scatter x=age y=logcort / group = group markerattrs=(symbol=CircleFilled);
	series x=age y=Pred / group=group lineattrs=(thickness=2);
    xaxis label="Age (months)";
    yaxis label="fGCM (Log10)";
    title "Predicted Values of fGCM (Log10) by Age and Group";
run;
***********************************************************
*** fT3 - Thyroid ***;
proc mixed data=pilot;
	class season id group;
	model logthyroid = group age season age*group age2 age2*group group*season / 
		solution cl; 
	Repeated season/subject= ID;
		lsmeans group / at (age age2)=(4 16) pdiff=all cl; *Min;	
		lsmeans group / at (age age2)=(46.2 2134) pdiff=all cl; *1 SD below;
		lsmeans group / pdiff=all cl; * Mean;
		lsmeans group / at (age age2)=(137.4 18879) pdiff=all cl; * 1 SD above;
		lsmeans group / at (age age2)=(186 34596) pdiff=all cl; *Max;	
		lsmeans group|season / pdiff = all cl;
run;
*For the sake of making a PLOT we don't need season or random effects;
ods exclude all; *Suppress output from this model;
proc sort data=pilot;
	by descending Group Season ;
run;
proc mixed data=pilot;
	class season id group;
	model logthyroid = group age age*group age2 age2*group/ 
		solution outpred=predThyroid cl;
run;
ods select all; *Turn output on again;
proc sort data=predThyroid; by age Season; run;
proc sgplot data=predThyroid;
	styleattrs
   	datacontrastcolors= (teal orange);
    scatter x=age y=logthyroid / group = group markerattrs=(symbol=CircleFilled);
	series x=age y=Pred / group=group lineattrs=(thickness=2);
    xaxis label="Age (months)";
    yaxis label="fT3 (Log10)";
    title "Predicted Values of fT3 (Log 10) by Age and Group";
run;
** BCS **;
proc sort data=pilot; by id season; run;
proc mixed data=pilot;
	class season id group;
	model BCS = group age season group*season/ solution outpred=PredBCS cl;
	Repeated season/subject= ID;
	lsmeans group|season / at age=91.66 pdiff = all cl;
	lsmeans group / at age=0 pdiff;
		lsmeans group / at age=1 pdiff;
		lsmeans group / at age=50 pdiff;
run;
* For the PLOT *;
ods exclude;
proc mixed data=pilot;
	class season id group;
	model BCS = group age / solution outpred=PredBCS cl;
run;
ods select all;

proc sort data=PredBCS; by age; run;
proc sgplot data=PredBCS;
styleattrs
   	datacontrastcolors= (teal orange);
    scatter x=age y=BCS / group=group markerattrs=(symbol=CircleFilled);
    series x=age y=Pred / group=group lineattrs=(thickness=2);
    xaxis label="Age (months)";
    yaxis label="BCS";
    title "Predicted Values of BCS by Age and Group";
run;
proc sgplot data=PredBCS;
	reg x=age y=BCS / group = group clm;
run;

*to see if sex is signifcant using only individuals where sex is known;
proc mixed data=pilot;
	where sex in (0,1);
	class season id group;
	model logcort = sex/ solution;
	Repeated season/subject= ID;
run;
*p=0.408;

proc mixed data=pilot;
	where sex in (0,1);
	class season id group sex;
	model logthyroid = sex/ solution;
	Repeated season/subject= ID;
	lsmeans sex / diff=all cl;
run;
*p=0.0678;

proc mixed data=pilot;
	where sex in (0,1);
	class season id group;
	model BCS = sex/ solution;
	Repeated season/subject= ID;
run;
*p=0.122;
*Boxplots;
ods graphics on;
data myfill;
	length id $4 value FillColor $7;
	input id value FillColor;
	datalines;
	MyID Orphans orange
	MyID Control teal
;
proc sort data=pilot;
	by = group;
run;
proc sgplot data=pilot dattrmap=myfill;
styleattrs datacontrastcolors=(black) datasymbols = (circle);
vbox logcort/ category=season group=group attrid=MyID;
yaxis label = 'log fGCM';
title 'Boxplot of Log fGCM for Group by Season';
 run;
 proc sgplot data=pilot dattrmap=myfill;
 styleattrs
   	datacontrastcolors= (black)datasymbols = (circle);
vbox logthyroid/ category=season group=group attrid=MyID;
yaxis label="log fT3";
title 'Boxplot of Log fT3 for Group by Season';
 run;