Ferramentas do usuário

Ferramentas do site


Barra lateral

CMQ
Centro de Métodos Quantitativos


USP ESALQ
Depto. de Ciências Florestais
ESALQ
UNIVERSIDADE de SÃO PAULO
Av. Pádua Dias, 11
Caixa Postal 09
13418-900 - Piracicaba - SP
BRASIL
biometria:verossim:03b-model

Essa é uma revisão anterior do documento!


Modelos com Parâmetros Constantes

Tutoriais

Poisson e Binomial Negativa

Vamos retomar o tutorial sobre distribuições discretas, onde simulamos uma distribuição agregada de planats em uma parcela:

set.seed(42)
cx <- runif(20,0,20)
cy <- runif(20,0,20)
px <- rnorm(2000)
py <- rnorm(2000)
x1 <- cx+px
y1 <- cy+py
x2 <- x1[x1>0&x1<20&y1>0&y1<20]
y2 <- y1[x1>0&x1<20&y1>0&y1<20]

O numero de plantas por parcela é obtido com:

cont2 <- data.frame(table(x2.parc,y2.parc))$Freq

A avaliação visual feita no tutorial anterior indica que a variável binomial negativa é uma descrição mais adequada destes dados. Vamos ajustar estes dois modelos e fazer a comparação.

Primeiro definimos funções de log-verossimilhança negativa para cada modelo:

LL.pois <- function(lam){
  -sum(dpois(cont2,lambda=lam,log=T))
}

LL.nbin <- function(media,k){
  -sum(dnbinom(cont2,mu=media,size=k,log=T))
}

Em seguida minimizamos estas funções. Para isto carregue o pacote bbmle e use a função mle2. É preciso fornecer valores iniciais razoáveis, no argumento start:

library(bbmle) # basta uma vez por seção 
mod1 <- mle2(LL.pois,start=list(lam=mean(cont2)))
mod2 <- mle2(LL.nbin,start=(list(media=mean(cont2),k=0.1)))

Como esperado, a binomial negativa é um modelo muito mais plausível:

> logLik(mod1)
'log Lik.' -1786.502 (df=1)
> logLik(mod2)
'log Lik.' -1014.969 (df=2)

Você pode obter um resumo dos modelo com o comando summary:

summary(mod1)
summary(mod2)

E podemos fazer um gráfico dos valores previsto pelos dois modelos com:

##MLEs de cada modelo
cf1 <- coef(mod1)
cf2 <- coef(mod2)

##grafico
cont2.f <- factor(cont2, levels=0:max(cont2))
plot(table(cont2.f)/400, xlab="N de indivíduos na parcela", ylab="Proporção das parcelas")
points(x=0:max(cont2),y=dpois(0:max(cont2),lambda=cf1), type="b", col="blue", lty=2)
points(x=0:max(cont2),y=dnbinom(0:max(cont2),mu=cf2[1],size=cf2[2]), type="b", col="red", lty=2)

Por fim, avalie o perfil de verossimilhança dos dois parâmetros da binomial negativa 1):

p.mod2 <- profile(mod2)
par(mfrow=c(1,2))
plot.profmle(p.mod2)
par(mfrow=c(1,1))
1)
para isto você vai precisar da função plot.profmle
biometria/verossim/03b-model.1244446647.txt.gz · Última modificação: 2015/08/10 20:48 (edição externa)