Ajuste do modelo de von Bertalanffy

O modelo de crescimento de von Bertalanffy é muito utilizado para descrever a variação de comprimento de peixes, moluscos e crustáceos ao longo do tempo.

A seguir apresento um passo-a-passo para ajustar a curva aos dados de comprimento (Lt) na idade (t), analisar os parâmetros e fazer o gráfico.
O ajuste é feito de forma não linear pela função nls,  os intervalos de confiança das estimativas são calculados com a função confint e o coeficiente de determinação (R2) pela função Rsq do pacote qpcR. As funções expression e substitute são utilizadas para escrever as equações no gráfico.

t	Lt
1	102,0
2	167,0
3	219,4
4	260,7
5	294,9
6	323,2
7	343,0
8	369,5
9	401,7
10	410,0

 

# carrega pacote para cálculo do R2

library("qpcR")

# importa dados da área de transferência

dat.tL <- read.delim("clipboard",dec=",")

attach(dat.tL)

# ajuste do modelo

vb.pargo <- nls(Lt~Linf*(1-exp(-k*(t-t0))),start=list(Linf=500,k=0.2,t0=0))
summary(vb.pargo)

Formula: Lt ~ Linf * (1 - exp(-k * (t - t0)))

Parameters:
      Estimate Std. Error t value Pr(>|t|)   
Linf 501.51567   19.81444  25.311 3.84e-08 ***
k      0.16185    0.01541  10.504 1.55e-05 ***
t0    -0.46264    0.13937  -3.319   0.0128 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.587 on 7 degrees of freedom

Number of iterations to convergence: 5
Achieved convergence tolerance: 5.881e-07

# calcula intervalo de confiança dos parâmetros

confint(vb.pargo)

            2.5%       97.5%
Linf 461.8889677 562.3393409
k      0.1250699   0.1996879
t0    -0.8411141  -0.1621223


1-(deviance(vb.pargo)/((length(Lt)-1)*var(Lt))) # R2 "na mão"

[1] 0.9976615

Rsq(vb.pargo) # R2 pelo pacote qpcR

[1] 0.9976615

Rsq.ad(vb.pargo) # R2 ajustado pelo pacote qpcR
[1] 0.9969934 

# desenha o gráfico. Em Windows substituir ["\U221E"] por [infinity]

plot(Lt~t,xlab="idade (anos)",ylab="comprimento total (mm)",
xlim=range(0,10),ylim=range(0,500),cex.lab=1.2,
main=expression(L[i]==L["\U221E"]*"["*1-e^{-k(t-t[0])}*"]"),cex.main=1.5)
 

#desenha a curva

curve(coef(vb.pargo)[1]*(1-exp(-coef(vb.pargo)[2]*(x-coef(vb.pargo)[3]))),add=T,col="tomato1")
 

# coloca a legenda, deve-se clicar no gráfico para indicar o local da legenda

legend(locator(1),bty="n",legend=substitute(L[i]==Linf%*%"["*1-e^{-k%*%(t-t0)}*"]", list(Linf=round(coef(vb.pargo)[1],1),k=round(coef(vb.pargo)[2],2),t0=-round(coef(vb.pargo)[3],2))),cex=1.5)

detach(dat.tL)