I have a data set on sand martin population sizes along a stretch of river over 40 years. The river is split into sections and the number of birds per section was counted.

I have been trying to model changes in population sizes at these different sites in relation to sahel weather, north african temperature and english temperature among other factors but have run into a problem. There are obviously two main problems associated with my data, the first is autocorrelation of population sizes between years and the second is correlation between sites. I am not sure which is best to use and am unsure if I am able to combine the two into a single model. I'll talk you through my model and have posted the questions at the end.

Before I tried the autocorrelation structure I tried a normal model and there was a big cone-shaped pattern in the residuals. So I tried to create the model with autocorrelation, using gls in R. This is the model I came up with:

m1<-gls(pop~sahel.Index + Arrival.Date + NAO + Year + locaion + North.african.weather,
data=ndd,
method="REML",
correlation = corAR1(form=~Year|locaction))


From what I understand this means that population is correlated by year within location. The residuals outputs from this model look fine, I have used the normalised residuals to take into account the correlation structure. The residuals are completely normally distributed and have no pattern whatsoever in the residuals vs fitted. I would post them but can't because of my low "reputation."

So I then wanted to try and add spatial correlation to the model. I tried the different correlation structures: corLin, corGaus, corExp and corSpher. corExp had the lowest AIC so went with that.

My model for the spatial correlation thus looked like this:

mbestcor<-gls(pop~sahel.Index + Arrival.Date + NAO + Year + locaion + North.african.weather,
data=ndd,
method="ML",
correlation = corExp(form=~scYear|loc.ID))


And the residuals vs fitted and histogram look like this:

These plots look slightly odd, the residuals vs fitted is slightly cone-shaped and the histogram although normally distributed has quite a long right hand tail.

When I compare the AIC of these two models the spatial model is lower than the first:

 AIC(m2, bestm)
df      AIC
m2    23 2171.202
bestm 23 2157.828


But obviously the residual plots for the second one are also worse than the first.

I have a few questions: Which model should I go with?

Am I best basing this decision on the AIC or on the normality of residuals?

Ideally, I guess I would like to have a model that combines the two, but after lookign at Zuur et al. 2009 and online I can't seem to find a package that will let me do this unless I use a model builder.

Any help would be very much appreciated.

Thanks,

Thomas