/**********************************************************************************/
/* Name: logexpos.sas                                                             */
/* Author: Terry Shaffer (terry_shaffer@usgs.gov)                                 */
/* SAS Version: 8                                                                 */
/* Supporting files: aic_mac.sas                                                  */
/* Date: 4 Dec 2002                                                               */
/*                                                                                */
/* This program illustrates 1) fitting of logistic-exposure nest-survival         */
/* models with PROC GENMOD, 2) computation of AIC model-selection criteria,       */
/* and 3) computation of model-averaged regression coefficients and unconditional */
/* standard errors. The example involves evaluation of eight candidate models.    */
/*                                                                                */
/**********************************************************************************/

/* Include the file containing the macros for computing model-selection criteria  */
/* and model-averaged regression coefficients. The user will need to modify the   */
/* %Include statement to point to the location of the aic_mac.sas file on their   */
/* computer.  Detailed descriptions of the macros are given in the macro file.    */
%Include "i:\klett\macros\aic_mac.sas";

/* Read in one observation for each interval of exposure on each nest.    */
/* expos = interval length; survive=1 (if the nest survives the interval, */
/* survive=0 otherwise; nest_ht is a continuous explanatory variable;     */
/* parastat and patsize are categorical explanatory variables.            */
data chat;
 length patsize $ 5;
 input expos nest_ht parastat patsize survive;
 if nest_ht=. then delete;
 trials=1;
 cards;
1 0.8 1 small 1
3 0.8 1 small 1
6 0.8 1 small 1
3 0.8 1 small 1
2 0.8 1 small 1
3 0.8 1 small 1
2 0.8 1 small 1
1 0.9 0 large 1
1 0.9 0 large 1
2 0.9 0 large 1
1 0.9 0 large 1
1 0.9 0 large 1
1 0.9 0 large 1
1 0.9 0 large 1
5 0.9 0 large 1
1 0.9 0 large 1
1 0.9 0 large 0
1 0.6 0 small 1
2 0.6 0 small 0
2 . 1 small 1
3 . 1 small 1
3 . 1 small 1
2 . 1 small 1
2 . 1 small 1
1 . 1 small 0
1 0.3 1 small 1
1 0.3 1 small 1
2 0.3 1 small 1
2 0.3 1 small 1
4 0.3 1 small 1
2 0.3 1 small 1
1 0.3 1 small 0
1 . 0 large 1
1 . 0 large 1
1 . 0 large 1
1 1.7 1 large 1
1 1.7 1 large 1
1 1.7 1 large 1
1 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
2 0.5 1 small 1
2 0.5 1 small 1
5 0.5 1 small 1
3 0.5 1 small 1
5 0.5 1 small 1
3 0.5 1 small 1
1 0.6 1 large 1
5 0.6 1 large 1
3 0.6 1 large 1
1 0.6 1 large 1
3 0.6 1 large 1
4 0.6 1 large 1
3 0.6 1 large 1
2 0.5 0 large 1
1 0.5 0 large 1
2 0.5 0 large 1
1 0.5 0 large 1
1 0.5 0 large 1
1 0.5 0 large 1
1 0.5 0 large 1
4 0.5 0 large 1
3 0.5 0 large 0
1 0.5 1 large 1
1 0.5 1 large 1
1 0.45 1 small 1
3 0.45 1 small 1
1 0.45 1 small 1
1 0.45 1 small 1
1 0.45 1 small 1
1 0.45 1 small 1
4 0.45 1 small 1
1.5 0.45 1 small 0
4 0.5 1 large 1
4 0.5 1 large 1
1 0.5 1 large 1
1 0.5 1 large 1
4 0.5 1 large 1
6 0.5 1 large 1
1 0.5 1 large 1
1 0.5 1 large 1
1 0.5 1 large 1
6 0.7 0 small 1
1 0.7 0 small 1
1 0.7 0 small 0
2 0.7 1 large 1
3 0.7 1 large 1
1 0.7 1 large 1
2 0.7 1 large 1
1 0.7 1 large 1
1 0.7 1 large 1
4 0.61 0 large 1
2 0.61 0 large 1
2.5 0.61 0 large 0
2 1.52 1 small 1
2 1.52 1 small 1
1 1.52 1 small 1
2 1.52 1 small 1
1 1.52 1 small 0
2 0.5 1 small 1
1 0.5 1 small 1
1.5 0.5 1 small 0
4 0.55 0 large 1
4 0.55 0 large 1
4 0.55 0 large 1
1 0.55 0 large 1
1 0.55 0 large 1
2 0.55 0 large 1
1 0.55 0 large 1
1 0.55 0 large 1
0.5 0.55 0 large 0
2 0.77 0 large 1
6 0.77 0 large 1
1 0.77 0 large 1
2 0.77 0 large 1
1 0.77 0 large 1
1 0.77 0 large 1
0.5 0.77 0 large 0
1 0.82 0 large 1
2 0.82 0 large 1
3 0.82 0 large 1
2 0.82 0 large 1
2 0.82 0 large 1
1 0.82 0 large 0
1 0.45 1 small 1
1 0.45 1 small 1
1 0.45 1 small 1
3 0.45 1 small 1
1 0.45 1 small 1
1 0.45 1 small 1
1.5 0.45 1 small 0
2 0.6 1 large 1
1 0.6 1 large 1
1 0.6 1 large 1
1.5 0.6 1 large 0
1 0.8 1 large 1
2 0.8 1 large 1
1 0.8 1 large 1
1 0.8 1 large 1
1 0.8 1 large 1
2 0.8 1 large 1
1 0.8 1 large 1
1 0.8 1 large 1
1 0.8 1 large 1
1 0.8 1 large 1
1 0.8 1 large 1
1 0.8 1 large 1
1 0.8 1 large 1
1 0.28 1 small 1
3 0.28 1 small 1
1 0.28 1 small 1
1 0.28 1 small 1
1 0.28 1 small 1
1 0.28 1 small 0
2 0.42 1 small 1
2 0.42 1 small 1
2 0.42 1 small 1
2 0.42 1 small 1
2.5 0.42 1 small 0
2 0.5 1 small 1
2 0.5 1 small 1
2 0.5 1 small 1
2 0.5 1 small 1
5 0.5 1 small 1
2 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
2 0.5 1 small 1
2 0.7 1 small 1
4 0.7 1 small 1
1 0.7 1 small 1
2 0.7 1 small 1
3 0.7 1 small 1
1 0.7 1 small 1
2 0.7 1 small 1
1 0.7 1 small 1
1 0.8 0 large 1
3 0.8 0 large 1
1 0.8 0 large 1
0.5 0.8 0 large 0
1 0.53 1 small 1
1 0.53 1 small 1
2 0.53 1 small 1
1 0.53 1 small 1
1 0.53 1 small 1
1.5 0.53 1 small 0
1 0.9 1 large 1
1 0.9 1 large 1
1 0.9 1 large 1
1 0.9 1 large 1
2 0.9 1 large 1
1.5 0.9 1 large 0
2 1.2 1 small 1
2 1.2 1 small 1
2 1.2 1 small 1
1 1.2 1 small 1
1 1.2 1 small 1
1 0.7 1 large 1
1 0.7 1 large 1
2 0.7 1 large 1
2 0.7 1 large 1
3 0.7 1 large 1
2 0.7 1 large 1
1.5 0.7 1 large 0
1 1.05 1 small 1
1 1.05 1 small 1
0.5 1.05 1 small 0
2 0.9 1 large 1
2 1.3 1 small 1
0.5 1.3 1 small 0
1 0.9 0 large 1
1 0.9 0 large 1
1 0.9 0 large 1
1 0.9 0 large 1
2 0.9 0 large 1
1 0.9 0 large 1
1 0.9 0 large 1
2 0.9 0 large 1
2 0.9 0 large 1
2 0.9 0 large 1
2 0.9 0 large 1
2 0.9 0 large 1
1 0.9 0 large 1
1 0.85 0 large 1
1 0.85 0 large 1
1 0.85 0 large 1
1 0.85 0 large 1
1 0.85 0 large 1
1 0.85 0 large 1
1 0.85 0 large 1
2 0.85 0 large 1
4 0.85 0 large 1
1 0.85 0 large 1
1 0.85 0 large 1
2 0.85 0 large 1
2 0.85 0 large 1
3 0.85 0 large 1
2 0.85 0 large 1
1 0.85 0 large 1
1 0.45 1 large 1
1 0.45 1 large 1
1 0.45 1 large 1
1 0.45 1 large 1
1 0.45 1 large 1
3 0.45 1 large 1
3 0.45 1 large 1
2 0.45 1 large 1
2 0.45 1 large 1
2 0.45 1 large 1
1 0.45 1 large 1
2 0.45 1 large 1
1 0.45 1 large 1
1 0.45 1 large 1
1 0.45 1 large 1
1 0.45 1 large 1
1 0.75 1 large 1
2 0.75 1 large 1
2 0.75 1 large 1
2 0.75 1 large 1
1 0.85 1 large 1
1 0.85 1 large 1
1 0.85 1 large 1
1 0.85 1 large 1
1 0.85 1 large 1
3 0.85 1 large 1
0.5 0.85 1 large 0
1 0.5 1 small 1
1 0.5 1 small 1
2 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
2 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
2 0.5 1 small 1
1 0.5 1 small 1
3 0.5 1 small 1
3 0.5 1 small 1
1 0.5 1 small 1
2 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
1 0.5 1 small 1
1 0.3 1 small 1
1 0.3 1 small 1
1 0.3 1 small 1
1 0.3 1 small 1
1 0.3 1 small 0
1 0.45 1 small 1
1 0.45 1 small 1
1 0.45 1 small 1
1 0.45 1 small 1
1 0.45 1 small 0
1 0.45 0 small 1
1 0.45 0 small 1
0.5 0.45 0 small 0
1 . 1 small 1
1 . 1 small 1
1 . 1 small 1
1 . 1 small 1
2 . 1 small 1
2 . 1 small 1
2 . 1 small 1
3 . 1 small 1
2 . 1 small 1
1 . 1 small 1
;

/* following code computes the effective sample size for computing AICc */
/* Note: See Rotella at al. (2004) */
/* <http://www.montana.edu/rotella/nestsurv/Rotella%20et%20al-NestSurv.pdf> */
/* for a discussion of effective sample size */
Data N_eff(keep=n_eff);
  Set Chat end=lastobs;
  if survive=0 then n_eff+1;
  else if survive=1 then n_eff+expos;
  if lastobs then do;
    put "The effective sample size for computing AICc is " n_eff ;
    output;
  end;
run;      

/* following code fits a logistic-exposure constant-survival model    */
/* and creates three data sets containing information about the model */
/* and the results.                                                   */
proc genmod data=chat;
  a=1/expos;
  fwdlink link = log((_mean_**a)/(1-_mean_**a));
  invlink ilink = (exp(_xbeta_)/(1+exp(_xbeta_)))**expos;
  model survive/trials = / dist=bin;
  ods output modelfit=modelfit;
  ods output modelinfo=modelinfo;
  ods output ParameterEstimates=ParameterEstimates;
  title 'logistic-exposure, constant-survival model';
run;

/* following code invokes the aicc macro to read the data sets created by GENMOD */
/* and compute aic values.  AIC results will be stored in the data set named in  */
/* the dsn= statement.  Parameter estimates will be stored in the data set named */
/* in the estdsn= statement. model_dimension is set to 2 because the highest     */
/* order interaction in the suite of candidate models is 2 (parastat*patsize).   */ 
%aicc(dsn=chat_aic,estdsn=estimate,model=constant survival,model_dimension=2);

/* following code fits a logistic-exposure model with a parastat main effect */
proc genmod data=chat;
  class parastat;
  a=1/expos;
  fwdlink link = log((_mean_**a)/(1-_mean_**a));
  invlink ilink = (exp(_xbeta_)/(1+exp(_xbeta_)))**expos;
  model survive/trials = parastat / dist=bin;
  ods output modelfit=modelfit;
  ods output modelinfo=modelinfo;
  ods output ParameterEstimates=ParameterEstimates;
  title 'logistic-exposure, parasitism status main effect';
run;

/* Summarize and store the results from the above model */
%aicc(dsn=chat_aic,estdsn=estimate,model=parastat main effect,model_dimension=2);

/* following code fits a logistic-exposure model with a patsize main effect */
proc genmod data=chat;
  class patsize;
  a=1/expos;
  fwdlink link = log((_mean_**a)/(1-_mean_**a));
  invlink ilink = (exp(_xbeta_)/(1+exp(_xbeta_)))**expos;
  model survive/trials = patsize / dist=bin;
  ods output modelfit=modelfit;
  ods output modelinfo=modelinfo;
  ods output ParameterEstimates=ParameterEstimates;
  title 'logistic-exposure, patch size main effect';
run;

/* Summarize and store the results from the above model */
%aicc(dsn=chat_aic,estdsn=estimate,model=patsize main effect,model_dimension=2);

/* following code fits a logistic-exposure model with a parastat and patsize main effects */
proc genmod data=chat;
  class parastat patsize;
  a=1/expos;
  fwdlink link = log((_mean_**a)/(1-_mean_**a));
  invlink ilink = (exp(_xbeta_)/(1+exp(_xbeta_)))**expos;
  model survive/trials = parastat patsize / dist=bin;
  ods output modelfit=modelfit;
  ods output modelinfo=modelinfo;
  ods output ParameterEstimates=ParameterEstimates;
  title 'logistic-exposure, parastat and patsize main effect';
run;

/* Summarize and store the results from the above model */
%aicc(dsn=chat_aic,estdsn=estimate,model=parastat and patsize main effect,model_dimension=2);

/* following code fits a logistic-exposure model with a parastat and patsize main effects */
/* and their interaction                                                                  */
proc genmod data=chat;
  class parastat patsize;
  a=1/expos;
  fwdlink link = log((_mean_**a)/(1-_mean_**a));
  invlink ilink = (exp(_xbeta_)/(1+exp(_xbeta_)))**expos;
  model survive/trials = parastat patsize parastat*patsize / dist=bin;
  ods output modelfit=modelfit;
  ods output modelinfo=modelinfo;
  ods output ParameterEstimates=ParameterEstimates;
  title 'logistic-exposure, main effects and interaction';
run;

/* Summarize and store the results from the above model */
%aicc(dsn=chat_aic,estdsn=estimate,model=main effects and interaction,model_dimension=2);

/* following code fits a logistic-exposure model for nest_ht, a continuous covariate.     */
proc genmod data=chat;
  a=1/expos;
  fwdlink link = log((_mean_**a)/(1-_mean_**a));
  invlink ilink = (exp(_xbeta_)/(1+exp(_xbeta_)))**expos;
  model survive/trials = nest_ht / dist=bin;
  ods output modelfit=modelfit;
  ods output modelinfo=modelinfo;
  ods output ParameterEstimates=ParameterEstimates;
  title 'logistic-exposure, nest_ht';
run;

/* Summarize and store the results from the above model */
%aicc(dsn=chat_aic,estdsn=estimate,model=nest height only,model_dimension=2);

/* following code fits a logistic-exposure model with effects of parastat, patsize  */
/* and nest_ht.                                                                     */
proc genmod data=chat;
  class parastat patsize;
  a=1/expos;
  fwdlink link = log((_mean_**a)/(1-_mean_**a));
  invlink ilink = (exp(_xbeta_)/(1+exp(_xbeta_)))**expos;
  model survive/trials = parastat patsize nest_ht / dist=bin;
  ods output modelfit=modelfit;
  ods output modelinfo=modelinfo;
  ods output ParameterEstimates=ParameterEstimates;
  title 'logistic-exposure, main effects and nest-height';
run;

/* Summarize and store the results from the above model */
%aicc(dsn=chat_aic,estdsn=estimate,model=main effects and nest-height,model_dimension=2);

/* following code fits a logistic-exposure model with effects of parastat, patsize  */
/* parastat*patsize, and nest_ht.                                                                     */
proc genmod data=chat;
  class parastat patsize;
  a=1/expos;
  fwdlink link = log((_mean_**a)/(1-_mean_**a));
  invlink ilink = (exp(_xbeta_)/(1+exp(_xbeta_)))**expos;
  model survive/trials = parastat patsize parastat*patsize nest_ht / dist=bin;
  ods output modelfit=modelfit;
  ods output modelinfo=modelinfo;
  ods output ParameterEstimates=ParameterEstimates;
  title 'logistic-exposure, main effects, interaction, and nest-height';
run;

/* Summarize and store the results from the above model */
%aicc(dsn=chat_aic,estdsn=estimate,model=main effects + interaction + nest-height,model_dimension=2);

/* Compute AIC model selection criteria and Akaike weights*/
%DeltaAic(DataIn=chat_aic,DataOut=chat_delta_aicc,Vari=aicc);

/* Compute model-averaged estimates and unconditional standard errors */
%Modelavg(estimate=estimate, Akaike=chat_delta_aicc, DataOut=Model_Averaged_Estimates, model_dimension=2);