* Created by: Jay Rotella
* Date: June 24, 2003
* Purpose:
*   This file:
*	1. inputs a dataset containing nest survival information,
*	2. calculates the effective sample size for the dataset,
*	3. runs a variety of models of NEST SURVIVAL in NLMIXED,
*	4. creates an AICc table for model selection, & 
*	5. outputs the AICc table to HTML, RTF, and pdf files.


*  Step 1: calculate the effective sample size according to the methods
*  of Dinsmore et al. (2002).  Here, n-ess is incremented by 1 for
*  each day a nest was under observation and survived and by 1 for
*  each interval for which a nest was under observation and failed.;

* This step:
*   1. calculates the contribution to n-ess for each observation interval and
*		adds that contribution to the sum of ness, i.e., ness column is a running total
*	2. creates dummy/indicator variables for each of 4 habitat types;
data Mall; set Sasuser.mall2000nd;
 if ifate=0 then ness+1;
 else if ifate=1 then ness+t;
/* create indicator variables for different nesting habitats */
    if hab=1 then NatGr=1; else NatGr=0;  				/* Native Grassland */
    if hab=2 or hab=3 or hab=9 then CRP=1; else CRP=0;  /* CRP & similar    */
    if hab=7 or hab=22 then Wetl=1; else Wetl=0;        /* Wetland sites    */
    if hab=20 then Road=1; else Road=0;					/* Roadside sites   */
run;

* This step finds the actual n-ess for the dataset, 
*   which is the maximum value in the ness column.;
Proc Univariate data=Mall;
 var ness;
 output out=ness max=ness;
run;

* This step sorts the data by site, which is used as a random factor in some models below.
* PROC NLMIXED assumes that a new realization occurs whenever the SUBJECT= variable changes
*  from the previous observation, so your input data set should be clustered according to
*  this variable. One easy way to accomplish this is to run PROC SORT prior to calling 
*  PROC NLMIXED using the SUBJECT= variable as the BY variable. ;
Proc Sort data=Mall;
by site; run;

* This step reformats the Fit Statistics table of NLMIXED
*   so it displays more decimal places in the created tables;
proc template;                                                                
   define table Stat.Nlm.FitStatistics;                                       
      notes "Fit statistics";                                                 
      column Descr Value;                                                     
      header H1;                                                              
      define H1;                                                              
         text "Fit Statistics";                                               
         space = 1;                                                           
      end;                                                                    
      define Descr;                                                           
         header = "Description";                                              
         width = 30;                                                          
         print_headers = OFF;                                                 
         flow;                                                                
      end;                                                                    
      define Value;                                                           
         header = "Value";                                                    
         format = 12.4;                                                        
         print_headers = OFF;                                                 
      end;                                                                    
   end;                                                                       
run;

* This macro merges various pieces of information from NLMIXED output
*   and calculates the AICc value based on ness rather than on the number
*   of observations in the dataset;
%Macro AICtable(out=,dataset=,mod=);

	/****************************************/
	/*  Store information for output file	*/
	/*	1. AIC value from NLMIXED   		*/
	/*	2. Dataset name (from macro call)	*/
	/*	3. Model name (from macro call)		*/
	/*	4. Convergence status				*/
	/*	5. AICc value based on n-ess		*/
	/*	6. n = # obs used for model			*/
	/*	7. k = # parameters in model		*/
	/****************************************/

	data temp2;
	  set converge;
		length dataset $10.;
		dataset="&dataset";
		length mod $30.;
		mod="&mod";
		merge converge n k ness aic;
		if Reason='NOTE: GCONV convergence criterion satisfied.' 
			then Conv='Yes'; else Conv='No';
		length Conv $5.;
		aicc=aic+(2*k)*((k+1)/(ness-k-1));
		keep dataset Conv mod n k ness aic aicc;
		run;
	/********************************************************************/  
	/* Append the information (see previous step) to a file that    	*/
	/* holds this information for all models to be run in this batch.	*/
	/********************************************************************/  
 
	Proc append base=&out data=temp2 force;
	run;

proc print data=&out; run;

%Mend AICtable; 
	/********************************************************************/
	/* model statements made based on X-vars (1) patch (2) age (3) pred	*/
	/*	and Y = ifate													*/
	/*																	*/
	/*	Need to change info in 3 places	for each model run				*/
	/*		1. parms - add parameters (& start values if want/need to)	*/
	/*		2. logit equation - put in proper model statement			*/
	/*		3. macro call for "mod=" - put in text label for the model	*/
	/*																	*/
	/*																	*/	
	/*	Need to change info in 1 place the 1st time a model is run		*/
	/*		1. macro call for "data=" - give a text label for dataset 	*/
	/********************************************************************/

* Run a model of interest;
* This is where you can add, delete, or alter models to suit your needs;
Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=B0;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
* call the AICtable macro to append the key info to a results table;
%AICtable(out=results,dataset=Mall,mod=B0) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B1=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=B0+B1*(sdate+i);
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B1*Date) run; 
Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B2=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=B0+B2*(sage+i);
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B2*Age) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B3=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=B0+B3*Robel;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B3*Robel) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B4=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=B0+B4*PctGr4;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B4*PctGr4) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B5=0, B6=0, B7=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=B0+B5*NatGr+B6*CRP+B7*Wetl;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B5*NatGr+B6*CRP+B7*Wetl) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B2=0, B3=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=B0+B2*(sage+i)+B3*Robel;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B2*Age+B3*Robel) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B2=0, B4=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=B0+B2*(sage+i)+B4*PctGr4;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B2*Age+B4*PctGr4) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B2=0, B5=0, B6=0, B7=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=B0+B2*(sage+i)+B5*NatGr+B6*CRP+B7*Wetl;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B2*Age+B5*NatGr+B6*CRP+B7*Wetl) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B2=0, B4=0, B8=0;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=(B0+(B8*Ob))+B2*(sage+i)+B4*PctGr4;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B2*Age+B4*PctGr4+B8*Ob) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=2.42, B2=0.019, B4=0.38, vsite=0.1;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=(B0+u)+B2*(sage+i)+B4*PctGr4;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
random u~normal(0,vsite) subject=site;
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B2*Age+B4*PctGr4+b1*site) run; 

Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=2.42, B2=0.019, B4=0.38, B8=-1, vsite=0.1;
	p=1;
	   do i=0 TO t-1;
	   if i=0 then Ob=1;
	   else Ob=0;
	   	   logit=(B0+u)+B2*(sage+i)+B4*PctGr4+B8*Ob;
	      p=p*(exp(logit)/(1+exp(logit)));
	   end;
model ifate~binomial(1,p);
random u~normal(0,vsite) subject=site;
ods output ConvergenceStatus=converge;
ods output FitStatistics=aic 
	   (where=(Descr='AIC (smaller is better)') 
			   rename=(Value=AIC)) ;
ods output Dimensions=n
		(where=(Descr='Observations Used')
				rename=(Value=n));
ods output Dimensions=k
		(where=(Descr='Parameters')
				rename=(Value=k));
RUN;
%AICtable(out=results,dataset=Mall,mod=B0+B2*Age+B4*PctGr4+B8*Ob+b1*site) run; 


* The following statements and procedures calculate
*   1. minimumAICc values,
*   2. deltaAICc values, 
*   3. model weights, &
*   4. model likelihoods.;
Proc sort data=work.results;
 by AICc;
run;
Proc Univariate data=work.results;
 var AICc;
 output out=minaicc min=MinAICc;
run;
 Data AICctemp;
 if _n_=1 then set minaicc;
 set work.results;
 DeltaAICc=AICc-minaicc;
    Odds = exp(0.5*DeltaAICc);
    InvOdds = 1/Odds;
Run;
Proc Univariate data=work.AICctemp;
 var InvOdds;
 output out=invoddssum sum=sum;
run;
Data Weighttemp;
 if _n_=1 then set invoddssum;
 set work.AICctemp;
   Weight = InvOdds/sum;
Run;
Proc Univariate;
 var Weight;
 output out=BestWeight max=BestWeight;
run;
Data ModLiktemp;
 if _n_=1 then set BestWeight;
 set work.Weighttemp;
    ModLik=Weight/BestWeight;
Run;
 

	/********************************************************************/
	/* Create an output file (RTF) file that contains the key output    */
	/*	You can change:                                                 */
	/*																	*/
	/*	Can change the output's:										*/
	/*		1. formats so that variables appear as desired              */
	/*		2. file destination                                         */
	/*		3. file name                                             	*/
	/*		4. file type (can turn on or off 'RTF' to 'PDF' or 'HTML')  */
	/*																	*/
	/********************************************************************/

ods rtf File='C:\Documents and Settings\Rotella\My Documents\research projects\nest success ducks\FixMixEffOut.rtf';
* ods pdf File='C:\Documents and Settings\Rotella\My Documents\research projects\nest success ducks\output.pdf';
* ods html File='C:\Documents and Settings\Rotella\My Documents\research projects\nest success ducks\output.html';
Data Final;
	set ModLiktemp;
	drop sum MinAICc BestWeight Odds InvOdds;
	label ness='n-ess';
	label aicc='AICc';
   file print ods=(variables=(dataset
              	mod
				Conv
				n
				k
				ness
				AIC
				AICc
				DeltaAICc
				Weight
				ModLik));
      format dataset $5.;
	  format mod $28.;
	  format Conv $3.;
	  format n 6.0;
	  format k 2.0;
	  format ness 9.0;
	  format AIC 9.3;
	  format AICc 9.3;
      format DeltaAICc 8.3;
	  format Weight 5.3;
	  format ModLik 5.3;
   put _ods_;
run;
 ods rtf close;
 ods pdf close;
 ods html close;