Category: portuguese

Verossimilhança para um modelo de efeitos aleatórios

Modelos com termos de efeito aleatórios estão sendo cada vez mais usados na modelagem de dados. Nestes modelos, parte da variabilidade observada na resposta pode ser explicada por um conjunto de fatores de efeito fixo, outra parte pode estar associada a realização de variáveis latentes (ou não observáveis)  a qual é atribuída o efeito aleatório. […]

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