Assignment 5 — Solution

> slugs<-read.table('http://www.bio.ic.ac.uk/research/mjcraw/statcomp/data/slugsurvey.txt', header=TRUE)

The models that were fit in class

Poisson Models

#Common Mean Model

poi.1<-function(data,p) -sum(log(dpois(data$slugs, lambda=p)))
> nlm(function(p) poi.1(slugs,p),2)->out1

#Separate Means Model

poi.2<-function(data,p)
{
field.dummy<-as.numeric(data$field)-1
mylambda<-p[1]+p[2]*field.dummy
negloglike<- -sum(log(dpois(data$slugs,lambda=mylambda)))
negloglike
}
> nlm(function(p) poi.2(slugs,p),c(1.2,1),hessian=TRUE)->out2

ZIP Models

#Common lambda and theta

zip1<-function(data,p)
{
lambda<-p[1]
theta<-p[2]
zero.term<-sum(log(theta+(1-theta)*
dpois(data$slugs[data$slugs==0],lambda)))
other.term<-sum(log((1- theta)*dpois(data$slugs[data$slugs>0], lambda)))
negloglike<- -(zero.term+other.term)
negloglike
}

> nlm(function(p) zip1(slugs,p),c(2.2,.4))->out7

#Different lambda, common theta

zip2<-function(data,p)
{
field.dummy<-as.numeric(slugs$field)-1
mylambda<-p[1]+p[3]*field.dummy
theta<-p[2]
zero.term<-sum(ifelse(data$slugs==0,log(theta+(1-theta)*
dpois(data$slugs,lambda=mylambda)),0))
other.term<-sum(ifelse(data$slugs>0,log((1-theta)*
dpois(data$slugs,lambda=mylambda)),0))
negloglike<- -(zero.term+other.term)
negloglike
}

> nlm(function(p) zip2(slugs,p),c(3,.42,-.4))->out8

#Common lambda, different theta

zip3<-function(data,p)
{
field.dummy<-as.numeric(data$field)-1
mylambda<-p[1]
theta<-p[2]+p[3]*field.dummy
zero.term<-sum(ifelse(data$slugs==0,log(theta+(1-theta)*
dpois(data$slugs,lambda=mylambda)),0))
other.term<-sum(ifelse(data$slugs>0,log((1-theta)*
dpois(data$slugs,lambda=mylambda)),0))
negloglike<- -(zero.term+other.term)
negloglike
}

> nlm(function(p) zip3(slugs,p),c(3,.62,-.4))->out9

Log-normal models

#Common mean and variance

norm.neglike<-function(data,p)
{
t.y<-log(data$slugs+1)
negloglike<- -sum(log(dnorm(t.y,mean=p[1],sd=p[2])))
negloglike
}

> mean(log(slugs$slugs+1))
[1] 0.733247
> sd(log(slugs$slugs+1))
[1] 0.7414184

> nlm(function(p) norm.neglike(slugs,p),c(.73,.74))->out.norm

norm.like<-function(data,out)
{
t.y<-log(data$slugs+1)
negloglike<- -sum(log(dnorm(t.y,mean=out$estimate[1],
sd=out$estimate[2])*1/(data$slugs+1)))
out<-list(negloglike,out$estimate)
names(out)<-c("minimum","estimate")
out
}

> norm.like(slugs,out.norm)->out11

#Log-transform model with different means

norm.neglike2<-function(data,p)
{
t.y<-log(data$slugs+1)
field.dummy<-as.numeric(data$field)-1
mu<-p[1]+field.dummy*p[3]
negloglike<- -sum(log(dnorm(t.y,mean=mu,sd=p[2])))
negloglike
}

> tapply(log(slugs$slugs+1),slugs$field,mean)
Nursery Rookery
0.4967358 0.9697582

> nlm(function(p) norm.neglike2(slugs,p),c(.5,74,.5))->outnorm2

#new likelihood function

norm.like2<-function(data,out)
{
t.y<-log(data$slugs+1)
field.dummy<-as.numeric(data$field)-1
mu<-out$estimate[1]+field.dummy*out$estimate[3]
negloglike<- -sum(log(dnorm(t.y,mean=mu,
sd=out$estimate[2])*1/(data$slugs+1)))
out<-list(negloglike,out$estimate)
names(out)<-c("minimum","estimate")
out
}

> norm.like2(slugs,outnorm2)->out12

The new models for this assignment

Negative Binomial Models

#Common mean and theta

nbin1.nursery<-function(data,p)
{
field.dummy<-as.numeric(data$field)-1
mu<-p[1]
theta<-p[2]
negloglike<- -sum(log(dnbinom(data$slugs,mu=mu,size=theta)))
negloglike
}

> nlm(function(p) nbin1.nursery(slugs,p),c(1,.2))->out3

#check convergence
> out3$code
[1] 1
> out3$gradient
[1] -3.961435e-05 3.689138e-05

#Separate means, same theta

nbin2.nursery<-function(data,p)
{
field.dummy<-as.numeric(data$field)-1
mu<-p[1]+p[3]*field.dummy
theta<-p[2]
negloglike<- -sum(log(dnbinom(data$slugs,mu=mu,size=theta)))
negloglike
}

> nlm(function(p) nbin2.nursery(slugs,p),c(1,.2,.2))->out4

#check convergence
> out4$code
[1] 1
> out4$gradient
[1] -4.168520e-06 -1.136868e-06 -1.676881e-06

#Same mean, separate thetas

nbin3.nursery<-function(data,p)
{
field.dummy<-as.numeric(data$field)-1
mu<-p[1]
theta<-p[2]+p[3]*field.dummy
negloglike<- -sum(log(dnbinom(data$slugs,mu=mu,size=theta)))
negloglike
}

> nlm(function(p) nbin3.nursery(slugs,p),c(1,.2,.2))->out5

#check convergence
> out5$code
[1] 1
> out5$gradient
[1] -6.971642e-06 3.979039e-06 -1.328574e-06

#Separate means and thetas

nbin4.nursery<-function(data,p)
{
field.dummy<-as.numeric(data$field)-1
mu<-p[1]+p[3]*field.dummy
theta<-p[2]+p[4]*field.dummy
negloglike<- -sum(log(dnbinom(data$slugs,mu=mu,size=theta)))
negloglike
}

> out5$estimate
[1] 2.0913779 0.2506013 1.6472337
> out4$estimate
[1] 1.2749991 0.7859308 1.0000000
> nlm(function(p) nbin4.nursery(slugs,p),c(1.2,.25,.8,1))->out6

#check convergence
> out6a$code
[1] 1
> out6a$gradient
[1] 1.174757e-05 2.600586e-05 -3.370815e-05 -3.768655e-05

New ZIP Model

#Different lambda and theta

zip4<-function(data,p)
{
field.dummy<-as.numeric(data$field)-1
mylambda<-p[1]+p[3]*field.dummy
theta<-p[2]+p[4]*field.dummy
zero.term<-sum(ifelse(data$slugs==0,log(theta+(1-theta)* dpois(data$slugs,lambda=mylambda)),0))
other.term<-sum(ifelse(data$slugs>0,log((1-theta)* dpois(data$slugs,lambda=mylambda)),0))
negloglike<- -(zero.term+other.term)
negloglike
}

> table(slugs$slugs,slugs$field)/40

 Nursery Rookery
 0 0.625 0.225
 1 0.125 0.225
 2 0.050 0.200
 3 0.050 0.125
 4 0.050 0.050
 5 0.025 0.100
 6 0.025 0.025
 7 0.025 0.000
 8 0.000 0.025
 9 0.000 0.025
10 0.025 0.000

#use elements from first row for initial estimate
> nlm(function(p) zip4(slugs,p),c(3,.62,-.2,-.4))->out10

#check convergence
> out10$code
[1] 1
> out10$gradient
[1] 7.385899e-07 -1.040235e-05 1.762146e-06 -1.054445e-05

Lognormal models

#Same mean, different variances

norm.neglike3<-function(data,p)
{
t.y<-log(data$slugs+1)
field.dummy<-as.numeric(data$field)-1
mu<-p[1]
sd<-p[2]+field.dummy*p[3]
negloglike<- -sum(log(dnorm(t.y,mean=mu,sd=sd)))
negloglike
}

> tapply(log(slugs$slugs+1),slugs$field,sd)
  Nursery Rookery
0.7345233 0.6776645

> nlm(function(p) norm.neglike3(slugs,p),c(.73,.73,-.05))->outnorm3

norm.like3<-function(data,out)
{
t.y<-log(data$slugs+1)
field.dummy<-as.numeric(data$field)-1
mu<-out$estimate[1]
sd<-out$estimate[2]+field.dummy*out$estimate[3]
negloglike<- -sum(log(dnorm(t.y,mean=mu, sd=sd)*1/(data$slugs+1)))
out<-list(negloglike,out$estimate)
names(out)<-c("minimum","estimate")
out
}

> norm.like3(slugs,outnorm3)->out13

#Different means and different variances

norm.neglike4<-function(data,p)
{
t.y<-log(data$slugs+1)
field.dummy<-as.numeric(data$field)-1
mu<-p[1]+field.dummy*p[3]
sd<-p[2]+field.dummy*p[4]
negloglike<- -sum(log(dnorm(t.y,mean=mu,sd=sd)))
negloglike
}

> nlm(function(p) norm.neglike4(slugs,p),c(.5,.73,.5,-.05))->outnorm4

#check convergence
> outnorm4$code
[1] 1
> outnorm4$gradient
[1] 2.032152e-06 -1.480771e-05 -2.600586e-06 1.841727e-05
> outnorm4$estimate
[1] 0.49673535 0.72528297 0.47302233 -0.05614327

norm.like4<-function(data,out)
{
t.y<-log(data$slugs+1)
field.dummy<-as.numeric(data$field)-1
mu<-out$estimate[1]+field.dummy*out$estimate[3]
sd<-out$estimate[2]+field.dummy*out$estimate[4]
negloglike<- -sum(log(dnorm(t.y,mean=mu, sd=sd)*1/(data$slugs+1)))
out<-list(negloglike,out$estimate)
names(out)<-c("minimum","estimate")
out
}

> norm.like4(slugs,outnorm4)->out14

Square root transformed models

We need to add 1 to the response to account for zero values, not because the square root function would be undefined, but because our back-transformation of the likelihood would be undefined. In the back transformation we end up dividing by the square root of the response.

norm.sqrt<-function(data,p)
{
t.y<-sqrt(data$slugs+1)
mu<-p[1]
my.sd<-p[2]
negloglike<- -sum(log(dnorm(t.y,mean=mu,sd=my.sd)) )
negloglike
}

#obtain initial estimates of mean and standard deviation
mean(sqrt(log(slugs$slugs+1)))
[1] 0.6365834
> sd(sqrt(log(slugs$slugs+1)))
[1] 0.5763338
> nlm(function(p) norm.sqrt(slugs,p),c(.63,.57))->out.sq

#check convergence
> out.sq$code
[1] 1
> out.sq$gradient
[1] 1.447808e-05 7.258905e-05
> out.sq$estimate
[1] 1.5488749 0.6131752

#Likelihood for square root transformation
norm.sq.like<-function(data,out)
{
t.y<-sqrt(data$slugs+1)
negloglike<- -sum(log(dnorm(t.y,mean=out$estimate[1], sd=out$estimate[2])*1/(2*sqrt(data$slugs+1))))
out<-list(negloglike,out$estimate)
names(out)<-c("minimum","estimate")
out
}

> norm.sq.like(slugs,out.sq)->out15

Obtain table of AIC values

model.names<-c('Pois.common', 'Pois.mean', 'NB.common', 'NB.mean', 'NB.theta', 'NB.mean.theta', 'Zip.common', 'Zip.mean', 'Zip.theta', 'Zip.mean.theta', 'lognormal', 'lognormal.mean', 'lognormal.sd', 'lognormal.meansd', 'sqrt model')

AIC.func<-function(model.list,n,modelnames)
{
output<-NULLfor (i in 1:length(model.list))
{
cur.model<-model.list[[i]]
LL<- -cur.model$minimum
K<-length(cur.model$estimate)
AIC<- -2*LL + 2*K
AICc<-AIC + 2*K*(K+1)/(n-K-1)
output<-rbind(output,c(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
}

> AIC.func(models,80,model.names)
         modelnames      LogL K      AIC     AICc      deltai    wi
1       Pois.common -176.8383 1 355.6766 355.7279 72.94402255 0.000
2         Pois.mean -171.1275 2 346.2551 346.4109 63.62702501 0.000
3         NB.common -144.3980 2 292.7961 292.9519 10.16805531 0.003
4           NB.mean -142.6750 3 291.3500 291.6658  8.88187640 0.006
5          NB.theta -138.2398 3 282.4796 282.7953  0.01147261 0.491
6     NB.mean.theta -137.1253 4 282.2505 282.7839  0.00000000 0.493
7        Zip.common -150.4711 2 304.9422 305.0980 22.31417029 0.000
8          Zip.mean -150.4209 3 306.8417 307.1575 24.37362949 0.000
9         Zip.theta -143.7118 3 293.4236 293.7394 10.95554452 0.002
10   Zip.mean.theta -143.2979 4 294.5958 295.1291 12.34522757 0.001
11        lognormal -147.7365 2 299.4730 299.6288 16.84492612 0.000
12   lognormal.mean -143.3864 3 292.7727 293.0885 10.30463606 0.003
13     lognormal.sd -147.6051 3 301.2101 301.5259 18.74204291 0.000
14 lognormal.meansd -143.2567 4 294.5133 295.0467 12.26280308 0.001
15       sqrt model -159.1684 2 322.3368 322.4927 39.70880290 0.000

From the table we would conclude that only two of the models, the negative binomial common mean and separate theta model and the negative binomial separate mean and separate theta model, are reasonable models for these data. The Akaike weights are essentially equal for these two and approximately zero for the rest. The interpretation of the Akaike weights is that under repeated sampling we'd expect half the time one model and the other half the time the other model would be ranked best by AIC_c.

Goodness of Fit

If you adhere to the AIC_c results religiously then you would choose the negative binomial model with separate mean and thetas as the best model. Alternatively since the two best models are essentially equivalent, it might make sense to choose the simpler of the two, the common mean and separate theta negative binomial model. I carry out a goodness of fit test for both, first using the randomization test and then again using the parametric version of the Pearson chi-squared test.

Goodness of Fit for best model

Randomization Test

First I obtain the expected frequencies under this model. I combine the last category with the remaining probabilities in the tail. I do the calculation separately for each field type using the estimated parameters from the model. For field type nursey the mean and size parameters are p[1] and p[2]. For the rookery slugs the mean and size parameters are p[1]+p[3] and p[2]+p[4] respectively.

> exp.nurse<-c(dnbinom(0:9, mu=out6$estimate[1], size=out6$estimate[2]), 1-pnbinom(9, mu=out6$estimate[1], size=out6$estimate[2]))
> sum(exp.nurse)
[1] 1
> exp.rook<-c(dnbinom(0:9, mu=out6$estimate[1]+out6$estimate[3], size=out6$estimate[2]+out6$estimate[4]), 1-pnbinom(9, mu=out6$estimate[1]+out6$estimate[3], size=out6$estimate[2]+out6$estimate[4]))
> sum(exp.rook)
[1] 1

To obtain a common probability estimate for both field types I calculate a weighted average of the individual field probabilities using their proportional sample sizes as the weights.

> model6.probs<-(40*exp.nurse+40*exp.rook)/(40+40)
> chisq.test(table(slugs$slugs),p=model6.probs,simulate.p.value=TRUE)

Chi-squared test for given probabilities with simulated p-value (based
on 2000 replicates)

data: table(slugs$slugs)
X-squared = 1.8488, df = NA, p-value = 0.9995

The randomization test fails to find evidence of lack of fit.

Parametric Test

First I obtain the expected frequencies using the expected probabilities calculated above.

> model6.freqs<-40*exp.nurse+40*exp.rook
> model6.freqs
[1] 33.5611500 15.0203648 10.3714943 7.0915699 4.7374935 3.1142344 2.0301618
[8]  1.3213344  0.8635899  0.5695965 1.3190106

Next I determine at what point to merge the tail frequencies.
> sum(model6.freqs[8:11])
[1] 4.073531

We should merge from position 7 forward. I also merge counts in postion 5 and 6 so that their combined total exceeds 5.
> Ei.6<-c(model6.freqs[1:4],sum(model6.freqs[5:6]),sum(model6.freqs[7:11]))
> obs<-table(slugs$slugs)

I carry out exactly the same merges on the observed counts and then carry out the Pearson test.
> Oi.6<-c(obs[1:4],sum(obs[5:6]),sum(obs[7:11]))
> Ei.6
[1] 33.561150 15.020365 10.371494 7.091570 7.851728 6.103693
> Oi.6
 0  1  2 3
34 14 10 7 9 6
> pearson.6<-sum((Oi.6-Ei.6)^2/Ei.6)
> 1-pchisq(pearson.6,length(Oi.6)-1-length(out6$estimate))
[1] 0.6106476

So we fail to reject the null hypothesis of no significant lack of fit. How many degrees of freedom did we have?
> length(Oi.6)-1-length(out6$estimate)
[1] 1

Clearly we're close to overfitting these data. We have only one degree of freedom left. If we had included another parameter in our model we would have not been able to carry out this test.

Goodness of Fit for the second-best but more parsimonious model

Randomization test

I obtain the expected probabilities for rookery and nursery slugs separately under this model and then combine the results by averaging. In both I combine the last category with the remaining probabilities in the tail. I do the calculation separately for each field type using the estimated parameters from the model. For field type nursey the mean and size parameters are p[1] and p[2]. For the rookery slugs the mean and size parameters are p[1] and p[2]+p[3] respectively.

> exp.nurse5<-c(dnbinom(0:9,mu=out5$estimate[1],size=out5$estimate[2]),
1-pnbinom(9,mu=out5$estimate[1],size=out5$estimate[2]))
> sum(exp.nurse)
[1] 1
> exp.rook5<-c(dnbinom(0:9,mu=out5$estimate[1],size=out5$estimate[2]+out5$estimate[3]),
1-pnbinom(9,mu=out5$estimate[1],size=out5$estimate[2]+out5$estimate[3]))
> sum(exp.rook)
[1] 1

> model5.probs<-(40*exp.nurse5+40*exp.rook5)/(40+40)
> chisq.test(table(slugs$slugs),p=model5.probs,simulate.p.value=TRUE)

Chi-squared test for given probabilities with simulated p-value (based
on 2000 replicates)

data: table(slugs$slugs)
X-squared = 2.769, df = NA, p-value = 0.9875

The randomization test fails to find evidence of lack of fit.

Parametric Test

> model5.freqs<-40*exp.nurse5+40*exp.rook5
> model5.freqs
[1] 32.6139148 14.8305972 10.2366468 6.9407027 4.6156646 3.0496809 2.0264146
[8]  1.3681602  0.9463188  0.6742496 2.6976499

I examine how many of the tail probabilities I need to combine.
> sum(model5.freqs[8:11] )
[1] 5.686378

The expected count for X = 4 is 4.6 and is very close to 5. Since 5 is only a guideline and up to 20% of the counts can be under 5, I elect to leave this cell alone to preserve degrees of freedom. I do merge cells 6 and 7, X=5 and X=6, however.
> Ei.5<-c(model5.freqs[1:5],sum(model5.freqs[6:7]),sum(model5.freqs[8:11]))

I repeat the merging for the observed counts.
> obs<-table(slugs$slugs)
> Oi.5<-c(obs[1:5],sum(obs[6:7]),sum(obs[8:11]))
> Ei.5
[1] 32.613915 14.830597 10.236647 6.940703 4.615665 5.076095 5.686378
> Oi.5
 0  1  2 3 4
34 14 10 7 4 7 4

The observed and expected cell counts look very close. I carry out the parametric test.

> pearson.5<-sum((Oi.5-Ei.5)^2/Ei.5)
> 1-pchisq(pearson.5,length(Oi.5)-1-length(out5$estimate))
[1] 0.7001923

So we have no evidence of a significant lack of fit. This test is a bit more powerful than was the test of the more complicated model in that we have three degrees of freedom.
> length(Oi.5)-1-length(out5$estimate)
[1] 3

Course Home Page