*================================================================*
* SAS programmer: James Nguyen, USEPA                            *
*                                                                *
* Project: OP inhibition studies                                 *
*                                                                *
* Chemical: methamidophos                                        *
*                                                                *
* Descriptions:                                                  *
*     - run the selected (Exponent) final models                 *
*     - perform regression diagnostics                           *
*     - calculated ki values from human vs. rat analysis         *
*                                                                *
* Date: March 2020                                               *
*===============================================================*;
option formdlim="=" ls=100 ps=100 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 methamidophos;
     set epa.methamidophos;
run;
proc sort data=methamidophos;
     by id;
run;

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




*===> Human vs. Rat;

*===> the model below might be better, but it was not selected by Exponent;

title "Human vs. Rat: model 2";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=methamidophos  tech=newrap lis=3;
      PARMS ha=.8450		hb=.7677
			ra=0.5998		rb=.4709
			ga=0.1316	
			gab=.1163
			gb=.1175
			S2E=0.8187/best=1000;

            rha=ha+ua;
			rra=ra+ua;
            rhb=hb+ub;
			rrb=rb+ub;
            PRED = (rha*H_I__M + rra*R_I__M)/(rhb*Human*1e-3 + H_I__M + rrb*Rat*1e-3 + R_I__M);
      MODEL KAPP ~ NORMAL(PRED, S2E*1e-4);
      RANDOM  ua ub ~ normal([0, 0],  [ga, gab, gb]) 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;

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


/*
title1 'Methamidophos Modelling - Exponent Final Model 3';
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=Methamidophos method=firo tech=newrap lis=3 maxiter=10000 gconv=0;
      PARMS ha=0.7562
			hb=0.6833	
			ra=0.3921
			rb=0.20442	
			gha=0.2042	ghab=.2015		ghb=.2206
			gra=.1690	grab=0.1949	grb=.2461
			S2E=0.8182;

            rha=ha+uha;
			rra=ra+ura;
            rhb=hb+uhb;
			rrb=rb+urb;
            PRED = (rha*H_I__M + rra*R_I__M)/(rhb*Human*1e-3 + H_I__M + rrb*Rat*1e-3 + R_I__M);
      MODEL KAPP ~ NORMAL(PRED, S2E*1e-4);
      RANDOM  uha uhb ura urb ~ normal([0, 0, 0, 0],	[gha, 
												 		 ghab, ghb,
														 0,    	  0,   gra*1e-1,
														 0,       0,   grab*1e-1, grb*1e-1]) 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;

ods rtf file="C:\Users\JNguyen\Desktop\OP Inhibition\Diagnostics\Methamidophos human-rat (Exponent model).rtf" startpage=no;
%DiagnosticsHR(chemical=Methamidophos, model=3)
ods rtf close;
*/


*===> 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=methamidophos(where=(Human=1))  tech=nrridg lis=3 maxiter=1000 method=firo;
      PARMS ada=0.7067	adB=6.4950
			ia=.9285 iB=7.8211
			gha=0.1906 	ghab=1.9501	ghb=21.6014
			S2E=0.7164;

			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-4 + H_I__M*adult + rib*1e-4*infant + H_I__M*infant);
      MODEL KAPP ~ NORMAL(PRED, S2E*1e-4);
      RANDOM  uha uhb ~ normal([0, 0],  [gha, ghab, ghb]) 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\Methamidophos Adult-Infant.rtf" startpage=no;
%DiagnosticsAI(chemical=Methamidophos, model=2);
ods rtf close;


*===> Male vs. Female;

*===> the model below might be better, but it was not selected by Exponent;
title "final model 2b - male vs. female";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=methamidophos(where=(Human=1)) method=firo  tech=newrap lis=3 maxiter=100000;
      PARMS ma=.6603 mb= 5.7324
			fa=.8599 fB=7.8420
			gha=2.1590 	ghab=20.8888	ghb=21.7683
			S2E=.07074;

            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*1e-1, ghab*1e-1, 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\Methamidophos Male-Female.rtf" startpage=no;
%DiagnosticsMF(chemical=Methamidophos, model=2);
ods rtf close;

/*
title "Final model 3 - male vs. female";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=methamidophos(where=(Human=1)) method=firo  tech=newrap lis=3 maxiter=100000;
      PARMS ma=0.5767	mB=5.2115
			fa=0.9513 fB=8.5214
			gma=0.1533 	gmab=3.0989 	gmb=7.2074
			gfa=4.8876 	gfab=44.2878	gfb=40.3962
			S2E=.07045/ best=1000;
         	rma=ma+uma;
            rmb=mB+umb;
			rfa=fa+ufa;
			rfb=fB+ufb;
            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  uma umb ufa ufb ~ normal([0, 0, 0, 0],   [gma*1e-1, 
														gmab*1e-1, gmb, 
														0,    0,     gfa*1e-1, 
														0,    0,     gfab*1e-1,   gfb]) 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\Methamidophos Male-Female (exponent model).rtf" startpage=no;
%DiagnosticsMF(chemical=Methamidophos, model=3);
ods rtf close;
*/

*===> Caucasian vs. Other races;

title "Final model 2 - Caucasian vs. Other Race";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=methamidophos(where=(Human=1)) method=firo  tech=newrap lis=3 maxiter=100000;
      PARMS ca=0.8034	cB=7.7190
			ora=.6992 orB=5.6120
			gha=2.0131	ghab=19.1207	ghb=19.6158
			S2E=.07137;

			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-4 + H_I__M*caucasian + rorb*otherrace*1e-4 + H_I__M*otherrace);
      MODEL KAPP ~ NORMAL(PRED, S2E*1e-3);
      RANDOM  uha uhb ~ normal([0, 0],  [gha*1e-1, ghab*1e-1, ghb]) 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\Methamidophos Caucasian-Other Races.rtf" startpage=no;
%DiagnosticsCO(chemical=Methamidophos, model=2);
ods rtf close;

