*====================================================================*
* SAS programmer: James Nguyen, USEPA                                *
*                                                                    *
* Project: OP inhibition studies                                     *
*                                                                    *
* Chemical: 16 OPs                                                   *
*                                                                    *
* Descriptions:                                                      *
*   - calculate the coverage ratios of ki distributions (i.e., the   *
*     90%-tile/Geometric Mean, 95%-tile/GM, 97.5%-tile/GM, etc.) for *
*     each chemical.  Note that the ki values were obtained from     *
*     NLMIXED models in the human vs. rat analyses.                  *
*                                                                    *
*   - create Q-Q normal quantile plots for ki values                 *
*                                                                    *
* Date: March 2020                                                   *
*===================================================================*;
option FormDlim="=" nodate nonumber ls=100 ps=100;

libname KI "C:\Users\JNguyen\Desktop\OP Inhibition\NLMIXED ki";

Data All_KI;
	set KI.bensulideoxon_ki
		KI.chlorethoxyfos_ki
		KI.ddvp_ki
		KI.dicrotophos_ki
		KI.ethoprop_ki
		KI.fenamiphos_ki
		KI.malaoxon_ki
		KI.methamidophos_ki
		KI.naled_ki
		KI.omethoate_ki
		KI.phorateoxonsulfone_ki
		KI.phorateoxonsulfoxide_ki
		KI.phosmetoxon_ki
		KI.tebuprimphosoxon_ki
		KI.terbufosoxonsulfone_ki
		KI.terbufosoxonsulfoxide_ki;
	lnki=log(ki);
run;


*===> calculate the coverage ratios of ki distributions;

Proc SQL;
	create table Coverage_ratios as
		select chemical, species, MEAN_ki, SD_ki, 
			(MEAN_ki+probit(0.975)*SD_ki)/(MEAN_ki+probit(0.025)*SD_ki) as P975_p025, 
			(MEAN_ki+probit(0.950)*SD_ki)/(MEAN_ki+probit(0.050)*SD_ki) as P950_p050,
			(MEAN_ki+probit(0.900)*SD_ki)/(MEAN_ki) as P900_mean,
			(MEAN_ki+probit(0.950)*SD_ki)/(MEAN_ki) as P950_mean,
			(MEAN_ki+probit(0.975)*SD_ki)/(MEAN_ki) as P975_mean,
			(MEAN_ki+probit(0.990)*SD_ki)/(MEAN_ki) as P990_mean,

			exp(MEAN_lnki) as GM, exp(SD_lnki) as GSD,
			exp(MEAN_lnki+probit(0.975)*SD_lnki)/exp(MEAN_lnki+probit(0.025)*SD_lnki) as P975_p025l, 
			exp(MEAN_lnki+probit(0.950)*SD_lnki)/exp(MEAN_lnki+probit(0.050)*SD_lnki) as P950_p050l,
			exp(MEAN_lnki+probit(0.900)*SD_lnki)/exp(MEAN_lnki) as P900_meanl,
			exp(MEAN_lnki+probit(0.950)*SD_lnki)/exp(MEAN_lnki) as P950_meanl,
			exp(MEAN_lnki+probit(0.975)*SD_lnki)/exp(MEAN_lnki) as P975_meanl,
			exp(MEAN_lnki+probit(0.990)*SD_lnki)/exp(MEAN_lnki) as P990_meanl

		From 	(select chemical, species, avg(ki) as MEAN_ki, std(ki) as SD_ki, avg(lnki) as MEAN_lnki, std(lnki) as SD_lnki
				 from All_KI
				 group by chemical, species)
		order by species, chemical;
quit;

Proc export data = Coverage_ratios dbms=xlsx
	outfile="C:\Users\JNguyen\OneDrive - Environmental Protection Agency (EPA)\Animal studies\50773501-2-3 In Vitro Inhibition Study\from James\Summarized results ratios.xlsx";
run;

ods noptitle;

*===> create Q-Q plots for ki values;

%Macro QQPLOT(chemical);
	ods select none;

	Data temp;
		set All_KI(where=(chemical="&chemical" and Species="Human"));
	run;

	Proc means data = temp mean std;
		var ki;
		output out=MeanSD mean=mean std=SD;
	run;

	Proc means data = All_KI(where=(chemical="&chemical" and Species="Human")) mean std;
		var lnki;
		output out=MeanSDln mean=meanln std=SDln;
	run;

	Proc SQL noprint;
		select mean, SD
		into :loc, :scale
		from MeanSD;

		select meanln, SDln
		into :locln, :scaleln
		from MeanSDln;
	quit;

	Proc rank data = temp normal=blom out=temp;
		var ki;
		ranks ki_qq;
	run;
	Proc rank data = temp normal=blom out=temp;
		var lnki;
		ranks lnki_qq;
	run;

	ods select default;
	title;
	title1 "&chemical:  Human ki";
	Proc sgplot data = temp noautolegend;
	 	styleattrs	datacolors        =(blue)
					datacontrastcolors=(blue);
		scatter x = ki_qq y = ki/markerattrs=(symbol=CircleFilled size=10);
		lineparm x=0 y=&loc slope=&scale;
		xaxis label= "Normal Quantiles";
		yaxis label= "Human ki";
	run;
	title1 "&chemical:  Human ln(ki)";
	Proc sgplot data = temp noautolegend;
	 	styleattrs	datacolors        =(blue)
					datacontrastcolors=(blue);
		scatter x = lnki_qq y = lnki/markerattrs=(symbol=CircleFilled size=10);
		lineparm x=0 y=&locln slope=&scaleln;
		xaxis label= "Normal Quantiles";
		yaxis label= "Human ln(ki)";
	run;
%Mend;
ods rtf file = "C:\Users\JNguyen\Desktop\OP Inhibition\NLMIXED ki\Q-Q Plot of ki.rtf" startpage=no;
%QQPLOT(BensulideOxon);
%QQPLOT(Chlorethoxyfos);
%QQPLOT(DDVP);
%QQPLOT(Dicrotophos);
%QQPLOT(Ethoprop);
%QQPLOT(Fenamiphos);
%QQPLOT(Malaoxon);
%QQPLOT(Methamidophos);
%QQPLOT(Naled);
%QQPLOT(Omethoate);
%QQPLOT(PhorateOxonSulfone);
%QQPLOT(PhorateOxonSulfoxide);
%QQPLOT(PhosmetOxon);
%QQPLOT(TebuprimphosOxon);
%QQPLOT(TerbufosOxonSulfone);
%QQPLOT(TerbufosOxonSulfoxide);
ods rtf close;
