*SAS programme code for modelling ASSOCIATIONS WITH ARCTIC SKUA POPULATION TRENDS, using GLMMs with Poisson errors*;

*INTER-CORRELATIONS BETWEEN PREDICTORS*;
proc corr data=work.astrendvprod;var ASprod4yrlag lnGS1km lnGS1kmperAS;run;

*MODEL 1 - EFFECT OF ARCTIC SKUA PRODUCTIVITY FOUR YEARS PREVIOUSLY*;
proc glimmix data=work.astrendvprod plots=(PearsonPanel);where lnGS1km is not null;*to make datasets the same for all three models*;
class site;
model AScount/ASsum = ASprod4yrlag / D=b SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*MODEL 2 - AS ABOVE, BUT WITH THE ADDITION OF GREAT SKUA POPULATION SIZE*;
proc glimmix data=work.astrendvprod plots=(PearsonPanel);
class site;
model AScount/ASsum = ASprod4yrlag lnGS1km / D=b SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*MODEL 3 - AS MODEL 1, BUT WITH THE ADDITION OF GREAT SKUA:ARCTIC SKUA RATIO*;
proc glimmix data=work.astrendvprod plots=(PearsonPanel);
class site;
model AScount/ASsum = ASprod4yrlag lnGS1kmperAS / D=b SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*FITTING MODELS WITH LAPLACE TO GET MODEL FIT VALUES FOR COMPARISON BETWEEN MODELS*;
*MODEL 1*;
proc glimmix data=work.astrendvprod method=laplace;where lnGS1km is not null;*to make datasets the same for all three models*;
class site;
model  AScount/ASsum = ASprod4yrlag / solution dist=bin link=logit;
random site;
run;*AICC = 729.55, -2LogLikelihood = 725.45*;

*MODEL 2*;
proc glimmix data=work.astrendvprod method=laplace;
class site;
model  AScount/ASsum = ASprod4yrlag lnGS1km / solution dist=bin link=logit;
random site;
run;*AICC = 730.90, -2LogLikelihood = 724.71*;

*MODEL 3*;
proc glimmix data=work.astrendvprod method=laplace;
class site;
model  AScount/ASsum = ASprod4yrlag lnGS1kmperAS / solution dist=bin link=logit;
random site;
run;*AICC = 731.40, -2LogLikelihood = 725.21*;

*NULL MODEL*;
proc glimmix data=work.astrendvprod method=laplace;where lnGS1km is not null;*to make datasets the same for all three models*;
class site;
model  AScount/ASsum =  / solution dist=bin link=logit;
random site;
run;*NULL MODEL - AICC = 889.91, -2LogLikelihood = 887.89*;


*TESTING FOR OVERALL ASSOCIATION BETWEEN ARCTIC SKUA POPULATION TREND AND GREAT SKUA POPULATION SIZE*;
proc glimmix  data=work.astrend plots=(PearsonPanel);
class sitecode;
model AS2 = lnS2000 lat2 year2 year2*year2 lngs1km year2*lngs1km / D=p SOLUTION DDFM=KENWARDROGER;
random sitecode;
random _resid_;
run;