O R é um programa livre multiplataforma para análises estatísticas que pode ser baixado em seu site ou adicionado na lista de repositórios de máquinas linux. Suas possibilidades de aplicação em diversas áreas são praticamente ilimitadas.
Neste blog postarei o resultado de minha experiência em sua utilização nas áreas de dinâmica de populações de peixes, ciência pesqueira e ecologia.
As postagens deste blog se destina, além de mim mesmo, a iniciantes no R e alunos da minha área de atuação.
Aprendi muito em livros e nas listas de discussão R-help e a R_STAT, mas ainda tenho muito pela frente. Agradeço desde já qualquer contribuição.

quarta-feira, 4 de abril de 2012

Gráfico de barras com erro padrão

São muitas as opções para se fazer um gráfico de barras com a indicação do erro padrão (ou outra medida de dispersão) no R. Com uma rápida busca na rede pode-se ter dezenas de sugestões. Eu resolvi contribuir com mais uma.
Meu objetivo era representar a variação da densidade de algumas espécies medida em períodos noturnos e diurnos.
A partir de três tabelas, uma contendo as médias, outra o erro padrão e uma terceira com o número de observações, utilizo os comandos barplot, points, arrows e text para gerar o gráfico desejado.
Os dados das tabelas são copiados para o R pelo comando read.delim (clique aqui para mais detalhes) e transformados em matrizes pelo comando as.matriz.

Médias


SP1 SP2 SP3 SP4
D 10 12 15 20
N 20 15 12 10

Erros padrões


SP1 SP2 SP3 SP4
D 1 2 3 4
N 3 2 1 5

Número de observações


SP1 SP2 SP3 SP4
D 30 40 35 10
N 20 50 20 27

# copia os dados das tabelas para o R
means <- as.matrix(read.delim("clipboard",row.names=1))
se <- as.matrix(read.delim("clipboard",row.names=1))
n <- as.matrix(read.delim("clipboard",row.names=1))
 

# calcula a amplitude da variação do erro
se.sup<-means+se
se.inf<-means-se


# desenha o gráfico

bp<-barplot(means,beside=T,ylim=c(0,max(se.sup*1.15)),
ylab="densidade",legend.text=c("Diurno","Noturno"),
args.legend=list(x = "topleft", bty="n"))
points(bp,means,pch=19)
arrows(bp,se.sup,bp,se.inf, code=3,angle=90,length=0.05)
text(bp,se.sup+1,n)


 




domingo, 25 de março de 2012

Amostrar linhas em um data.frame

Algumas vezes necessitamos obter um subconjunto de um conjunto de dados. Se pudermos estabelecer um critério a partir dos valores deste conjunto de dados podemos utilizar o comando subset. Para obter uma amostra aleatória das linhas pensei na seguinte solução:

# cria um conjuto de dados com duas colunas e 500 linhas
dados <- data.frame(rnorm(500,5,3),rnorm(500,10,6))
dim(dados)
summary(dados)

# faz a amostra de 100 linhas sem reposição
dados.amostra<-dados[sample(1:nrow(dados),100,replace=F),]
dim(dados.amostra)
summary(dados.amostra)


sábado, 3 de março de 2012

Identificando grupos em um dendrograma

As técnicas de agrupamento (cluster analysis) são muito úteis para organizar, ou classificar, dados observados em estruturas de fácil interpretação. Uma ótima referência é o livro Numerical Ecology with R (Boccard et al. 2011). O tuturial do pacote vegan (clique aqui) também trás explicações e exemplos interessantes.
Os dendrogramas são como mobiles, seus grupos são sempre os mesmos, mas podem mudar de posição. O grupo 1 não é necessariamente o primeiro à esquerda, nem o grupo 2 vem a seguir. Isto pode complicar na hora de interpretar os fatores relacionados na formação dos grupos. Um bom exemplo está na página 38 do tutorial do vegan. Temos a figura de um dendrograma seguida de um boxplot. Se não prestarmos atenção poderemos ser levados a crer, por exemplo, que o grupo 2, onde encontramos a maior mediada é formado pelos objetos 1, 2, 10, 5, 6 e 7. No entanto, este conjunto de objetos forma o grupo 1, como veremos a seguir, que tem a menor mediana. A interpretação do dendrograma poder ser facilitada se o número do grupo puder ser visualizado em conjunto. Vamos ao exemplo do tutorial

# carrega a biblioteca e os dados
library(vegan)
data(dune)

# calcula a matriz de distância, as ligações, plota o dendrograma e identifica os grupos
dis <- vegdist(dune)
cluc <- hclust(dis, "complete")
plot(cluc)
rect.hclust(cluc, 3)

# verifica que objetos foram classificados em cada grupo
grp<-cutree(cluc, 3)
grp

# substitui no dendrograma o nome do objeto pelo número do grupo
plot(cluc, labels = as.character(grp))

# desenha no dendrograma os retângulos de cada grupo e os numera
# agradeço a dica de Elias T. Krainski (lista R-Br)
plot(cluc)
r <- rect.hclust(cluc, 3)
text(cumsum(sapply(r,length)),
     rep(mean(tail(unique(cluc$hei),2)), length(r)),
     paste(unique(grp[cluc$ord])))

Agora fica mais fácil relacionar o dendrograma ao boxplot no tutorial do vegan.

quarta-feira, 14 de dezembro de 2011

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")

terça-feira, 25 de outubro de 2011

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))


segunda-feira, 3 de outubro de 2011

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


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)