*================================================================*
* SAS programmer: James Nguyen, USEPA                            *
*                                                                *
* Project: OP inhibition studies                                 *
*                                                                *
* Chemical: Chlorethoxyfos                                       *
*                                                                *
* 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";

/* 
*============================================================================*
* Programmer: Kelly Higgins, Exponent                                        *
* Templates provided by James Nguyen, US EPA                                 *
* Project: OP In Vitro Inhibition studies                                    *
*                                                                            *
* Study/Chemical: chlorethoxyfos                                             *
*                                                                            *
* Purpose: to show the use of PROC NLMIXED                                   *
*                                                                            *
* Note: instead of reanalyzing the data analyses in steps 1 and 2 to         *
*       obtain the data for the analysis in step 3, the data (derived        *
*       from steps 1 and 2 by Exponent) in the EXCEL (provided               *
*       by Exponent) are used to demonstrate the use of PROC NLMIXED         * 
*       for the data analysis in step 3.                                     *
*                                                                            *
* This compound has no replicate data.                                       *     
*                                                                            *
* Date: 2-13-20                                                              *
*===========================================================================*;

%Macro readdata(file=,out=,range=);
	Proc import datafile="C:\Users\JNguyen\OneDrive - Environmental Protection Agency (EPA)\Animal studies\50773501-2-3 In Vitro Inhibition Study\Chlorethoxyfos\&file"
		dbms=xlsx replace out=&out;
		range="&range";
		getnames=yes;
	quit;
%Mend;

%readdata(file=KineticsFemale Rat 1.xlsx, out=Female_Rat_1, range=Summary$P16:S22);
%readdata(file=KineticsFemale Rat 2.xlsx, out=Female_Rat_2, range=Summary$P16:S22);
%readdata(file=KineticsFemale Rat 3.xlsx, out=Female_Rat_3, range=Summary$P16:S22);
%readdata(file=KineticsH10yrCF1.xlsx,     out=H10yrCF1,     range=Summary$P16:S22);
%readdata(file=KineticsH10yrCF2.xlsx,     out=H10yrCF2,     range=Summary$P16:S22);
%readdata(file=KineticsH10yrCM.xlsx,      out=H10yrCM,      range=Summary$P16:S22);
%readdata(file=KineticsH11yrCF.xlsx,      out=H11yrCF,      range=Summary$P16:S22);
%readdata(file=KineticsH13yrCM.xlsx,      out=H13yrCM,      range=Summary$P16:S22);
%readdata(file=KineticsH16yrCF.xlsx,      out=H16yrCF,      range=Summary$P16:S22);
%readdata(file=KineticsH23yrHM.xlsx,      out=H23yrHM,      range=Summary$P16:S22);
%readdata(file=KineticsH27yrCF.xlsx,      out=H27yrCF,      range=Summary$P16:S22);
%readdata(file=KineticsH30yrAAF.xlsx,     out=H30yrAAF,     range=Summary$P16:S22);
%readdata(file=KineticsH31yrAAM.xlsx,     out=H31yrAAM,     range=Summary$P16:S22);
%readdata(file=KineticsH35yrHF.xlsx,      out=H35yrHF,      range=Summary$P16:S22);
%readdata(file=KineticsH38wks5dCM.xlsx,    out=H38wk5dCM,    range=Summary$P16:S22);
%readdata(file=KineticsH39wk4dCM.xlsx,    out=H39wk4dCM,    range=Summary$P16:S22);
%readdata(file=KineticsH39wks1dAAF.xlsx,  out=H39wks1dAAF,  range=Summary$P16:S22);
%readdata(file=KineticsH41wkCF.xlsx,      out=H41wkCF,      range=Summary$P16:S22);
%readdata(file=KineticsH46yrCM.xlsx,      out=H46yrCM,      range=Summary$P16:S22);
%readdata(file=KineticsH51yrCF.xlsx,      out=H51yrCF,      range=Summary$P16:S22);
%readdata(file=KineticsH60yrCM.xlsx,      out=H60yrCM,      range=Summary$P16:S22);
%readdata(file=KineticsMale Rat 1.xlsx,   out=Male_Rat_1,   range=Summary$P16:S22);
%readdata(file=KineticsMale Rat 2.xlsx,   out=Male_Rat_2,   range=Summary$P16:S22);
%readdata(file=KineticsMale Rat 3.xlsx,   out=Male_Rat_3,   range=Summary$P16:S22);

Data chlorethoxyfos;
	set H10yrCF1	(in=s1) 	
		H10yrCF2	(in=s2)		
		H10yrCM		(in=s3)	
		H11yrCF		(in=s4) 		
		H13yrCM		(in=s5)		
		H16yrCF		(in=s6)	
		H23yrHM		(in=s7)	
		H27yrCF		(in=s8)	
		H30yrAAF	(in=s9)
		H31yrAAM	(in=s10)	
		H35yrHF		(in=s11)	
		H38wk5dCM	(in=s12)
		H39wk4dCM	(in=s13)	
		H39wks1dAAF	(in=s14)
		H41wkCF		(in=s15)
		H46yrCM		(in=s16)		
		H51yrCF		(in=s17)		
		H60yrCM		(in=s18)

		Male_Rat_1	(in=s19)	
		Male_Rat_2	(in=s20)	
		Male_Rat_3	(in=s21)
		Female_Rat_1(in=s22)	
		Female_Rat_2(in=s23)
		Female_Rat_3(in=s24) ;

	length ID $15.;
if s1 =1  then do ;  ID = "H10yrCF1" ;  AnimalIdicator = 0; female=1; male=0; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=1;  end;
if s2 =1  then do ;  ID = "H10yrCF2" ;    AnimalIdicator = 0; female=1; male=0; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=1;  end;
if s3 =1  then do ;  ID = "H10yrCM" ;      AnimalIdicator = 0; female=0; male=1; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=1;  end;
if s4 =1  then do ;  ID = "H11yrCF" ;    AnimalIdicator = 0; female=1; male=0; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=1;  end;
if s5 =1  then do ;  ID = "H13yrCM" ;    AnimalIdicator = 0; female=0; male=1; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=1;  end;
if s6 =1  then do ;  ID = "H16yrCF" ;    AnimalIdicator = 0; female=1; male=0; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=1;  end;
if s7 =1  then do ;  ID = "H23yrHM" ;    AnimalIdicator = 0; female=0; male=1; AA=0; Hispanic=1; Caucasian=0; OtherRace=1; adult=1; infant=0; age2=2;  end;
if s8 =1  then do ;  ID = "H27yrCF" ;    AnimalIdicator = 0; female=1; male=0; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=2;  end;
if s9 =1  then do ;  ID = "H30yrAAF" ; AnimalIdicator = 0; female=1; male=0; AA=1; Hispanic=0; Caucasian=0; OtherRace=1; adult=1; infant=0; age2=2;  end;
if s10 =1 then do ;  ID = "H31yrAAM" ;    AnimalIdicator = 0; female=0; male=1; AA=1; Hispanic=0; Caucasian=0; OtherRace=1; adult=1; infant=0; age2=2;  end;
if s11 =1 then do ;  ID = "H35yrHF" ;    AnimalIdicator = 0; female=1; male=0; AA=0; Hispanic=1; Caucasian=0; OtherRace=1; adult=1; infant=0; age2=2;  end;
if s12 =1 then do ;  ID = "H38wk5dCM" ; AnimalIdicator = 0; female=0; male=1; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=0; infant=1; age2=0;  end;
if s13 =1 then do ;  ID = "H39wk4dCM" ; AnimalIdicator = 0; female=0; male=1; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=0; infant=1; age2=0;  end;
if s14 =1 then do ;  ID = "H39wks1dAAF" ; AnimalIdicator = 0; female=1; male=0; AA=1; Hispanic=0; Caucasian=0; OtherRace=1; adult=0; infant=1; age2=0;  end;
if s15 =1 then do ;  ID = "H41wkCF" ; AnimalIdicator = 0; female=1; male=0; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=0; infant=1; age2=0;  end;
if s16 =1 then do ;  ID = "H46yrCM" ;    AnimalIdicator = 0; female=0; male=1; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=2;  end;
if s17 =1 then do ;  ID = "H51yrCF" ;    AnimalIdicator = 0; female=1; male=0; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=2;  end;
if s18 =1 then do ;  ID = "H60yrCM" ; AnimalIdicator = 0; female=0; male=1; AA=0; Hispanic=0; Caucasian=1; OtherRace=0; adult=1; infant=0; age2=2;  end;

if s19 =1 then do ;  ID = "Male_Rat_1" ;  AnimalIdicator = 1 ; end; 
if s20 =1 then do ;  ID = "Male_Rat_2" ; AnimalIdicator = 1 ; end; 
if s21 =1 then do ;  ID = "Male_Rat_3" ; AnimalIdicator = 1 ; end; 
if s22 =1 then do ;  ID = "Female_Rat_1" ; AnimalIdicator = 1 ; end; 
if s23 =1 then do ;  ID = "Female_Rat_2" ; AnimalIdicator = 1 ; end; 
if s24 =1 then do ;  ID = "Female_Rat_3" ; AnimalIdicator = 1 ; end; 
	rename VAR4 = Inversed_I_M;
run;

Data epa.chlorethoxyfos;
	set chlorethoxyfos;
	if AnimalIdicator = 0 then do; H_I__M = _I__M; Human=1; end;
	else do; H_I__M = 0; Human=0;end;

	if AnimalIdicator = 1 then do; R_I__M = _I__M; Rat = 1; end;
	else do; R_I__M = 0; Rat = 0; end;
run;

*/


/****************************************************************************/
/* Final NLMIXED Model                                                      */
/****************************************************************************/
data Chlorethoxyfos;
     set epa.Chlorethoxyfos;
run;
proc sort data=Chlorethoxyfos;
     by id;
run;

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


*===> Human vs. Rat;

title1 'Human vs. rat: Chlorethoxyfos Modelling - Final Model 2';
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=Chlorethoxyfos /*itdetails*/ tech=newrap lis=3;
      PARMS ha=0.9236
			hb=3.6790
			ra=0.6739		
			rb=1.2808
			ga=4.3744	
			gab=7.6246
			gb=1.9
			S2E=1.6695/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-8 + H_I__M + rrb*Rat*1e-8 + R_I__M);
      MODEL KAPP ~ NORMAL(PRED, S2E*1e-4);
      RANDOM  ua ub ~ normal([0, 0],  [ga/100, gab/100, gb/10]) 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\Chlorethoxyfos Human-Rat.rtf" startpage=no;
%DiagnosticsHR(chemical=Chlorethoxyfos, model=2)
ods rtf close;
%Calculate_Ki(chemical=Chlorethoxyfos, Hadj=1e-8, Radj=1e-8);


*===> 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=chlorethoxyfos(where=(Human=1)) tech=nrridg lis=3 maxiter=1000 method=firo;
      PARMS ada=.8714	adB=.3706
			ia=1.2998 iB=.4712
			gha=.1229 	ghab=.6217 	ghb=.3666
			S2E=0.9714;

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


*===> Males vs. females;


title "final model 2b - male vs. female";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=chlorethoxyfos(where=(Human=1)) method=firo tech=newrap lis=3 maxiter=100000;
      PARMS ma=.90 mb= .37
			fa=1.05 fB=.43
			gha=.165 	ghab=.77 	ghb=.44
			S2E=.9542;

            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-7 + H_I__M*male + rfb*female*1e-7 + H_I__M*female);
      MODEL KAPP ~ NORMAL(PRED, S2E*1e-4);
      RANDOM  uha uhb ~ normal([0, 0],  [gha, ghab*1e-1, ghb*1e-1]) 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\Chlorethoxyfos Male-Female.rtf" startpage=no;
%DiagnosticsMF(chemical=Chlorethoxyfos, model=2)
ods rtf close;


*===> Caucasian vs. Others;


title "model 2 - Caucasian vs. Other Race";
*ods select  ParameterEstimates  FitStatistics ConvergenceStatus  AdditionalEstimates;
ods output ParameterEstimates=parm;
PROC NLMIXED data=chlorethoxyfos(where=(Human=1)) method=firo tech=newrap lis=3 maxiter=100000;
      PARMS ca=1.0034	cB=.4047
			ora=.9018 orB=.3885
			gha=1.6280	ghab=.7633	ghb=.4335
			S2E=.9533;

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