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

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