Gráfico de círculos

Elaborei um gráfico para a representação do número relativo de observações por categoria de determinadas variáveis registradas em diferentes pontos de coleta. No exemplo abaixo, eu tenho o número de embarcações observadas por classe de comprimento total nos municípios de São Paulo (dados são fictícios).

Nesta rotina a relativização é feita por ponto de coleta (município).

Dados (fictícios)

Município	0 ˫ 6	6 ˫ 9	9 ˫ 12	12 ˫ 15	15 ˫ 18	≥ 18
Ubatuba	37	46	7	7	0	0
Caraguatatuba	41	41	15	2	0	0
Ilhabela	50	39	9	0	0	0
São Sebastião	43	40	14	1	0	0
Bertioga	0	52	28	18	0	0
Santos/Guarujá	6	27	16	16	4	28
São Vicente	75	12	12	0	0	0
Praia Grande	60	39	0	0	0	0
Mongaguá	16	83	0	0	0	0
Itanhaém	0	85	14	0	0	0
Peruíbe	31	43	25	0	0	0
Iguape	38	60	0	0	0	0
Ilha Comprida	66	34	0	0	0	0
Cananéia	34	40	10	11	2	0
Total	497	641	150	55	6	28

# importa dados da área de transferência

dat.graf <- t(read.delim("clipboard",dec=",",row.names=1))
dat.graf

# indica o valor da variáveis que serão utilizadas

NCATY<-nrow(dat.graf)
NCATX<-ncol(dat.graf)
MARX<-10 # margem da abcissa
MARY<-5 # margem da ordenada
CIRC<-13 # tamanho máximo do círculo
COR<-"blue" # cor do círculo

# plota o gráfico

par(mar=c(MARX,MARY,2,2))
plot(c(1:NCATX),rep(0,NCATX),xlab="",ylab="",
xaxp=c(1,NCATX,NCATX-1),yaxp=c(1,NCATY+1,NCATY),
ylim=c(0,NCATY+1),xaxt="n",yaxt="n",type="n")
axis(1,at=c(1:NCATX),labels=colnames(dat.graf),las=2,cex.axis=1.5)
axis(2,at=c(1:NCATY),labels=rownames(dat.graf),las=2,cex.axis=1.5)
for (x in 1:NCATX) {
  for(y in 1:NCATY) {
    points(x,y,col=COR,pch=19,cex=dat.graf[y,x]*CIRC/sum(dat.graf[,x]))
  }
}
# caso a última linha da tabela seja uma totalização e deva ficar em destaque.
abline(v=NCATX-0.5,lty=2)

-------------

Algumas pessoas tiveram problemas com a visualização, compreensão e inserção do símbolos matemáticos de indicação dos intervalos de classe.

6 ˫ 9, 6 |- 9, [6,9[ ou [6,9) indica um intervalo de classe que inclui o limite inferior (6) e exclui o superior (9)

No Linux Ubuntu a forma correta de indicar o símbolo "|-" (˫) é Shift+Ctrl+U 02EB e o código para ">=" (≥) é Shift+Ctrl+U 2265, , como podemos ver na figura.

Dependendo no navegador e do sistema operacional estes símbolos podem não ser visualizados corretamente no blog.

Caso os símbolos matemáticos e as acentuações não sejam corretamente importados pelo R através área de transferência (clipboard) podemos gravar a tabela de dados e importa-la como comando read.csv.

Se mesmo assim os probelmas continuarem podemos re-escrever os nomes de colunas e linhas com colnames e rownames:

colnames(dat.graf)
colnames(dat.graf)[4]<-"São Sebastião"
colnames(dat.graf)[6]<-"Santos/Guarujá"

rownames(dat.graf)

rownames(dat.graf)<-c("0 ˫ 6","6 ˫ 9","9 ˫ 12","12 ˫ 15","15 ˫ 18","≥ 18")

Avaliação do direcionamento de pescarias

Em uma pescaria é importante determinar que espécies são alvo e quais podem ser consideradas fauna acompanhante. A análise da variação do direcionamento, ou da importância relativa de uma espécie na composição das capturas, ao longo do tempo pode indicar padrões na dinâmica das frotas pesqueiras ou na abundância do recurso. O trabalho  Biseau, A. 1998. Definition os a directed fishing effort in a mixed-species trawl fishery, and its impacts on stock assessments. Aquat. Living Resour. 11(3):119-136 apresenta um método simples e objetivo para esta avaliação.

Abaixo apresento um script para a aplicação do método a partir de dados de captura total por viagem (Ti) e da captura da espécie por viagem (Tis).

# simulação de dados para rodar o script

Ti <- rnorm(1000,4000,300)

Tis <- Ti*runif(1000,0,1)

dat.biseau<-data.frame(Ti,Tis)

rm(Ti,Tis)

# aplicação do método de Biseau

# Ti é a captura total na vigem e Tis a captura da espécie na viagem

summary(dat.biseau)
dat.biseau$C<-dat.biseau$Tis/dat.biseau$Ti
dat.biseau$j<-as.integer(dat.biseau$C*100)
agg.j<-aggregate(dat.biseau$Tis,list(j=dat.biseau$j),FUN=sum)
names(agg.j)<-c("j","TC")
dat.direc <- data.frame(0:100,rep(0,101))
names(dat.direc)<-c("j","TC")
for (i in 0:100) {
ifelse(nrow(subset(agg.j,j==i,TC))==0,
dat.direc$TC[dat.direc$j==i]<-0,
dat.direc$TC[dat.direc$j==i]<-agg.j$TC[agg.j$j==i])
}
dat.direc$P <- cumsum(dat.direc$TC)/sum(dat.direc$TC)*100
dat.direc
plot(dat.direc$P~dat.direc$j,type="l",xlab="proporção da espécie na descarga da viagem (%)",ylab="descargas acumuladas (%)",xlim=c(0,100),ylim=c(0,100))

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)  

Utilização do RKWard

O RKWard é uma órima interface para o R. É a que mais utilizo.

Além tornar algumas funções básicas acessíveis através de menu (algo que nunca utilizei), possui um ótimo editor de script e outras funcionalidades. Por exemplo o editor de script indica os parâmetros de cada função e verifica os parêntesis. Também torna-se mais fácil exportar gráficos para diversos formatos.

O RKWard está nos repositórios de diversas distribuições Linux e pode ser facilmente instalado por um gerenciado de pacotes como o Synaptic.

No entanto, dependendo da versão do R instalada, o RKWard pode passar a indicar ao final de cada linha o aviso "Erro: unprotect(): somente 1 itens protegidos" ou outras mensagens erros. Isto ocorre se a versão do R não for compatível com a do RKWard.

O RKWard está em inglês e roda no Linux. Em seu site é indicada uma forma de utiliza-lo no Windows.

Para resolver este problema a versão do RKWard deve ser atualizada pela versão mais nova disponível. Detalhes pode ser obtidos em http://rkward.sourceforge.net/

Por exemplo, no Ubuntu, pode-se adicionar o repositório indicado no PPA em https://launchpad.net/~rkward-devel/+archive/rkward-stable-cran (versão estável) ou, em último caso, em https://launchpad.net/~rkward-devel/+archive/rkward-devel-cran (versão em desenvolvimento).

Análise de Variância e teste Kruskal–Wallis

É frequente termos a necessidade de verificar se as diferenças de uma variável medida em diversos grupos são significativas. Por exemplo, podemos querer testar a diferença dos comprimentos de peixes capturados em diferentes locais  ou da potência do motor de embarcações de diferentes grupos identificados em um cluster.

A seguir apresento uma sequencia básica para a realização da análise, incluindo testes à posteriori paramétricos e não paramétricos. Para a comparação múltipla não paramétrica é utilizado o pacote npmc. Os dados utilizados no exemplo são gerados no início da rotina.

# gera dados

GRP <- as.factor(rep(c("A","B","C","D"), c(30,50,40,20)))
UA <- round(rnorm(30,360,75),1)
UB <- round(rnorm(50,500,90),1)
UC <- round(rnorm(40,150,30),1)
UD <- round(rnorm(20,385,60),1)
U <- c(UA,UB,UC,UD)
dat.GRP <- data.frame(GRP,U)
rm(GRP,UA,UB,UC,UD,U)

# explora dados
attach(dat.GRP)
dat.GRP
summary(dat.GRP)
aggregate(U,by = list(GRP=GRP), FUN = "length")
aggregate(U,by = list(GRP=GRP), FUN = "min")
aggregate(U,by = list(GRP=GRP), FUN = "max")
aggregate(U,by = list(GRP=GRP), FUN = "mean")
aggregate(U,by = list(GRP=GRP), FUN = "sd")
aggregate(U,by = list(GRP=GRP), FUN = "median")
boxplot(U~GRP,notch=T,at=rank(tapply(U,GRP, median)))

# --- adaptado de Dalgaard 2008 Introductory Statistics with R
xbar <- tapply(U, GRP, mean)
s <- tapply(U, GRP, sd)
n <- tapply(U, GRP, length)
sem <- s/sqrt(n)
stripchart(U~GRP, method="jitter", jitter=0.05, pch=16, vert=T)
arrows(1:length(levels(GRP)),xbar+sem,1:length(levels(GRP)),xbar-sem,angle=90,code=3,length=.1)
lines(1:length(levels(GRP)),xbar,pch=4,type="b",cex=2)
rm(xbar,s,n,sem)
# ---

# ANOVA e teste de comparação múltipla de Tukey
aov.GRP<-aov(U~GRP)
summary(aov.GRP)
TukeyHSD(aov.GRP, ordered = TRUE)
plot(TukeyHSD(aov.GRP, ordered = TRUE))

# teste de Kruskal–Wallis e teste de comparação múltipla não paramétrica
kru.GRP<-kruskal.test(U~GRP)
kru.GRP

ori.names<-names(dat.GRP)
names(dat.GRP)<-c("class","var") # o pacote npmc requer que o nome da coluna que identifica o grupo seja class e a de valores var

library(npmc)
npmc(dat.GRP)

names(dat.GRP)<-ori.names
rm(ori.names)
detach(dat.GRP)

Gráfico de dipersão com boxplots

O gráfico de dispersão e uma forma bastante ilustrativa de representar a correlação entre duas variáveis. Informações sobre a distribuição das variáveis da ordenada e da abscissa podem ser representadas por boxplots colocados na margem do gráfico de dispersão.

Para o exemplo abaixo utilizei os dados e as curvas de regressão do tópico Relação comprimento-peso de peixes

# Plotagem do gráfico de dispersão

zones=matrix(c(2,0,1,3),ncol=2,byrow=TRUE)
layout(zones,widths=c(20,5),heights=c(5,20))
par(mar=c(8,8,0,0))
plot(P~C,subset=G=="M",pch=20,xlab="Ct (mm)",ylab="Pt (g)",cex.axis=1.3,cex.lab=1.3)

# para plotagem da regressão linear
curve(exp(coef(lm.CPm)[1])*x^(coef(lm.CPm)[2]),col="blue",add=T)
posicao<-locator(1) # clique no gráfico para indicar onde colocar a legenda
legend(posicao,bty="n",cex=1.5,text.col="blue",legend=substitute(Pt==a%*%Ct^b,list(a=exp(coef(lm.CPm)[1]),b=round(coef(lm.CPm)[2],3))))
posicao$y <- posicao$y-500
legend(posicao,bty="n",cex=1.5,text.col="blue",legend=substitute(R^2==r,list(r=round(summary(lm.CPm)$adj.r.squared,3))))

# para plotagem da regressão não linear
curve(coef(nls.CPm)[1]*x^coef(nls.CPm)[2], col="darkgreen",add=T)
posicao<-locator(1) # clique no gráfico para indicar onde colocar a legenda
legend(posicao,bty="n",cex=1.5,text.col="darkgreen",legend=substitute(Pt==a%*%Ct^b,list(a=coef(nls.CPm)[1],b=round(coef(nls.CPm)[2],3))))
posicao$y <- posicao$y-500
legend(posicao,bty="n",cex=1.5,text.col="darkgreen",legend=substitute(R^2==r,list(r=round(Rsq.ad(nls.CPm),3))))

# para plotagem dos boxplots
par(mar=c(0,8,0,0))
boxplot(C,subset=G=="M",axes=FALSE, space=0, horizontal=TRUE)
par(mar=c(8,0,0,0))
boxplot(P,subset=G=="M",axes=FALSE, space=0)