Assignment 6 — Solution

> isles<-read.table( 'http://www.unc.edu/courses/2006spring/ecol/145/001/data/lab6/islands.txt', header=TRUE, sep=',')

Fitting the models

#Gleason Model

> gleason<-lm(sp.richness~log(island.area),data=isles)
> coef(gleason)
(Intercept) log(island.area)
  -171.0272         109.4431

The estimated model is μ = –171.03 + 109.44 × log(area), where μ is the mean of species richness.

#log-Arrhenius Model

> log.arrhen<-lm(log(sp.richness)~log(island.area),data=isles)
> coef(log.arrhen)
(Intercept) log(island.area)
  3.8389781       0.3283347

The estimated model is μ = 3.84 + 0.33 × log(area), where μ is the mean of log(species richness). Alternatively we can interpret μ as the geometric mean of species richness (see lecture 23).

#Arrhenius model

I exponentiate the estimates obtained from the log-Arrhenius model for use as initial estimates to the nls function for the Arrhenius model.

> arrhen<-nls(sp.richness~b0*island.area^b1, data=isles, start=list(b0=exp(coef(log.arrhen)[1]), b1=coef(log.arrhen)[2]))
> arrhen
Nonlinear regression model
model: sp.richness ~ b0 * island.area^b1
data: isles
        b0        b1
83.4789648 0.2549076
residual sum-of-squares: 451062.6

The estimated model is μ = 83.48 × area^0.26, where μ is the mean of species richness.

#Poisson model

> poisson.glim<-glm(sp.richness~log(island.area),data=isles,family=poisson)
> coef(poisson.glim)
(Intercept) log(island.area)
  4.1925024        0.2835591

#Negative binomial model

> library(MASS)
> negbin.glim<-glm.nb(sp.richness~log(island.area),data=isles)
> coef(negbin.glim)
(Intercept) log(island.area)
  3.9729777        0.3172363

Graphing the Models

The log-Arrhenius model, since it models the log of species richness, stands alone and I plot it by itself. The Gleason model and Arrhenius model are also graphically incompatible. But the Poisson model and negative binomial model are of the same form. They can be plotted together with the Gleason model, or together with the Arrhenius model, as I show below.

#Log-Arrhenius model

> plot(log(isles$island.area), log(isles$sp.richness), xlab='log(Area)', ylab='log(Species)')
> abline(log.arrhen,col=2,lty=2)
> mtext('log-Arrhenius model', side=3, line=.5)

#Arrhenius model

> arrhenius.func<-function(x) coef(arrhen)[1]*x^(coef(arrhen)[2])
> plot(isles$island.area,isles$sp.richness,xlab='Area',ylab='Species')
> lines(seq(0,25000,100),arrhenius.func(seq(0,25000,100)),lty=2,col=2)
> mtext('Arrhenius model',side=3,line=.5)

Fig. 1   Mean log(richness) and mean richness predicted by the log-Arrhenius and Arrhenius models, respectively

#Gleason, Poisson, and negative binomial models

> plot(log(isles$island.area), isles$sp.richness, xlab='log(Area)', ylab='Species')
> abline(gleason, col=2, lty=2, lwd=2)
> generic.func<-function(coef,x) exp(coef[1]+coef[2]*x)
> lines(seq(1,10,.1), generic.func(coef(poisson.glim), seq(1,10,.1)), lty=3, col=3, lwd=2)
> lines(seq(1,10,.1), generic.func(coef(negbin.glim), seq(1,10,.1)), lty=4, col=4, lwd=2)
> legend(1.1, 1200, c('Gleason', 'Poisson', 'negative binomial'), col=c(2,3,4), lty=c(2,3,4), bty='n',lwd=c(3,3,3))

Fig. 2   Mean richness as a function of log area as predicted by the Gleason, Poisson, and negative binomial models

Comparing the Models

The loglikelihood for the Arrhenius model needs to be calculated separately since the response is log-transformed in this model.

> norm.loglike<-function(data,model)
{
t.y<-log(data$sp.richness)
sigma2<-(sum(residuals(model)^2))/dim(data)[1]
loglike<-sum(log(dnorm(t.y, mean=predict(model), sd=sqrt(sigma2))*1/(data$sp.richness)))
out<-list(loglike, c(coef(model), sigma2))
out
}
> norm.loglike(isles,log.arrhen)->logarrhen.vals

The following function calculates the various statistics we need for model comparison.

> AIC.func<-function(LL,K,n,modelnames)
{
#LL is loglikelihood,
#K is number of estimated parameters
#n is the sample size
AIC<- -2*LL + 2*K
AICc<-AIC + 2*K*(K+1)/(n-K-1)
output<-cbind(LL,K,AIC,AICc)
colnames(output)<-c('LogL','K','AIC','AICc')
minAICc<-min(output[,"AICc"])
deltai<-output[,"AICc"]-minAICc
rel.like<-exp(-deltai/2)
wi<-round(rel.like/sum(rel.like),3)
out<-data.frame(modelnames,output,deltai,wi)
out
}

I create a list of model names, string together the loglikelihoods, and create a vector listing the number of estimated parameters in each model. Three parameters are estimated in each model (for maximum likelihood), except the Poisson model in which only two parameters were estimated.

> model.names<-c('Gleason', 'Arrhenius', 'Log-Arrhenius', 'Poisson', 'Neg Binom')
> loglike <-c(logLik(gleason), logLik(arrhen), logarrhen.vals[[1]], logLik(poisson.glim), logLik(negbin.glim))
> numparms<-c(3,3,3,2,3)

I carry out the model comparisons.

> AIC.func(loglike, numparms, dim(isles)[1], model.names)
     modelnames      LogL K       AIC      AICc     deltai    wi
1       Gleason -141.4406 3  288.8813  290.2146  16.586909 0.000
2     Arrhenius -140.4282 3  286.8563  288.1896  14.561930 0.000
3 Log-Arrhenius -133.1472 3  272.2944  273.6277   0.000000 0.632
4       Poisson -542.1834 2 1088.3668 1088.9984 815.370651 0.000
5     Neg Binom -133.6894 3  273.3788  274.7121   1.084434 0.367

From the model summary only the log-Arrhenius and negative binomial models are reasonable models for these data. Weight of evidence makes the log-Arrhenius almost twice as likely to be the best model as the negative binomial, but in truth the two models are very close (being within 1 of each other on the AIC scale).

It's worth reflecting on the fact that in Fig. 2, the graphs of the worst model (the Poisson) and one of the best models (negative binomial) are virtually identical. In a least squares sense these two models would be found to be nearly indistinguishable, but of course we didn't use least squares to fit these models. We used maximum likelihood. In maximum likelihood we're not evaluating fit by how close the mean trajectory comes to the observations. Instead we're trying to maximum the probability of obtaining the observed data. Thus when we compare models using loglikelihood or AIC we're evaluating probability generating mechanisms. From the loglikelihood (or the AIC) reported in the above table we would conclude that it is very unlikely that a Poisson model generated these data. It is far more likely for a negative binomial model to have generated these data.

Course Home Page