• Variáveis aleatórias teóricas como modelos
• Parâmetros e estimativas
• Simulação de amostras de variáveis aleatórias
• Ajuste de modelos por otimização
• Diagnóstico de modelos

Vamos retomar o tutorial sobre distribuições discretas, onde simulamos uma distribuição agregada de plantas em uma área dividida em quadrículas (parcelas):

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]
x2.parc <- cut(x2,breaks=0:20, labels=1:20)
y2.parc <- cut(y2,breaks=0:20, labels=LETTERS[1: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 ¹⁾:

p.mod2 <- profile(mod2)
par(mfrow=c(1,2))
plot.profmle(p.mod2)
par(mfrow=c(1,1))

• Um roteiro de ajuste de modelos na página do Laboratório de Ecologia Teórica do IB-USP.
• Outro roteiro, bem mais resumido, da disciplina de introdução à modelagem com verossimilhança de Chad Brassil.
• Excelentes exercícios de simulação e ajustes de distribuições no site de apoio de Bolker (2008). ²⁾

Indique os tutoriais e partes da leitura básica que merecem mais atenção na discussão aqui.