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, 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.