[R] nls - can't get published AICc and parameters

[Ben Bolker]

Roland Sookias <r.sookias <at> gmail.com> writes:

> 
> Hi
> 
> I'm trying to replicate Smith et al.'s
> (http://www.sciencemag.org/content/330/6008/1216.abstract) findings by
> fitting their Gompertz and logistic models to their data (given in
> their supplement). I'm doing this as I want to then apply the
> equations to my own data.
> 
> Try as a might, I can't quite replicate them. Any thoughts why are
> much appreciated. I've tried contacting the authors but I think
> they're away.
> 
> The equations as I've used them are:
> 
> log(mass)~log(K)-log(K/M)*exp(-a*t)
> 
> and
> 
> log(mass)~C*t^G
> 
> The data file I've been using is attached. It starts at the K-Pg
> boundary with a body size of 3.3kg, following their description in the
> supplement.
> 
> Their parameters and AICc values are given in their paper. I get
> something quite close for some of them (~0.245 G, their gamma, K
> estimates reasonable etc.), but in the logistic model C is more like 3
> than 1.5, and the AICc values differ by ~10 rather than ~1.

  

Here's my attempt (your attachment didn't come through;
I typed in the data from the supplement of Smith et al.).

 I don't get quite the same answers they do, either.

## X <- read.table("mammals.txt",header=TRUE)
## 
## X <- transform(X,time=65.5-age,
##               LBM=log(maxbodymass/3.3))


X <- structure(list(age = c(0.013, 0.9035, 2.703, 4.465, 8.47, 13.79, 
19.5, 25.715, 31.15, 35.55, 42.9, 52.2, 57.25, 60.2, 63.6, 70.6, 
105.5), maxbodymass = c(10000, 15000, 17450, 17450, 17450, 6568, 
5917, 15000, 15000, 5907, 4500, 700, 700, 54, 54, 3.3, 5), 
order = structure(c(6L, 
6L, 6L, 6L, 6L, 6L, 6L, 5L, 5L, 5L, 2L, 4L, 4L, 1L, 1L, 3L, 1L
), .Label = c("Condylarthra", "Dinocerata", "Multituberculata", 
"Pantodonta", "Perissodactyla", "Proboscidea"), class = "factor"), 
    time = c(65.487, 64.5965, 62.797, 61.035, 57.03, 51.71, 46, 
    39.785, 34.35, 29.95, 22.6, 13.3, 8.25, 5.3, 1.9, -5.09999999999999, 
    -40), LBM = c(8.01641790350375, 8.42188301161191, 8.57317245915814, 
    8.57317245915814, 8.57317245915814, 7.59604218265983, 7.49166237420426, 
    8.42188301161191, 8.42188301161191, 7.4899708988348, 7.21791020728598, 
    5.35715786657097, 5.35715786657097, 2.79506157809184, 2.79506157809184, 
    0, 0.415515443961666)), .Names = c("age", "maxbodymass", 
"order", "time", "LBM"), row.names = c(NA, -17L), class = "data.frame")

X0 <- subset(X,time>0,select=c(LBM,time))
X0 <- rbind(X0,c(log10(3.3),time=0))
n1 <- nls(LBM~c0*time^gamma,data=X0,
    start=list(c0=1.5,gamma=0.25))

n2 <- nls(LBM~logK-(logK-logM0)*exp(-alpha*time),data=X0,
    start=list(logK=log(13183),logM0=log(6.92),alpha=0.08))

coef(n2)
library(bbmle)
AICctab(n1,n2,nobs=nrow(subset(X,time>0)))
## 9.3-8.2 = 1.1.

par(bty="l",las=1)
with(X,plot(log10(maxbodymass)~time,
                           pch=16,cex=2,
                           axes=FALSE,
            xlim=c(-9.5,65),ylim=c(0,5)))
axis(side=1,at=seq(-9.5,60.5,by=10),
     labels=seq(75,5,by=-10))
abline(v=0,lty=3)
tvec <- -5:70
pred1 <- 1/log(10)*(log(3.3)+predict(n1,newdata=data.frame(time=tvec)))
pred2 <- 1/log(10)*(log(3.3)+predict(n2,newdata=data.frame(time=tvec)))
lines(tvec,pred1,lty=2)
lines(tvec,pred2,lty=1)