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)  

Relação comprimeto-peso em peixes

A relação comprimento-peso (comprimento-massa) é utilizada para descrever a variação de peso em função da variação de comprimento, ou seja, para estimar o peso de um peixe a partir de um dado comprimento. Normalmente esta relação é descrita por um modelo de potência, P=a×C^b (veja Sparre e Venema, 1997 item 2.6).

A relação comprimento-peso também é utilizada para indicar a condição física ou higidez do organismo e seu padrão de crescimento como isométrico ou alométrico. Diz-se que o crescimento é isométrico quando o organismo cresçe na mesma taxa em todas as dimensões. Neste caso o valor de "b" será 3. A condição de isometria é um pressuposto para alguns métodos de avaliação de estoque.

A seguir posto uma sequencia para o ajuste dos parâmetros da relação comprimento-peso (a e b) pelos métodos linearizado e não linearizado, para a avaliação da significância da diferença entre os parâmetros ajustados para machos e fêmeas, e para a determinação da significância da diferença do parâmetro "b" do valor 3 (teste de isometria).

O juste linear é realizado com a função ln sobre os dados logaritmizados.  O ajuste não linear é realizado pela função nls. Os intervalos de confiança das estimativas são estimados pela função confint e o R2 do ajuste não linear pela função Rsq.ad do pacote qpcR. A isometria é verificada pelo teste t para a comparação de duas inclinações (Zar, Biostatistical Analysis). A comparação entre os modelos lineares ajustados é realizado através da ANCOVA (Faraway, 2002 pg 160). A comparação dos modelos não lineares é realizada por máxima verossimilhança, adaptado do método proposto por Kimura, 1980. Neste á utilizado a estatística qui-quadrado (pchisq)

No início da sequencia os dados de comprimento e peso de machos e fêmeas são gerados para termos um conjunto de dados para nos permitir seguir o exemplo.
Para fazer gráficos mais elaborados veja o tópico "Gráfico de Dispersão com Boxplots". 

# gera dados (C, P)
# =======================
# rotina para gerar um conjunto de dados
# Ma, Mb, Fa e Fb são os parâmetros do modelo de potência que
# que descrevem a relação comprimento-peso. MC, MP, FC e FP
# são os comprimentos e pesos de machos e fêmeas

Ma <- 2.45E-05
Fa <- 2.45E-05
Mb <- 3.0152
Fb <- 2.9164
MC <- round(rnorm(100,400,70),0)
FC <- round(rnorm(100,300,65),0)
MP <- round((Ma*MC^Mb)+((Ma*MC^Mb)*rnorm(100,0,0.08)),1)
FP <- round((Fa*FC^Fb)+((Fa*FC^Fb)*rnorm(100,0,0.08)),1) 

# prepara conjunto de dados
# =========================
# organiza os dados gerados em um conjunto (data.frame),
# indicando quais dados são de machos e quais são de fêmeas

dat.CPm <- data.frame(MC,MP)
names(dat.CPm)<-c("C","P")
dat.CPf <- data.frame(FC,FP)
names(dat.CPf)<-c("C","P")
dat.CP <- rbind(dat.CPm,dat.CPf)
dat.CP$G<-as.factor(rep(c("M","F"),each=100,1))

# visualiza dados
# ===============

attach(dat.CP)
summary(dat.CP)
boxplot(C~G,ylab="Lt mm")
plot(P~C,subset=G=="M")
plot(P~C,subset=G=="F")
plot(P~C,subset=G=="M",col="blue",ylim=c(0,6000))
points(P~C,subset=G=="F",col="red")

# ajusta curva de potência pelo método de linearização
# ====================================================
# a transformação logarítmica das variáveis P e C é
# aplicadas para a linearização na relação.

lm.CPm <- lm(log(P)~log(C),subset=G=="M")
summary(lm.CPm)
confint(lm.CPm)
lm.CPf <- lm(log(P)~log(C),subset=G=="F")
summary(lm.CPf)
confint(lm.CPf)

plot(log(P)~log(C),subset=G=="M",col="blue",ylim=c(log(35),log(6000)))
points(log(P)~log(C),subset=G=="F",col="red")
abline(lm.CPm,col="blue")
abline(lm.CPf,col="red")

# verifica isometria pelo teste t
# ===============================

tm<-(coef(summary(lm.CPm))[2,1]-3)/coef(summary(lm.CPm))[2,2]
dt(tm,nrow(dat.CP)-2)

tf<-(coef(summary(lm.CPf))[2,1]-3)/coef(summary(lm.CPf))[2,2]
dt(tf,nrow(dat.CP)-2)

# ANCOVA - análise de covariância
# ===============================
 
# testa inclinação
lm.b <- lm(log(P)~log(C)*G)
summary(lm.b)

# testa intercepto
lm.a <- lm(log(P)~log(C)+G)
summary(lm.a)

# ajusta curva de potência pelo método não linear
# ===============================================

require(qpcR)
nls.CPm<-nls(P~a*C^b,subset=G=="M",start=list(a=1E-05,b=3))
summary(nls.CPm)
confint(nls.CPm)
Rsq.ad(nls.CPm)

nls.CPf<-nls(P~a*C^b,subset=G=="F",start=list(a=1E-05,b=3))
summary(nls.CPf)
confint(nls.CPf)
Rsq.ad(nls.CPf)

plot(P~C,subset=G=="M",col="blue",ylim=c(0,6000))
points(P~C,subset=G=="F",col="red")
curve(coef(nls.CPm)[1]*x^coef(nls.CPm)[2], col="blue",add=T)
curve(coef(nls.CPf)[1]*x^coef(nls.CPf)[2], col="red",add=T)

# compara modelos por verossimilhança
# ===================================
# adaptado de Kimura 1980 Likelihood methods for the von Bertalanffy growth curve
 
# número médio de observações
par.NMed <- mean(c(nrow(subset(dat.CP,G=="M")),
nrow(subset(dat.CP,G=="F"))))

# cria objetos com os parâmetros das curvas de machos e fêmeas
par.CPf <- c(coef(nls.CPf),deviance(nls.CPf))
names(par.CPf)<- c("a","b","SSQ")

par.CPm <- c(coef(nls.CPm),deviance(nls.CPm))
names(par.CPm)<- c("a","b","SSQ")

# calcula a soma dos quadrados residual total
par.SSQ <- par.CPm[c("SSQ")]+ par.CPf[c("SSQ")]

par.CPf
par.CPm
par.SSQ

# ajuste do modelo com coef. linear (a) comum
nls.CPa <- nls(P ~ a*C^(ifelse(G=="M",bm,bf)),
start=list(a=1E-05,bm=3,bf=3))
par.CPa <- c(coef(nls.CPa),deviance(nls.CPa))
names(par.CPa)<- c("a","bm","bf","SSQ")

# ajuste do modelo com coef. angular (b) comum
nls.CPb <- nls(P ~ (ifelse(G=="M",am,af))*C^b,
start=list(am=1E-05,af=1E-05,b=3))
par.CPb <- c(coef(nls.CPb),deviance(nls.CPb))
names(par.CPb)<- c("am","af","b","SSQ")

# ajuste do modelo com ambos coef. comuns
nls.CPab <- nls(P ~ a*C^b,
start=list(a=1E-05,b=3))
par.CPab <- c(coef(nls.CPab),deviance(nls.CPab))
names(par.CPab)<- c("a","b","SSQ")

# lista o conjunto de parâmetros calculadas
par.NMed
par.CPf
par.CPm
par.SSQ
par.CPa
par.CPb
par.CPab

# teste Chi2 (Qui quadrado) para coef. linear
par.Chia <- -2*log((par.SSQ/par.CPa[c("SSQ")])^par.NMed)
names(par.Chia) <- "Chi2"
par.Chia # valor calculado de Chi2
1-pchisq(par.Chia, 1) # valor p da estatística Chi2

# teste Chi2 para coef. angular
par.Chib <- -2*log((par.SSQ/par.CPb[c("SSQ")])^par.NMed)
names(par.Chib) <- "Chi2"
par.Chib # valor calculado de Chi2
1-pchisq(par.Chib, 1) # valor p da estatística Chi2

# teste Chi2 para o modelo conjunto
par.Chiab <- -2*log((par.SSQ/par.CPab[c("SSQ")])^par.NMed)
names(par.Chiab) <- "Chi2"
par.Chib # valor calculado de Chi2
1-pchisq(par.Chiab, 2) # valor p da estatística Chi2

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)