Skip to content

nls - can't get published AICc and parameters

4 messages · Ben Bolker, Roland Sookias

Ben Bolker
# 
Formula: LBM ~ c0 * time^gamma

Parameters:
      Estimate Std. Error t value Pr(>|t|)
c0     2.58478    0.36767   7.030 8.93e-06 ***
gamma  0.29257    0.03781   7.738 3.21e-06 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

Residual standard error: 0.7124 on 13 degrees of freedom

Number of iterations to convergence: 5
Achieved convergence tolerance: 6.138e-06
Formula: LBM ~ logK - (logK - logM0) * exp(-alpha * time)

Parameters:
      Estimate Std. Error t value Pr(>|t|)
logK   8.41177    0.24971  33.686 2.98e-13 ***
logM0  1.55057    0.63230   2.452 0.030467 *
alpha  0.07176    0.01504   4.770 0.000456 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

Residual standard error: 0.5359 on 12 degrees of freedom

Number of iterations to convergence: 4
Achieved convergence tolerance: 6.094e-06

  Note that I am fitting logK and logM0, not K and M0 ...
dAICc df
n2 0.0   4
n1 5.9   3

  I wasn't quite sure what they meant by starting at the smallest body
mass.  It looked from their figure in the supplementary info that they
might have constrained the curve to start from (0,0) [on the log scale].
I definitely don't get exactly their answers -- on the other hand, the
points don't even look like they're in exactly the same place as in
their figure.  I would definitely keep trying to contact the authors.
On 11-07-27 08:43 AM, Roland Sookias wrote:
4 days later
Roland Sookias
# 
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.

Cheers.

Roland
Ben Bolker
# 
Roland Sookias <r.sookias <at> gmail.com> writes:
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)
Roland Sookias
# 
I also see you divided the body mass by the mass at time zero - why was that?

Thanks so much for trying to help.
On Wed, Jul 27, 2011 at 1:08 AM, Ben Bolker <bbolker at gmail.com> wrote: