*================================================================*
* SAS programmer: James Nguyen, USEPA                            *
*                                                                *
* Project: OP inhibition studies                                 *
*                                                                *
* Chemical: 16 OPs                                               *
*                                                                *
* Descriptions:                                                  *
*     - perform regression diagnostics for final models of human *
*       vs. rat, adult vs. infant, male vs. female, and          *
*       Caucasian vs. non-Caucasian analyses                     *
*     - calculated ki values from human vs. rat analysis         *
*                                                                *
* Date: March 2020                                               *
*===============================================================*;

*===> human vs. rat;
title;
%macro DiagnosticsHR(chemical=, model=);

	*===> create scatter-plot of random effects Ai and Bi of each subject;

	ods select none;
	Proc SQL;
		create table ID as
		select distinct ID, Rat
		from Resid_HR
		order by ID;
	quit;
	Proc sort data=RandID_HR; by ID; run;

	data RandID_HR;
		length effect $4.;
		merge RandID_HR ID;
		by ID;
		effect=strip(upcase(effect));
		if index(id,'Rat')>0 and effect in ('UHA', 'UHB') then delete;
		if index(id,'Rat')=0 and effect in ('URA', 'URB') then delete;
		keep Rat ID effect estimate; 
	run;

	%if &model=2 %then %do;
		Proc means data = RandID_HR(where=(effect="UA")) mean std;
			var estimate;
			output out=A_parameters mean=A_mean std=A_std;
		run;

		Proc means data = RandID_HR(where=(effect="UB")) mean std;
			var estimate;
			output out=B_parameters mean=B_mean std=B_std;
		run;

		Proc SQL noprint;
			select A_mean, A_std
			into :A_loc, :A_scale
			from A_parameters;

			select B_mean, B_std
			into :B_loc, :B_scale
			from B_parameters;
		quit;

		Proc sort data = RandID_HR;
			by effect;
		run;

		Proc rank data = RandID_HR normal=blom out=RandID_HR;
			by effect;
			var estimate;
			ranks estimate_qq;
		run;

		Proc sort data = RandID_HR; by Rat; run;

		ods select default;

		title;
		title1 "&chemical: Normal Q-Q Plot of UA";
		Proc sgplot data = RandID_HR(where=(effect='UA')) ;
	 		styleattrs	datacolors        =(blue red)
						datacontrastcolors=(blue red);
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10) group=Rat;
			lineparm x=0 y=&A_loc slope=&A_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UB";
		Proc sgplot data = RandID_HR(where=(effect='UB')) ;
	 		styleattrs	datacolors        =(blue red)
						datacontrastcolors=(blue red);
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10) group=Rat;
			lineparm x=0 y=&B_loc slope=&B_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UB random effects";
		run;
		ods select none;
	%end;
	%if &model = 3 %then %do;
		Proc means data = RandID_HR(where=(effect="UHA")) mean std;
			var estimate;
			output out=HA_parameters mean=HA_mean std=HA_std;
		run;

		Proc means data = RandID_HR(where=(effect="UHB")) mean std;
			var estimate;
			output out=HB_parameters mean=HB_mean std=HB_std;
		run;

		Proc means data = RandID_HR(where=(effect="URA")) mean std;
			var estimate;
			output out=RA_parameters mean=RA_mean std=RA_std;
		run;

		Proc means data = RandID_HR(where=(effect="URB")) mean std;
			var estimate;
			output out=RB_parameters mean=RB_mean std=RB_std;
		run;

		Proc SQL noprint;
			select HA_mean, HA_std
			into :HA_loc, :HA_scale
			from HA_parameters;

			select HB_mean, HB_std
			into :HB_loc, :HB_scale
			from HB_parameters;

			select RA_mean, RA_std
			into :RA_loc, :RA_scale
			from RA_parameters;

			select RB_mean, RB_std
			into :RB_loc, :RB_scale
			from RB_parameters;
		quit;

		Proc sort data = RandID_HR; by effect; run;


		Proc rank data = RandID_HR normal=blom out=RandID_HR;
			by effect;
			var estimate;
			ranks estimate_qq;
		run;

		ods select default;
		title1 "&chemical: Normal Q-Q Plot of UHA";
		Proc sgplot data = RandID_HR(where=(effect="UHA")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=blue) group=Rat;
			lineparm x=0 y=&HA_loc slope=&HA_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_HR(where=(effect="UHB")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=blue) group=Rat;
			lineparm x=0 y=&HB_loc slope=&HB_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of URA";
		Proc sgplot data = RandID_HR(where=(effect="URA")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=red) group=Rat ;
			lineparm x=0 y=&RA_loc slope=&RA_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "URA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of URB";
		Proc sgplot data = RandID_HR(where=(effect="URB")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=red)  group=Rat;
			lineparm x=0 y=&RB_loc slope=&RB_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "URB random effects";
		run;
	%end;

	title;

*===> create plots of Residuals (=Observed Values - Subject Predicted Values);

	ods select none;
	data Resid_HR;
		set Resid_HR;
		resid=kapp-pred;
	run;
	proc sort data=Resid_HR;
		by Rat;
	run;

	Proc means data = Resid_HR mean std;
		var resid;
		output out=parameters mean=mean std=std;
	run;

	Proc SQL noprint;
		select mean, std
		into :loc, :scale
		from parameters;
	quit;

	Proc rank data = Resid_HR normal=blom out=Resid_HR;
		var resid;
		ranks resid_qq;
	run;

	Proc sort data = Resid_HR; by Rat; run;

	ods select default;

	title1 "&Chemical: Residuals vs. Predicted Values";
	proc sgplot data=Resid_HR;
	 	styleattrs	datacolors        =(blue red)
					datacontrastcolors=(blue red);
		scatter x=pred y=resid/markerattrs=(symbol=Circle size=10) group=Rat;
		yaxis label= "Residuals";
		refline 0/ axis=y;
	run;
	
	title "&Chemical: Normal Q-Q Plot of Residuals";
	Proc sgplot data = Resid_HR;
	 	styleattrs	datacolors        =(blue red)
					datacontrastcolors=(blue red);
		scatter x = resid_qq y = resid/markerattrs=(symbol=Circle size=10) group=Rat;
		lineparm x=0 y=&loc slope=&scale;
		xaxis label= "Normal Quantiles";
		yaxis label= "Residuals";
	run;
	title;
	Proc datasets nolist; save sasmacr &chemical: RandID_HR Parm; quit;

	option formdlim="=" ls=100 ps=100 nodate nonumber;
	ods noptitle;
	ods graphics off;
%Mend;*DiagnosticsHR;


*==> adult vs. infant;
title;
%macro DiagnosticsAI(chemical=, model=);
	ods select none;

	*===> create scatter-plot of random effects Ai and Bi of each subject;
	Proc SQL;
		create table ID as
		select distinct ID, infant
		from Resid_AI
		order by ID;
	quit;
	Proc sort data=RandID_AI; by ID; run;

	data RandID_AI;
		length effect $4.;
		merge RandID_AI ID;
		by ID;

		effect=strip(upcase(effect));

		if effect in ('UIA','UIB') and infant = 0 then delete;
		if effect in ('UADA','UADB') and infant = 1 then delete;

		keep ID infant effect estimate; 
	run;

	%if &model=2 %then %do;
		Proc means data = RandID_AI(where=(effect="UHA")) mean std;
			var estimate;
			output out=A_parameters mean=A_mean std=A_std;
		run;

		Proc means data = RandID_AI(where=(effect="UHB")) mean std;
			var estimate;
			output out=B_parameters mean=B_mean std=B_std;
		run;

		Proc SQL noprint;
			select A_mean, A_std
			into :A_loc, :A_scale
			from A_parameters;

			select B_mean, B_std
			into :B_loc, :B_scale
			from B_parameters;
		quit;

		Proc sort data = RandID_AI;
			by effect;
		run;

		Proc rank data = RandID_AI normal=blom out=RandID_AI;
			by effect;
			var estimate;
			ranks estimate_qq;
		run;

		Proc sort data = RandID_AI; by Infant; run;

		ods select default;

		title;
		title1 "&chemical: Normal Q-Q Plot of UHA";
		Proc sgplot data = RandID_AI(where=(effect='UHA')) ;
	 		styleattrs	datacolors        =(blue red)
						datacontrastcolors=(blue red);
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10) group=Infant;
			lineparm x=0 y=&A_loc slope=&A_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_AI(where=(effect='UHB')) ;
	 		styleattrs	datacolors        =(blue red)
						datacontrastcolors=(blue red);
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10) group=Infant;
			lineparm x=0 y=&B_loc slope=&B_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
		ods select none;
	%end;
	%if &model = 3 %then %do;
		Proc means data = RandID_AI(where=(effect="UADA")) mean std;
			var estimate;
			output out=AA_parameters mean=AA_mean std=AA_std;
		run;

		Proc means data = RandID_AI(where=(effect="UADB")) mean std;
			var estimate;
			output out=AB_parameters mean=AB_mean std=AB_std;
		run;

		Proc means data = RandID_AI(where=(effect="UIA")) mean std;
			var estimate;
			output out=IA_parameters mean=IA_mean std=IA_std;
		run;

		Proc means data = RandID_AI(where=(effect="UIB")) mean std;
			var estimate;
			output out=IB_parameters mean=IB_mean std=IB_std;
		run;

		Proc SQL noprint;
			select AA_mean, AA_std
			into :AA_loc, :AA_scale
			from AA_parameters;

			select AB_mean, AB_std
			into :AB_loc, :AB_scale
			from AB_parameters;

			select IA_mean, IA_std
			into :IA_loc, :IA_scale
			from IA_parameters;

			select IB_mean, IB_std
			into :IB_loc, :IB_scale
			from IB_parameters;
		quit;

		Proc sort data = RandID_AI; by effect; run;


		Proc rank data = RandID_AI normal=blom out=RandID_AI;
			by effect;
			var estimate;
			ranks estimate_qq;
		run;

		ods select default;
		title1 "&chemical: Normal Q-Q Plot of UHA";
		Proc sgplot data = RandID_AI(where=(effect="UADA")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=blue) group=Infant;
			lineparm x=0 y=&AA_loc slope=&AA_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_AI(where=(effect="UADB")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=blue) group=Infant;
			lineparm x=0 y=&AB_loc slope=&AB_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of HUA";
		Proc sgplot data = RandID_AI(where=(effect="UIA")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=red) group=Infant;
			lineparm x=0 y=&IA_loc slope=&IA_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_AI(where=(effect="UIB")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=red) group=Infant;
			lineparm x=0 y=&IB_loc slope=&IB_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
	%end;

	title;

*===> create plots of Residuals (=Observed Values - Subject Predicted Values);

	ods select none;
	data Resid_AI;
		set Resid_AI;
		resid=kapp-pred;
	run;
	proc sort data=Resid_AI;
		by Infant;
	run;

	Proc means data = Resid_AI mean std;
		var resid;
		output out=parameters mean=mean std=std;
	run;

	Proc SQL noprint;
		select mean, std
		into :loc, :scale
		from parameters;
	quit;

	Proc rank data = Resid_AI normal=blom out=Resid_AI;
		var resid;
		ranks resid_qq;
	run;

	Proc sort data = Resid_AI; by Infant; run;

	ods select default;

	title1 "&Chemical:  Residuals vs. Predicted Values";
	proc sgplot data=Resid_AI;
	 	styleattrs	datacolors        =(blue red)
					datacontrastcolors=(blue red);
		scatter x=pred y=resid/markerattrs=(symbol=Circle size=10) group=Infant;
		yaxis label= "Residuals";
		refline 0/ axis=y;
	run;
	
	title "&Chemical:  Normal Q-Q Plot of Residuals";
	Proc sgplot data = Resid_AI;
	 	styleattrs	datacolors        =(blue red)
					datacontrastcolors=(blue red);
		scatter x = resid_qq y = resid/markerattrs=(symbol=Circle size=10) group=Infant;
		lineparm x=0 y=&loc slope=&scale;
		xaxis label= "Normal Quantiles";
		yaxis label= "Residuals";
	run;
	title;
	Proc datasets nolist; save sasmacr &chemical:; quit;

	option formdlim="=" ls=100 ps=100 nodate nonumber;
	ods noptitle;
	ods graphics off;
%Mend;*DiagnosticsAI;



*===> male vs. female;
title;
%macro DiagnosticsMF(chemical=, model=);
	ods select none;

	*===> create scatter-plot of random effects Ai and Bi of each subject;
	Proc SQL;
		create table ID as
		select distinct ID, male
		from Resid_MF
		order by ID;
	quit;
	Proc sort data=RandID_MF; by ID; run;

	data RandID_MF;
		length effect $4.;
		merge RandID_MF ID;
		by ID;

		effect=strip(upcase(effect));

		if effect in ('UMA','UMB') and male = 0 then delete;
		if effect in ('UFA','UFB') and male = 1 then delete;

		keep ID male effect estimate; 
	run;

	%if &model=2 %then %do;
		Proc means data = RandID_MF(where=(effect="UHA")) mean std;
			var estimate;
			output out=A_parameters mean=A_mean std=A_std;
		run;

		Proc means data = RandID_MF(where=(effect="UHB")) mean std;
			var estimate;
			output out=B_parameters mean=B_mean std=B_std;
		run;

		Proc SQL noprint;
			select A_mean, A_std
			into :A_loc, :A_scale
			from A_parameters;

			select B_mean, B_std
			into :B_loc, :B_scale
			from B_parameters;
		quit;

		Proc sort data = RandID_MF;
			by effect;
		run;

		Proc rank data = RandID_MF normal=blom out=RandID_MF;
			by effect;
			var estimate;
			ranks estimate_qq;
		run;

		Proc sort data = RandID_MF; by Male; run;

		ods select default;

		title;
		title1 "&chemical: Normal Q-Q Plot of UHA";
		Proc sgplot data = RandID_MF(where=(effect='UHA')) ;
	 		styleattrs	datacolors        =(blue red)
						datacontrastcolors=(blue red);
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10) group=Male;
			lineparm x=0 y=&A_loc slope=&A_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_MF(where=(effect='UHB')) ;
	 		styleattrs	datacolors        =(blue red)
						datacontrastcolors=(blue red);
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10) group=Male;
			lineparm x=0 y=&B_loc slope=&B_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
		ods select none;
	%end;
	%if &model = 3 %then %do;
		Proc means data = RandID_MF(where=(effect="UMA")) mean std;
			var estimate;
			output out=MA_parameters mean=MA_mean std=MA_std;
		run;

		Proc means data = RandID_MF(where=(effect="UMB")) mean std;
			var estimate;
			output out=MB_parameters mean=MB_mean std=MB_std;
		run;

		Proc means data = RandID_MF(where=(effect="UFA")) mean std;
			var estimate;
			output out=FA_parameters mean=FA_mean std=FA_std;
		run;

		Proc means data = RandID_MF(where=(effect="UFB")) mean std;
			var estimate;
			output out=FB_parameters mean=FB_mean std=FB_std;
		run;

		Proc SQL noprint;
			select MA_mean, MA_std
			into :MA_loc, :MA_scale
			from MA_parameters;

			select MB_mean, MB_std
			into :MB_loc, :MB_scale
			from MB_parameters;

			select FA_mean, FA_std
			into :FA_loc, :FA_scale
			from FA_parameters;

			select FB_mean, FB_std
			into :FB_loc, :FB_scale
			from FB_parameters;
		quit;

		Proc sort data = RandID_MF; by effect; run;


		Proc rank data = RandID_MF normal=blom out=RandID_MF;
			by effect;
			var estimate;
			ranks estimate_qq;
		run;

		ods select default;
		title1 "&chemical: Normal Q-Q Plot of UHA";
		Proc sgplot data = RandID_MF(where=(effect="UFA")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=blue) group=Male;
			lineparm x=0 y=&FA_loc slope=&FA_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_MF(where=(effect="UFB")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=blue) group=Male;
			lineparm x=0 y=&FB_loc slope=&FB_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHA";
		Proc sgplot data = RandID_MF(where=(effect="UMA")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=red) group=Male;
			lineparm x=0 y=&MA_loc slope=&MA_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_MF(where=(effect="UMB")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=red) group=Male;
			lineparm x=0 y=&MB_loc slope=&MB_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
	%end;

	title;

*===> create plots of Residuals (=Observed Values - Subject Predicted Values);

	ods select none;
	data Resid_MF;
		set Resid_MF;
		resid=kapp-pred;
	run;
	proc sort data=Resid_MF;
		by Male;
	run;

	Proc means data = Resid_MF mean std;
		var resid;
		output out=parameters mean=mean std=std;
	run;

	Proc SQL noprint;
		select mean, std
		into :loc, :scale
		from parameters;
	quit;

	Proc rank data = Resid_MF normal=blom out=Resid_MF;
		var resid;
		ranks resid_qq;
	run;

	Proc sort data = Resid_MF; by Male; run;

	ods select default;

	title1 "&Chemical:  Residuals vs. Predicted Values";
	proc sgplot data=Resid_MF;
	 	styleattrs	datacolors        =(blue red)
					datacontrastcolors=(blue red);
		scatter x=pred y=resid/markerattrs=(symbol=Circle size=10) group=Male;
		yaxis label= "Residuals";
		refline 0/ axis=y;
	run;
	
	title "&Chemical:  Normal Q-Q Plot of Residuals";
	Proc sgplot data = Resid_MF;
	 	styleattrs	datacolors        =(blue red)
					datacontrastcolors=(blue red);
		scatter x = resid_qq y = resid/markerattrs=(symbol=Circle size=10) group=Male;
		lineparm x=0 y=&loc slope=&scale;
		xaxis label= "Normal Quantiles";
		yaxis label= "Residuals";
	run;
	title;
	Proc datasets nolist; save sasmacr &chemical:; quit;

	option formdlim="=" ls=100 ps=100 nodate nonumber;
	ods noptitle;
	ods graphics off;
%Mend;*DiagnosticsMF;





*===> Caucasian vs. Other;
title;
%macro DiagnosticsCO(chemical=, model=);
	ods select none;

	*===> create scatter-plot of random effects Ai and Bi of each subject;
	Proc SQL;
		create table ID as
		select distinct ID, Caucasian
		from Resid_CO
		order by ID;
	quit;
	Proc sort data=RandID_CO; by ID; run;

	data RandID_CO;
		length effect $4.;
		merge RandID_CO ID;
		by ID;

		effect=strip(upcase(effect));

		if effect in ('UCA','UCB') and Caucasian = 0 then delete;
		if effect in ('UORA','UORB') and Caucasian = 1 then delete;

		keep ID Caucasian effect estimate; 
	run;

	%if &model=2 %then %do;
		Proc means data = RandID_CO(where=(effect="UHA")) mean std;
			var estimate;
			output out=A_parameters mean=A_mean std=A_std;
		run;

		Proc means data = RandID_CO(where=(effect="UHB")) mean std;
			var estimate;
			output out=B_parameters mean=B_mean std=B_std;
		run;

		Proc SQL noprint;
			select A_mean, A_std
			into :A_loc, :A_scale
			from A_parameters;

			select B_mean, B_std
			into :B_loc, :B_scale
			from B_parameters;
		quit;

		Proc sort data = RandID_CO;
			by effect;
		run;

		Proc rank data = RandID_CO normal=blom out=RandID_CO;
			by effect;
			var estimate;
			ranks estimate_qq;
		run;

		Proc sort data = RandID_CO; by Caucasian; run;

		ods select default;

		title;
		title1 "&chemical: Normal Q-Q Plot of UHA";
		Proc sgplot data = RandID_CO(where=(effect='UHA')) ;
	 		styleattrs	datacolors        =(blue red)
						datacontrastcolors=(blue red);
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10) group=Caucasian;
			lineparm x=0 y=&A_loc slope=&A_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_CO(where=(effect='UHB')) ;
	 		styleattrs	datacolors        =(blue red)
						datacontrastcolors=(blue red);
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10) group=Caucasian;
			lineparm x=0 y=&B_loc slope=&B_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
		ods select none;
	%end;
	%if &model = 3 %then %do;
		Proc means data = RandID_CO(where=(effect="UCA")) mean std;
			var estimate;
			output out=CA_parameters mean=CA_mean std=CA_std;
		run;

		Proc means data = RandID_CO(where=(effect="UCB")) mean std;
			var estimate;
			output out=CB_parameters mean=CB_mean std=CB_std;
		run;

		Proc means data = RandID_CO(where=(effect="UORA")) mean std;
			var estimate;
			output out=OA_parameters mean=OA_mean std=OA_std;
		run;

		Proc means data = RandID_CO(where=(effect="UORB")) mean std;
			var estimate;
			output out=OB_parameters mean=OB_mean std=OB_std;
		run;

		Proc SQL noprint;
			select CA_mean, CA_std
			into :CA_loc, :CA_scale
			from CA_parameters;

			select CB_mean, CB_std
			into :CB_loc, :CB_scale
			from CB_parameters;

			select OA_mean, OA_std
			into :OA_loc, :OA_scale
			from OA_parameters;

			select OB_mean, OB_std
			into :OB_loc, :OB_scale
			from OB_parameters;
		quit;

		Proc sort data = RandID_CO; by effect; run;

		Proc rank data = RandID_CO normal=blom out=RandID_CO;
			by effect;
			var estimate;
			ranks estimate_qq;
		run;

		ods select default;
		title1 "&chemical: Normal Q-Q Plot of UHA";
		Proc sgplot data = RandID_CO(where=(effect="UORA")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=blue) group=Caucasian;
			lineparm x=0 y=&OA_loc slope=&OA_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_CO(where=(effect="UORB")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=blue) group=Caucasian;
			lineparm x=0 y=&OB_loc slope=&OB_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHA";
		Proc sgplot data = RandID_CO(where=(effect="UCA")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=red) group=Caucasian;
			lineparm x=0 y=&CA_loc slope=&CA_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHA random effects";
		run;
		title1 "&chemical: Normal Q-Q Plot of UHB";
		Proc sgplot data = RandID_CO(where=(effect="UCB")) ;
			scatter x = estimate_qq y = estimate/markerattrs=(symbol=CircleFilled size=10 color=red) group=Caucasian;
			lineparm x=0 y=&CB_loc slope=&CB_scale;
			xaxis label= "Normal Quantiles";
			yaxis label= "UHB random effects";
		run;
	%end;

	title;

*===> create plots of Residuals (=Observed Values - Subject Predicted Values);

	ods select none;
	data Resid_CO;
		set Resid_CO;
		resid=kapp-pred;
	run;
	proc sort data=Resid_CO;
		by Caucasian;
	run;

	Proc means data = Resid_CO mean std;
		var resid;
		output out=parameters mean=mean std=std;
	run;

	Proc SQL noprint;
		select mean, std
		into :loc, :scale
		from parameters;
	quit;

	Proc rank data = Resid_CO normal=blom out=Resid_CO;
		var resid;
		ranks resid_qq;
	run;

	Proc sort data = Resid_CO; by Caucasian; run;

	ods select default;

	title1 "&Chemical:  Residuals vs. Predicted Values";
	proc sgplot data=Resid_CO;
	 	styleattrs	datacolors        =(blue red)
					datacontrastcolors=(blue red);
		scatter x=pred y=resid/markerattrs=(symbol=Circle size=10) group=Caucasian;
		yaxis label= "Residuals";
		refline 0/ axis=y;
	run;
	
	title "&Chemical:  Normal Q-Q Plot of Residuals";
	Proc sgplot data = Resid_CO;
	 	styleattrs	datacolors        =(blue red)
					datacontrastcolors=(blue red);
		scatter x = resid_qq y = resid/markerattrs=(symbol=Circle size=10) group=Caucasian;
		lineparm x=0 y=&loc slope=&scale;
		xaxis label= "Normal Quantiles";
		yaxis label= "Residuals";
	run;
	title;
	Proc datasets nolist; save sasmacr &chemical:; quit;

	option formdlim="=" ls=100 ps=100 nodate nonumber;
	ods noptitle;
	ods graphics off;
%Mend;*DiagnosticsCO;



*===> create plots of individual sample curves;
%Macro scatter_plot(chemical=);

	Proc sort data = &chemical; by Human; run;
	title "&chemical: individual sample curves";
	Proc sgpanel data = &chemical;
		panelby Human;
		scatter x = _I__M y = kapp/group=ID;
		series x = _I__M y = kapp/group=ID;
		colaxis label = "Concentration";
	run;

	Proc sort data = &chemical; by Infant; run;
	Proc sgpanel data = &chemical(where=(Human=1));
		panelby Infant;
		scatter x = _I__M y = kapp/group=ID;
		series x = _I__M y = kapp/group=ID;
		colaxis label = "Concentration";
	run;

	Proc sort data = &chemical; by male; run;
	Proc sgpanel data = &chemical(where=(Human=1));
		panelby male;
		scatter x = _I__M y = kapp/group=ID;
		series x = _I__M y = kapp/group=ID;
		colaxis label = "Concentration";
	run;

	Proc sort data = &chemical; by Caucasian; run;
	Proc sgpanel data = &chemical(where=(Human=1));
		panelby Caucasian;
		scatter x = _I__M y = kapp/group=ID;
		series x = _I__M y = kapp/group=ID;
		colaxis label = "Concentration";
	run;
%Mend;



*===> calculate ki from NLMIXED models;
libname KI "C:\Users\JNguyen\Desktop\OP Inhibition\NLMIXED ki";
%Macro Calculate_Ki(chemical=, Hadj=, Radj=);
	data RandID_HR;
		length effect $4.;
		set RandID_HR;
		effect=strip(upcase(effect));

		*==> if model = 3;
		if index(id,'Rat')>0 and effect in ('UHA', 'UHB') then delete;
		if index(id,'Rat')=0 and effect in ('URA', 'URB') then delete;

		if index(id,'Rat')=0 and index(effect, 'A')>0 then effect = 'UHA';
		if index(id,'Rat')=0 and index(effect, 'B')>0 then effect = 'UHB';
		if index(id,'Rat')>0 and index(effect, 'A')>0 then effect = 'URA';
		if index(id,'Rat')>0 and index(effect, 'B')>0 then effect = 'URB';
		RandEffect=estimate;

		keep ID effect RandEffect;

	run;

	Data Parm;
		set Parm;
		if      strip(upcase(Parameter))='HA' then effect = 'UHA';
		else if strip(upcase(Parameter))='HB' then effect = 'UHB';
		else if strip(upcase(Parameter))='RA' then effect = 'URA';
		else if strip(upcase(Parameter))='RB' then effect = 'URB';
		else delete;
		FixEffect=estimate;
		keep effect FixEffect;
	run;

	Proc sort data = parm; by effect; run;
	Proc sort data = RandID_HR; by effect; run;

	Data Ind_parm;
		length effect $4.;
		merge parm RandID_HR;
		by effect;
		estimate= FixEffect + RandEffect;
	run;

	Proc sort data = Ind_parm;
		by ID effect;
	run;

	Proc transpose data = Ind_parm out = Ind_parm;
		by ID;
		ID effect;
		var estimate;
	run;

	Data Ind_parm;
		set Ind_parm;
		length Species $6.;
		if index(ID,'Rat')>0 then do; Species = "Rat"; ki=URA/(URB*&Radj); end;
		if index(ID,'Rat')=0 then do; Species = "Human"; ki=UHA/(UHB*&Hadj); end;
 		keep ID Species ki;
	run;

	Data KI.&chemical._ki;
		length Chemical $40.;
		set Ind_parm;
		Chemical = "&chemical";
	run;

	option formdlim="=" ls=100 ps=100 nodate nonumber;
	ods noptitle;
	ods graphics off;
%Mend;
title;


