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.

quinta-feira, 22 de setembro de 2011

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

Nenhum comentário:

Postar um comentário