Estimativa da Mortalidade Total (Z) por Curvas de Captura - Composição de Comprimentos

Esta postagem está relacionada a Estimativa da Mortalidade Total (Z) por Curvas de Captura - Composição de Idades. 

A seguir está o procedimento para o cálculo de Z a partir da composição de comprimentos das capturas. Para aplica-lo é necessário, além dos dados de captura (C) por classe de comprimento (L), é necessário que indiquemos os parâmetros da curva de crescimento de von Bertalanffy (Linf, k e t0).

O comprimento indicado (L1) é o limite inferior da classe. Os parâmetros de crescimento a serem utilizados serão Linf = 460 mm, k = 0,198 / ano e t0 = -0,271 / ano.

L1	C
30	18
60	28
90	43
120	51
150	50
180	42
210	40
240	33
270	29
300	18
330	15
360	8
390	2
420	3
450	1

# importa e visualizadados
dat.CL <- read.delim("clipboard",dec=",")
dat.CL

# entra parâmetros da curva de crescimento
Linf <- 460
k <- 0.198
to <- -0.271

# calcula o intervalo de classe utilizado
IC <- dat.CL[2,1]-dat.CL[1,1]
IC

# estima a idade de quando o espécime entra na classe de comprimento
attach(dat.CL)
dat.CL$t <- to-(1/k)*log(1-L1/Linf)

# calcula o tempo que um espécime fica na classe de comprimento
dat.CL$dt <- (1/k)*log((Linf-L1)/(Linf-(L1+IC)))

# estima a idade média do espécime na classe
dat.CL$tm <- (to-(1/k)*log(1-L1/Linf)+to-(1/k)*log(1-(L1+IC)/Linf))/2

# estima o log da captura por unidade de tempo
dat.CL$cdt <- log(C/dat.CL$dt)
dat.CL
detach(dat.CL)

# desenha o gráfico e indica os pontos selecionados
attach(dat.CL)
par(mar=c(5,5,4,2))
plot(cdt~tm,xlab="idade (anos)",ylab=expression(log(C["L,L+1"] / Delta* t)),
cex.lab=1.5,pch=19,col="blue",xlim=c(0,max(na.omit(tm))*1.1))
# clique no gráfico para selecionar os pontos inicial e final.
# Depois de selecionar os dois pontos dê um clique direito para sair.
idn.CL<-identify(cdt~tm,plot=F)
idn.CL
dat.CLs <- dat.CL[c(idn.CL[1]:idn.CL[2]),]
dat.CLs
points(dat.CLs$tm,dat.CLs$cdt,cex=1.5,col="red")
detach(dat.CL)

# ajusta o modelo linear para o cálculo de Z
attach(dat.CLs)
lm.Zl <- lm(cdt~tm)
summary(lm.Zl)
lines(c(min(tm)-0.5,max(tm)+0.5),
  c(coef(lm.Zl)[1]+coef(lm.Zl)[2]*(min(tm)-0.5),coef(lm.Zl)[1]+coef(lm.Zl)[2]*(max(tm)+0.5)),
  col="red")
detach(dat.CLs)

# calcula Z, seus intervalos de confiança
# Z = -b; % de indivíduos que morrem anualmente = 1-exp(-Z)
Zl <- -coef(lm.Zl)[2]
Zl
-confint(lm.Zl)[2,]

# Mortalidade expressa em porcentagem
1-exp(-Zl)

# Porcentagem de sobreviventes (S)
exp(-Zl)

Estimativa da Mortalidade Total (Z) por Curvas de Captura - Composição de Idades

A mortalidade é um parâmetro importante para a compreensão da dinâmica de uma coorte. Através desta taxa pode-se estimar o número (absoluto ou relativo) de indivíduos que morrem na população em um dado período.

Existem diversos métodos para a estimativa da taxa instantânea de mortalidade total (Z). O detalhamento destes podem ser encontrados em referências como Chapman e Robson (1960), Pauly (1983, 1984a e 1984b), Pauly e Morgan (1987), Gulland e Rosenberg (1992) e Sparre e Venema (1997).

No R há o pacote FSA - Fisheries stock assessment methods, que trás várias funções para o estudo de dinâmica de populações e avaliação de estoques.

O cálculo da taxa de mortalidade total é simples e vale fazê-lo "na mão". Abaixo indico o procedimento para o cálculo de Z a partir de dados de captura (C) na idade (I). Z é uma taxa instantânea. No entanto, pode-se calcular os valores percentuais de mortos e sobreviventes (S) com facilidade.

  

I	C
0	1867,1
1	5090,5
2	5304,7
3	2186,6
4	947,5
5	380,7
6	176,8
7	75,3
8	31,6

# importa e visualiza dados
dat.CI <- read.delim("clipboard",dec=",")
dat.CI

# desenha o gráfico e indica os pontos selecionados
attach(dat.CI)
par(mar=c(5,5,4,2))
plot(log(C)~I,xlab="idade (anos)",ylab=expression(log(C["t,t+1"])),
cex.lab=1.5,pch=19,col="blue",ylim=c(0,max(log(C)*1.1)),xlim=c(0,max(I)*1.1))
# clique no gráfico para selecionar os pontos inicial e final.
# Depois de selecionar os dois pontos dê um clique direito para sair.
idn.CI<-identify(log(C)~I,plot=F)
idn.CI
dat.CIs <- dat.CI[c(idn.CI[1]:idn.CI[2]),]
dat.CIs
points(dat.CIs$I,log(dat.CIs$C),cex=1.5,col="red")
detach(dat.CI)

# ajusta o modelo linear para o cálculo de Z
attach(dat.CIs)
lm.Zi <- lm(log(C)~I)
summary(lm.Zi)
lines(c(min(I)-0.5,max(I)+0.5),
  c(coef(lm.Zi)[1]+coef(lm.Zi)[2]*(min(I)-0.5),coef(lm.Zi)[1]+coef(lm.Zi)[2]*(max(I)+0.5)),
  col="red")
detach(dat.CIs)

# calcula Z (Z = -b) e seus intervalos de confiança
Zi <- -coef(lm.Zi)[2]
Zi
- confint(lm.Zi)[2,]

# Mortalidade expressa em porcentagem
1-exp(-Zi)

# Porcentagem de sobreviventes (S)
exp(-Zi)  

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)