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)
 

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)