*================================================================*
* SAS programmer: James Nguyen, USEPA                            *
*                                                                *
* Project: OP inhibition studies                                 *
*                                                                *
* Chemical: DDVP                                                 *
*                                                                *
* Descriptions:                                                  *
*     - run the selected (Exponent) final models                 *
*     - perform regression diagnostics                           *
*     - calculated ki values from human vs. rat analysis         *
*                                                                *
* Date: March 2020                                               *
*===============================================================*;
option formdlim="=" ls=120 ps=120 nodate nonumber;

%include "C:\Users\JNguyen\OneDrive - Environmental Protection Agency (EPA)\Animal studies\50773501-2-3 In Vitro Inhibition Study\from James\SAS code Final Models\Diagnostics.SAS";
libname EPA "C:\Users\JNguyen\OneDrive - Environmental Protection Agency (EPA)\Animal studies\50773501-2-3 In Vitro Inhibition Study\from James\kapp data";
/****************************************************************************/
/* Final NLMIXED Model                                                      */
/****************************************************************************/
data DDVP;
     set epa.DDVP;
run;
proc sort data=DDVP;
     by id;
run;

ods rtf file="C:\Users\JNguyen\Desktop\OP Inhibition\Diagnostics\DDVP - individual sample curves.rtf" startpage=no;
%scatter_plot(chemical=DDVP);
ods rtf close;


*===> Human vs. Rat;

title "model 3: a model with different random effects between species";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=DDVP /*itdetails*/ tech=nrridg lis=3 gconv=0;
      PARMS 	ha=0.7261	HB=1.4717 	ra=1.2460  RB=2.8938

			GHA2=0.3503 ghab=0.9259 GHB2=0.2999
			GRA2=1 to 10 GRAB2=0.1 to 2 by 0.1 GRB2=1 to 10		S2E=0.3645/best=1000 ;

  	 bounds ha>0, HB>0, ra>0, RB>0, GHA2>0, GHB2>0, GRA2 > 0, GRB2 > 0, S2E>0;

          	rha=ha+uha;
            rhb=HB*10+uhb;
			rra=ra+ura;
			rrb=RB*10+urb;
            PRED = (rha*H_I__M + rra*R_I__M)/(rhb*Human*1e-6 + H_I__M + rrb*Rat*1e-6 + R_I__M);

      MODEL KAPP ~ NORMAL(PRED, S2E*0.001);

      RANDOM  uha uhb ura urb ~ normal([0, 0, 0, 0],   [GHA2/10, 
														ghab, GHB2*100, 
														0,    0,     GRA2/10, 
														0,    0,     GRAB2*10,   GRB2*100]) SUBJECT=ID out=RandID_HR; 
	  ESTIMATE "human/rat" (ha/hb)/(ra/rb);
	  ESTIMATE 'human ki' (ha/hb);
      ESTIMATE 'rat ki' (ra/rb);
	  predict pred out=Resid_HR;
run;

Data Parm;
	set Parm;
	if strip(upcase(Parameter)) in ("HB","RB") then estimate = estimate*10;
run;

ods rtf file="C:\Users\JNguyen\Desktop\OP Inhibition\Diagnostics\DDVP human-rat.rtf" startpage=no;
%DiagnosticsHR(chemical=DDVP, model=3)
ods rtf close;
%Calculate_Ki(chemical=DDVP, Hadj=1e-6, Radj=1e-6);


*===> Adult vs. Infant;

title "final model 2b - Adult (>=10 yr) vs. Infant";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=DDVP(where=(Human=1)) tech=nrridg lis=3 maxiter=1000 method=firo;
      PARMS ada=.7028	adB=14.6871
			ia=.7978 iB=14.3560
			gha=.3100 	ghab=.8459 	ghb=.2717
			S2E=0.3657;

			rada=ada+uha;
			ria=ia+uha;
            radb=adB+uhb;
			rib=iB+uhb;
            PRED = (rada*H_I__M*adult + ria*H_I__M*infant)/(radb*adult*1e-6 + H_I__M*adult + rib*1e-6*infant + H_I__M*infant);
      MODEL KAPP ~ NORMAL(PRED, S2E*1e-3);
      RANDOM  uha uhb ~ normal([0, 0],  [gha*1e-1, ghab, ghb*100]) SUBJECT=ID out=RandID_AI; 
	  ESTIMATE "Adult/Infant ki" (ada/adb)/(ia/ib);
	  ESTIMATE 'Adult ki' (ada/adb);
      ESTIMATE 'Infant ki' (ia/ib);
	  predict pred out=Resid_AI;
run;


ods rtf file="C:\Users\JNguyen\Desktop\OP Inhibition\Diagnostics\DDVP Adult-Infant.rtf" startpage=no;
%DiagnosticsAI(chemical=DDVP, model=2)
ods rtf close;

*===> Male vs. Female;

title "model 2b - male vs. female";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=DDVP(where=(Human=1)) method=firo tech=newrap lis=3 maxiter=100000;
      PARMS ma=.665 mb= .13
			fa=.77 fB=.16
			gha=.031 	ghab=.0082 	ghb=.00267
			S2E=.3617;

            rma=ma+uha;
			rfa=fa+uha;
            rmb=mB+uhb;
			rfb=fB+uhb;
 			PRED = (rma*H_I__M*male + rfa*H_I__M*female)/(rmb*male*1e-4 + H_I__M*male + rfb*female*1e-4 + H_I__M*female);
      MODEL KAPP ~ NORMAL(PRED, S2E*1e-3);
      RANDOM  uha uhb ~ normal([0, 0],  [gha, ghab, ghb]) SUBJECT=ID out=RandID_MF; 
	  ESTIMATE "male/female ki" (ma/mb)/(fa/fb);
	  ESTIMATE 'male ki' (ma/mb);
      ESTIMATE 'female ki' (fa/fb);
	  predict pred out=Resid_MF;
run;


ods rtf file="C:\Users\JNguyen\Desktop\OP Inhibition\Diagnostics\DDVP Male-Female.rtf" startpage=no;
%DiagnosticsMF(chemical=DDVP, model=2)
ods rtf close;

*===> Caucasian vs. Other Races;

title "model 2 - Caucasian vs. Other Race";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=DDVP(where=(Human=1)) method=firo tech=newrap lis=3 maxiter=100000;
      PARMS ca=.7001	cB=13.5
			ora=.7888 orB=19.43
			gha=.3262	ghab=.8162	ghb=.2319
			S2E=.3585;

            rca=ca+uha;
			rora=ora+uha;
            rcb=cB+uhb;
			rorb=orB+uhb;
		PRED = (rca*H_I__M*caucasian + rora*H_I__M*otherrace)/(rcb*caucasian*1e-6 + H_I__M*caucasian + rorb*otherrace*1e-6 + H_I__M*otherrace);
      MODEL KAPP ~ NORMAL(PRED, S2E*1e-3);
      RANDOM  uha uhb ~ normal([0, 0],  [gha*1e-1, ghab, ghb*100]) SUBJECT=ID out=RandID_CO; 
	  ESTIMATE "Caucasian/Other Race ki" (ca/cb)/(ora/orb);
	  ESTIMATE 'Caucasian ki' (ca/cb);
      ESTIMATE 'Other Race ki' (ora/orb);
	  predict pred out=Resid_CO;
run;

ods rtf file="C:\Users\JNguyen\Desktop\OP Inhibition\Diagnostics\DDVP Caucasian-Other Races.rtf" startpage=no;
%DiagnosticsCO(chemical=DDVP, model=2)
ods rtf close;

