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)