*SAS programme code for modelling individual species population trends using GLMMs with Poisson errors*;

*ARCTIC SKUA POPULATION TRENDS*;

*Relates to FIGURE 2a - overall population trend (solid line). Year fitted as a covariate*;

proc glimmix data=work.astrend plots=(PearsonPanel);
class sitecode;
model AS2 = lnS2000 year2 year2*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random sitecode;
random _resid_;
run;

*Relates to FIGURE 2a - annual population size estimates. Year fitted as a categorical variable*;

proc glimmix data=work.astrend plots=(PearsonPanel);
class sitecode year2;
model AS2 = lnS2000 year2 / D=p SOLUTION DDFM=KENWARDROGER;
random sitecode;
random _resid_;
lsmeans year2 / ilink cl plot=mean(ilink);
run;

*TESTING FOR GEOGRAPHICAL VARIATION IN ARCTIC SKUA POPULATION TRENDS*;

proc glimmix data=work.astrend plots=(PearsonPanel);
class sitecode;
model AS2 = lnS2000 lat2 year2 year2*year2 lat2*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random sitecode;
random _resid_;
run;

*Relates to FIGURE S1 - overall population trend for each colony type (solid lines). Year fitted as a covariate*;

proc glimmix data=work.astrend plots=(PearsonPanel);
class sitecode type;
model AS2 = lnS2000 lat2 type year2 year2*year2 type*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random sitecode;
random _resid_;
run;

*Relates to FIGURE S1 - annual population size estimates for each colony type. Year fitted as a categorical variable*;

proc glimmix data=work.astrend plots=(PearsonPanel);
class sitecode type year2;
model AS2 = lnS2000 lat2 type year2 type*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random sitecode;
random _resid_;
lsmeans type*year2 / ilink cl plot=mean(ilink);
run;


*GREAT SKUA POPULATION TRENDS*;

*Relates to FIGURE 3a - overall population trend (solid line). Year fitted as a covariate*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.gstrend plots=(PearsonPanel);
class site;
model AOT2 = lnS2000 year2 year2*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*Relates to FIGURE 3a - annual population size estimates. Year fitted as a categorical variable*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.gstrend plots=(PearsonPanel);
class site year2;
model AOT2 = lnS2000 year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
lsmeans year2 / ilink cl plot=mean(ilink);
run;

*TESTING FOR GEOGRAPHICAL VARIATION IN GREAT SKUA POPULATION TRENDS*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.gstrend plots=(PearsonPanel);
class site;
model AOT2 = lnS2000 lat2 year2 year2*year2 lat2*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*Relates to FIGURE S2 - overall population trend for each colony type (solid lines). Year fitted as a covariate*;

proc glimmix maxopt=200 pconv=1.11022E-5 data=work.gstrend plots=(PearsonPanel);
class site type;
model AOT2 = lnS2000 lat2 type year2 year2*year2 type*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*Relates to FIGURE S2 - annual population size estimates for each colony type. Year fitted as a categorical variable*;

proc glimmix maxopt=200 pconv=1.11022E-5 data=work.gstrend plots=(PearsonPanel);
class site type year2;
model AOT2 = lnS2000 lat2 type year2 type*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
lsmeans type*year2 / ilink cl plot=mean(ilink);
run;


*BLACK-LEGGED KITTIWAKE POPULATION TRENDS*;

*Relates to FIGURE 3b - overall population trend (solid line). Year fitted as a covariate*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.kitrend plots=(PearsonPanel);
class site;
model AOT2 = lnS2000 year2 year2*year2 / D=p solution ddfm=kenwardroger;
random site;
random _resid_;
run;

*Relates to FIGURE 3b - annual population size estimates. Year fitted as a categorical variable*;

proc glimmix  maxopt=200 pconv=1.11022E-6 data=work.kitrend plots=(PearsonPanel);
class site year2;
model AOT2 = lnS2000 year2 / D=p solution ddfm=kenwardroger;
random site;
random _resid_;
lsmeans year2 / ilink cl plot=mean(ilink);
run;

*TESTING FOR GEOGRAPHICAL VARIATION IN BLACK-LEGGED KITTIWAKE POPULATION TRENDS*;

proc glimmix  maxopt=200 pconv=1.11022E-6 data=work.kitrend plots=(PearsonPanel);
class site;
model AOT2 = lnS2000 lat2 year2 year2*year2 lat2*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*Relates to FIGURE S2 - overall population trend for each colony type (solid lines). Year fitted as a covariate*;

proc glimmix maxopt=200 pconv=1.11022E-5 data=work.kitrend plots=(PearsonPanel);
class type site;
model AOT2 = lnS2000 lat2 type year2 year2*year2 type*year2 type*year2*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*Relates to FIGURE S2 - annual population size estimates for each colony type. Year fitted as a categorical variable*;

proc glimmix  data=work.kitrend plots=(PearsonPanel);
class type site year2;
model AOT2 = lnS2000 lat2 type year2 type*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
lsmeans type*year2 / ilink cl plot=mean(ilink);
run;


*ARCTIC TERN POPULATION TRENDS*;

*Relates to FIGURE 3C - overall population trend (solid line). Year fitted as a covariate*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.attrend plots=(PearsonPanel);
class site;
model pairs2 = lnS2000 year2 year2*year2 / D=p solution ddfm=kenwardroger;
random site;
random _resid_;
run;

*Relates to FIGURE 3C - annual population size estimates. Year fitted as a categorical variable*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.attrend plots=(PearsonPanel);
class site year2;
model pairs2 = lnS2000 year2 / D=p solution ddfm=kenwardroger;
random site;
random _resid_;
lsmeans year2 / ilink cl plot=mean(ilink);
run;

*TESTING FOR GEOGRAPHICAL VARIATION IN ARCTIC TERN POPULATION TRENDS*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.attrend plots=(PearsonPanel);
class site;
model pairs2 = lnS2000 lat2 year2 year2*year2 lat2*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*Relates to FIGURE S2 - overall population trend for each colony type (solid lines). Year fitted as a covariate*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.attrend plots=(PearsonPanel);
class type site;
model pairs2 = lnS2000 lat2 type year2 year2*year2 type*year2  / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;

*Relates to FIGURE S2 - annual population size estimates for each colony type. Year fitted as a categorical variable*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.attrend plots=(PearsonPanel);
class type site year2;
model pairs2 = lnS2000 lat2 type year2 type*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
lsmeans type*year2 / ilink cl plot=mean(ilink);
run;


*COMMON GUILLEMOT POPULATION TRENDS*;

*Relates to FIGURE 3d - overall population trend (solid line). Year fitted as a covariate*;

proc glimmix maxopt=200 pconv=1.11022E-7 data=work.gutrend plots=(PearsonPanel);
class site;
model AOT = lnS2000 year2 year2*year2 / D=p solution ddfm=kenwardroger;
random site;
random _resid_;
run;

*Relates to FIGURE 3d - annual population size estimates. Year fitted as a categorical variable*;

proc glimmix maxopt=200 pconv=1.11022E-3 data=work.gutrend plots=(PearsonPanel);
class site year2;
model AOT = lnS2000 year2 / D=p solution ddfm=kenwardroger;
random site;
random _resid_;
lsmeans year2 / ilink cl plot=mean(ilink);
run;

*TESTING FOR GEOGRAPHICAL VARIATION IN COMMON GUILLEMOT POPULATION TRENDS*;

proc glimmix maxopt=200 pconv=1.11022E-6 data=work.gutrend plots=(PearsonPanel);
class site;
model AOT = lnS2000 lat2 year2 year2*year2 lat2*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;


*ATLANTIC PUFFIN POPULATION TRENDS*;

*Relates to FIGURE 3e - overall population trend (solid line). Year fitted as a covariate*;

proc glimmix maxopt=200 pconv=1.11022E-5 data=work.putrend plots=(PearsonPanel);
class site;
model AOT = lnS2000 year2 year2*year2 / D=p solution ddfm=kenwardroger;
random site;
random _resid_;
run;

*Relates to FIGURE 3e - annual population size estimates. Year fitted as a categorical variable*;

proc glimmix maxopt=200 pconv=1.11022E-5 data=work.putrend plots=(PearsonPanel);
class site year2;
model AOT = lnS2000 year2 / D=p solution ddfm=kenwardroger;
random site;
random _resid_;
lsmeans year2 / ilink cl plot=mean(ilink);
run;

*TESTING FOR GEOGRAPHICAL VARIATION IN ATLANTIC PUFFIN POPULATION TRENDS*;

proc glimmix maxopt=200 pconv=1.11022E-5 data=work.putrend plots=(PearsonPanel);
class site;
model AOT = lnS2000 lat2 year2 year2*year2 lat2*year2 / D=p SOLUTION DDFM=KENWARDROGER;
random site;
random _resid_;
run;