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

Fri Jul 22 19:17:32 CEST 2011

> summary(n1)

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

> summary(n2)

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 ...

> AICctab(n1,n2,nobs=nrow(subset(X,time>0)))
   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:
> 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:
>> 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)
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>