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

Dicas para Barplot

Apresento aqui mais algumas dicas para a barplots. Sobre barplots veja também esta postagem.

Tabela exemplo:

	G1	G2	G3
Local C	3	30	70
Local B	10	80	30
Local A	50	15	5

# gráfico de barras simples

dados<- read.delim("clipboard",row.names=1)
barplot(as.matrix(dados), ylim=c(0,140),
xlab="grupos",ylab="nº de viagens",
legend.text=row.names(dados),

args.legend=list(x = "topleft", bty="n"))
box()

para transpor a matriz:  t(as.matrix(dados))

para obter a frequência relativa: prop.table(as.matrix(dados),2)

para barras justapostas: beside=T

para barras horizontais: horiz=T

adiciona uma linha abaixo das barras (Opção ao box): axis.lty=1

Cálculo das fases da lua no R

A influência do ciclo lunar sobre os organismos marinhos é grande e, como consequência, sobre a atividade pesqueira. Muitas vezes ao analisar dados queremos relacionar a abundância das espécies ou a CPUE das pescarias à fase lunar.

Por esta razão, necessitei de uma função em R que indicasse a fase da lua a partir de uma data. Achei em http://www.paulsadowski.com/wsh/moonphase.htm um código em Visual Basic para este cálculo e o portei para R.
As fases da lua ficaram como:
nova -> crescente concava -> quarto crescente -> crescente convexa -> cheia -> minguante convexa -> quarto minguante -> minguante concava -> nova.

Verifiquei as respostas dadas pela função com o programa LunaBar e em um site Islâmico (http://islam.com.pt/), na seção Fases da Lua. Descobri que o calendário Islâmico é baseado no ciclo lunar e por isso eles oferecem a informação de maneira exata. Os resultados que obtive com a função foram bem próximos aos calculados no site e no LunaBar.

Acredito que a função pode ser utilizada sem inconvenientes, no entanto não posso garantir a exatidão dos cálculos.
A função pode ser facilmente alterada caso se deseje apenas a idade ou a fase da lua.
Achei interessante a utilização da função recode, do pacote car, que deve estar carregado.

require(car)

lua(2011,8,31)

Idade da lua (dias): 2
Fase da lua: crescente concava
Distância (raio): 56.67
Latitude Elíptica (graus): -4.8
Longitude Elíptica (graus): 190.11 

lua<- function(Y,M,D) {
# http://www.paulsadowski.com/wsh/moonphase.htm

# necessita do pacote cars
P2<-2*3.14159
YY<-Y-as.integer((12-M)/10)
MM=M+9
if (MM>=12) MM<-MM-12
K1=as.integer(365.25*(YY+4712))
K2=as.integer(30.6*MM+.5)
K3=as.integer(as.integer((YY/100)+49)*.75)-38

# J é a data às 12h UT do dia em questão

J<-K1+K2+D+59
if(J>2299160) J<-J-K3
# Calcula a fase sinódica
V<-(J-2451550.1)/29.530588853
V<-V-as.integer(V)
if(V<0) V<-V+1
IP<-V
# Idade da Lua em dias
AG<-IP*29.53
IP<-IP*P2 # Converte fase em radianos
# Calcula a distância
V<-(J-2451562.2)/27.55454988
V<-V-as.integer(V)
if(V<0) V<-V+1
DP<-V
DP<-DP*P2 # Converte em radianos

DI<-60.4-3.3*cos(DP)-.6*cos(2*IP-DP)-.5*cos(2*IP)
# Calcula a Latitude
V<-(J-2451565.2)/27.212220817
V<-V-as.integer(V)
if(V<0) V<-V+1
NP<-V
NP<-NP*P2 # Converte em radianos
LA<-5.1*sin(NP)
# Calcula a Longitude
V<-(J-2451555.8)/27.321582241
# Normaliza valores para o intervalo de 0 a 1
V<-V-as.integer(V)
if(V<0) V<-V+1
RP<-V
LO<-360*RP+6.3*sin(DP)+1.3*sin(2*IP-DP)+.7*sin(2*IP)
# fases em inglês
# http://home.hiwaay.net/~krcool/Astro/moon/moonphase/
Phase<-c("nova","crescente concava","quarto crescente","crescente convexa","cheia","minguante convexa","quarto minguante","minguante concava")
ThisPhase<-recode(as.integer(AG),"c(0,29)=1;c(1,2,3,4,5,6)=2;c(7)=3;c(8,9,10,11,12,13)=4;c(14)=5;c(15,16,17,18,19,20,21)=6;c(22)=7;c(23,24,25,26,27,28)=8")
message("Idade da lua (dias): ", as.integer(AG))
message("Fase da lua: ",Phase[ThisPhase])
message("Distância (raio): ",round(DI,2))
message("Latitude Elíptica (graus): ",round(LA,2))
message("Longitude Elíptica (graus): ",round(LO,2))
}

Preparação de uma matriz de dados biológicos

Ao prepararmos uma matriz de dados biológicos para análise multivariada temos que ter inicialmente dois cuidados: devemos fazer com que o identificador dos objetos (usualmente estações de coleta) sejam os nomes das linhas e devemos substituir NAs (células vazias) por zeros.

Normalmente os dados de abundância são submetidos a alguma transformação monotônica, como log(x+1), para tornar a distribuição normal, estabilizar a variância e fazer com que as medidas de distância trabalhem melhor.
Para a mudança dos nomes das linhas utilizamos a função rownames,  para substituição dos NAs por zeros is.na e, finalmente, para logaritimização log1p.

A seguir veremos um exemplo destas etapas iniciais de uma análise multivariada.

#dados

ST	SP1	SP2	SP3
ST1	4	2
ST2	8	4	1
ST3	1	3	5
ST4		3	7

# lê os dados
dat.bio <-read.delim("clipboard",row.names=1)

dat.bio 

    SP1 SP2 SP3
ST1   4   2  NA
ST2   8   4   1
ST3   1   3   5
ST4  NA   3   7

# substitui NAs por 0
dat.bio[is.na(dat.bio)]<-0
dat.bio

    SP1 SP2 SP3
ST1   4   2   0
ST2   8   4   1
ST3   1   3   5
ST4   0   3   7

# logaritimização  ln(x+1)
dat.biolog <- log1p(dat.bio)
dat.biolog
          SP1      SP2       SP3
ST1 1.6094379 1.098612 0.0000000
ST2 2.1972246 1.609438 0.6931472
ST3 0.6931472 1.386294 1.7917595
ST4 0.0000000 1.386294 2.0794415

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)

Como ordenar as categorias de um boxplot por suas medianas

A técnica gráfica de boxplot é muito utilizada na análise exploratória de dados, etapa essencial em um trabalho.

No entanto o R, por padrão, coloca os fatores em ordem alfabética. Uma forma de reordenar estes fatores é através da função ordered, como descrito neste tópico.

Outra forma muito útil é através do argumento at e das funções rank e tapply, exemplificado a baixo.

categ <- rep(c("c","b","a"),c(10,10,10))

valor<-c(rnorm(10,5,2),rnorm(10,12,4),rnorm(10,8,3))

dados1<-data.frame(categ,valor)

# boxplot normal, não ordenado

boxplot(dados1$valor~dados1$categ)

# boxplot ordenado

boxplot(valor~categ,

at=rank(tapply(dados1$valor,dados1$categ, median)))

A forma de ordenamento acima é a mais simples, mas se houver valores de mediana iguais as caixas sairão sobrepostas. Uma forma alternativa de ordenamento é:


categ2 <- with(dados1, factor(categ,
levels=levels(categ)[order(tapply(valor,categ,median))]))
boxplot(valor~categ2,dados1)
rm(categ2)


PS: encontrei esta dica no seguinte link, onde outras possibilidades também são apresentadas.

Gráfico ordenado com médias e intervalo de confiança

Este gráfico mostra os valores de b (coeficiente angular) e seus intervalos de confiança calculados para diversas espécies. Para melhor vizualização as espécies são ordenadas pelo valor de b e diferentes grupos (peixes, moluscos e crustáceos) são identificados com símbolos especícicos.

SP	GRP	b	ICi	ICs
S. brasiliensis	P	3,000	2,887	3,114
U. brasiliensis	P	2,993	2,931	3,055
U. mystacea	P	2,990	2,948	3,032
L. gastrophysus	P	2,897	2,841	2,952
C. jamaicensis	P	2,874	2,846	2,901
M. ancylodon	P	3,101	3,064	3,137
M. americanus	P	2,809	2,785	2,833
M. furnieri	P	2,946	2,935	2,957
P. punctatus	P	3,050	3,030	3,070
L. sanpaulensis	M	1,784	1,685	1,886
O. vulgaris	M	2,689	2,613	2,766
O. vulgaris F	M	2,651	2,539	2,764
O. vulgaris M	M	2,731	2,637	2,827
F. brasiliensis	C	2,409	2,352	2,466
F. brasiliensis F	C	2,350	2,262	2,439
F. brasiliensis M	C	2,201	2,011	2,392
F. paulensis	C	2,395	2,338	2,453
F. paulensis F	C	2,248	2,160	2,336
F. paulensis M	C	2,470	2,224	2,716
X. kroyeri	C	2,301	2,244	2,358
X. kroyeri F	C	2,305	2,228	2,383
X. kroyeri M	C	2,149	2,039	2,260

# carregando e verificando dados

rm(list=ls())

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

summary(dat.bvar)

                 SP     GRP         b              ICi             ICs      
 C. jamaicensis   : 1   C:9   Min.   :1.784   Min.   :1.685   Min.   :1.886 
 F. brasiliensis  : 1   M:4   1st Qu.:2.316   1st Qu.:2.232   1st Qu.:2.404 
 F. brasiliensis F: 1   P:9   Median :2.670   Median :2.576   Median :2.765 
 F. brasiliensis M: 1         Mean   :2.607   Mean   :2.527   Mean   :2.686 
 F. paulensis     : 1         3rd Qu.:2.934   3rd Qu.:2.877   3rd Qu.:2.956 
 F. paulensis F   : 1         Max.   :3.101   Max.   :3.064   Max.   :3.137 
 (Other)          :16                           

# reorganizando a tabela de dados
order(dat.bvar$b) # linhas em ordem crescente de b 

[1] 10 22 16 18 20 21 15 17 14 19 12 11 13  7  5  4  8  3  2  1  9  6

rank(dat.bvar$b) # posição de cada linha quando ordenada por b
[1] 20 19 18 16 15 22 14 17 21  1 12 11 13  9  7  3  8  4 10  5  6  2

dat.bvar$RAN <- rank(dat.bvar$b)
levels(dat.bvar$SP) # níveis de SP
 [1] "C. jamaicensis"    "F. brasiliensis"   "F. brasiliensis F"
 [4] "F. brasiliensis M" "F. paulensis"      "F. paulensis F"  
 [7] "F. paulensis M"    "L. gastrophysus"   "L. sanpaulensis" 
[10] "M. americanus"     "M. ancylodon"      "M. furnieri"     
[13] "O. vulgaris"       "O. vulgaris F"     "O. vulgaris M"   
[16] "P. punctatus"      "S. brasiliensis"   "U. brasiliensis" 
[19] "U. mystacea"       "X. kroyeri"        "X. kroyeri F"    
[22] "X. kroyeri M"

dat.bvar$SP<-ordered(dat.bvar$SP, levels=dat.bvar$SP[order(dat.bvar$b)]) # reordena SP pela ordem de b
levels(dat.bvar$SP) # níveis de SP reordenados

[1] "L. sanpaulensis"   "X. kroyeri M"      "F. brasiliensis M"
[4] "F. paulensis F"    "X. kroyeri"        "X. kroyeri F"    
[7] "F. brasiliensis F" "F. paulensis"      "F. brasiliensis" 
[10] "F. paulensis M"    "O. vulgaris F"     "O. vulgaris"     
[13] "O. vulgaris M"     "M. americanus"     "C. jamaicensis"  
[16] "L. gastrophysus"   "M. furnieri"       "U. mystacea"     
[19] "U. brasiliensis"   "S. brasiliensis"   "P. punctatus"    
[22] "M. ancylodon"
ord.sp<-dat.bvar$SP[order(dat.bvar$b)]

# criando o gráfico

attach(dat.bvar)
par(mar=c(8,4,2,2))
plot(c(1:nrow(dat.bvar)),rep(3,nrow(dat.bvar)),ylab="valores de b",xlab="",ylim=c(min(ICi)*0.97,max(ICs)*1.03),type="n",xaxt="n")
axis(1,1:nrow(dat.bvar),labels=ord.sp,las=2)
for (i in 1:length(ord.sp)){
segments(i,ICi[SP==ord.sp[i]],i,ICs[SP==ord.sp[i]])
}
points(RAN[GRP=="P"],b[GRP=="P"],pch=19)
points(RAN[GRP=="M"],b[GRP=="M"],pch=17)
points(RAN[GRP=="C"],b[GRP=="C"],pch=15)
legend(locator(1),bty="n",c("Peixes","Moluscos","Crustáceos"),pch=c(19,17,15)) # clique no gráfico para indicar o local da legenda

detach(dat.bvar)