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, 31 de agosto de 2011

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

quarta-feira, 8 de junho de 2011

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


terça-feira, 7 de junho de 2011

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)

sexta-feira, 3 de junho de 2011

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.

sexta-feira, 13 de maio de 2011

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)

sexta-feira, 3 de dezembro de 2010

Cálculo da moda

O R trás funções para o cálculo da média (mean) e da mediana (median) mas não para moda. A sugestão postada na lista R-Help por Dr. Brian D. Ripley para distribuições discretas foi:
 
statmod <- function(x) {
  z <- table(as.vector(x))
  names(z)[z == max(z)]
}
 
n <- c(1,2,2,2,3,3)
statmod(n)
[1] "2"

quinta-feira, 28 de outubro de 2010

Como colocar vírgulas como separador decimal em gráficos

Muitas vezes em um gráfico a abscissa e/ou a ordenada apresentam valores com casas decimais. O R por padrão utiliza ponto decimal e não vírgula. Para colocar a vírgula decimal podemos sobrescrever ou escrever os valores dos eixos com as funções axis, pretty, chartr, as.character e format. O eixo (axis) 1 é o inferior, o 2 da direita, o 3 o superior e o 4 da esquerda. O número de casas decimais é indicada em nsmall.
Se colocarmos na função plot a indicação axes=F nenhum dos eixos é desenhado. Abaixo estão exemplos de como colocar vírgulas decimais.

x<-c(1,2,3,4,5)
y<-c(0.1,0.3,0.5,0.7,0.9)
plot(x,y)
axis(2,at=pretty(y),labels=chartr(".", ",", as.character(format(pretty(y),nsmall=1))))

plot(x/100,y,axes=F)
axis(1,at=pretty(x/100),labels=chartr(".", ",", as.character(format(pretty(x/100),nsmall=2))))
axis(2,at=pretty(y),labels=chartr(".", ",", as.character(format(pretty(y),nsmall=1))))
box()