###################################################
### chunk number 1: 
###################################################
lista <- installed.packages()
lista[rownames(lista)=="vegan"|rownames(lista)=="bbmle"|
      rownames(lista)=="poilog",c(1,3)]


###################################################
### chunk number 2: 
###################################################
install.packages("vegan")
install.packages("bbmle")
install.packages("poilog")


###################################################
### chunk number 3: 
###################################################
library(vegan)
library(bbmle)
library(poilog)
library(MASS)


###################################################
### chunk number 4: 
###################################################
source("funcoes_logser.r")
source("funcoes_lognormal.r")
source("funcoes_niche_part.r")
source("rankplots.r")


###################################################
### chunk number 5: 
###################################################
##Fixa uma semente de n aleatorios
set.seed(51)
##Gera abundancias de especies com o modelo neutro de Caswell
com1 <- caswell(nu=5,lambda=0.1,tmax=1000)


###################################################
### chunk number 6: 
###################################################
radplot(com1,lwd=1.5)


###################################################
### chunk number 7: 
###################################################
## Alfa teorico
5/0.1
## Alfa obtido pelo ajuste a logserie
com1.alfa <- fit.logser(com1)$summary["alfa"]
com1.alfa


###################################################
### chunk number 8: 
###################################################
## total de individuos da comunidade
com1.N <- sum(com1)
com1.N
##Numero de especies da comunidade
com1.S <- length(com1)
com1.S
## Uma tabela com N de individuos por abundancia
com1.tab <- as.data.frame(table(com1))
names(com1.tab)[1] <- "Abund"
com1.tab$Abund <- as.integer(as.character(com1.tab$Abund))
##inspecionando os valores
com1.tab
## Probabilidades atribuidas pelo modelo a cada abundancia observada
com1.tab$p.ls <- dlr(com1.tab$Abund,N=com1.N,alfa=com1.alfa)
##inspecionando primeiros valores
head(com1.tab)


###################################################
### chunk number 9: 
###################################################
## Acrescentando o numero previsto de especies para cada abundancia
com1.tab$Pred.Ab.ls <- com1.tab$p.ls*com1.S
##inpecionando
head(com1.tab)


###################################################
### chunk number 10: 
###################################################
com1.LL.ls <- sum(log(com1.tab$p.ls)*com1.tab$Freq)


###################################################
### chunk number 11: 
###################################################
## Os mesmos procedimentos para lognormal, truncada em 0.5
## parametros da lognormal
com1.pln <- coef(lnormt.fit(com1))
com1.pln
## Probabilidades atribuidas a cada abundancia pela lognormal
com1.tab$p.ln <- plnormt.ab(com1.tab$Abund,mu=com1.pln[1],sigma=com1.pln[2])
## N previsto de especies em cada abundancia
com1.tab$Pred.Ab.ln <- com1.tab$p.ln*com1.S
##inspecionando primeiras linhas da tabela
head(com1.tab)
## Log-verossimilhanca do modelo log-normal
com1.LL.ln <- sum(log(com1.tab$p.ln)*com1.tab$Freq)


###################################################
### chunk number 12: 
###################################################
## comparando as duas verossimilhancas
com1.LL.ls
com1.LL.ln
## Razao de verossimilhancas
com1.LL.ls-com1.LL.ln
## Razao de verossimilhanças em escala aritmética
exp(com1.LL.ls-com1.LL.ln)


###################################################
### chunk number 13: 
###################################################
## AICS
##logserie
-2*com1.LL.ls+2
## lognormal
-2*com1.LL.ln+4


###################################################
### chunk number 14: 
###################################################
set.seed(42)
com2 <- caswell.ln(S0=200,n0=20,pbirth=0.105,pdeath=0.1,tmax=250)


###################################################
### chunk number 15: 
###################################################
radplot(com2)


###################################################
### chunk number 16: 
###################################################
##ajustes dos dois modelos
com2.fit.ls <- fit.logser(com2)
com2.fit.ln <- lnormt.fit(com2)


###################################################
### chunk number 17: 
###################################################
## Numero de individuos e especies
com2.N <- sum(com2)
com2.S <- length(com2)
com2.LL.bs <- sum(log(rad.bs(com2,com2.N,com2.S)/com2.S))
com2.AIC.bs <- -2*com2.LL.bs


###################################################
### chunk number 18: 
###################################################
## resultados
## serie logaritmica
com2.fit.ls$summary
## lognormal
c(N=com2.N,S=com2.S,coef(com2.fit.ln),
  loglik=logLik(com2.fit.ln),AIC=AIC(com2.fit.ln))
## brokenstick
c(N=com2.N,S=com2.S,coef(com2.fit.ln),
  loglik=com2.LL.bs,AIC=com2.AIC.bs)


###################################################
### chunk number 19: 
###################################################
##razoes de verossimilhanca log-normal x brokenstick
as.numeric(com2.LL.bs-logLik(com2.fit.ln))
## delta-AIC brokenstick x lognormal
AIC(com2.fit.ln)-com2.AIC.bs


###################################################
### chunk number 20: 
###################################################
## Um vetor que repete um código numérico para cada espécie
## o número de indivíduos que a espécie tem na comunidade
com1.ind <- rep(1:length(com1),com1)
## uma amostra de 10% da comunidade
set.seed(42)
com1.s <- table(sample(com1.ind,size=round(length(com1.ind)*0.1)))


###################################################
### chunk number 21: 
###################################################
w.plot(com1)
lines(rad.logser(com1))
points.w.plot(com1.s,col="red")


###################################################
### chunk number 22: 
###################################################
com1.fit.ls <- fit.logser(com1)
com1s.fit.ls <- fit.logser(com1.s)
com1.fit.ls$summary
com1s.fit.ls$summary


###################################################
### chunk number 23: 
###################################################
lines(rad.logser(com1.s),col="red",lwd=2)


###################################################
### chunk number 24: 
###################################################
lines(rad.logser(com1.s,alfa=com1s.fit.ls$summary["alfa"]),col="blue", lty=2)


###################################################
### chunk number 25: 
###################################################
## Um vetor que repete um código numérico para cada espécie
## o número de indivíduos que a espécie tem na comunidade
com2.ind <- rep(1:length(com2),com2)
## uma amostra de 10% da comunidade
set.seed(42)
com2.s <- table(sample(com2.ind,size=round(length(com2.ind)*0.1))) 


###################################################
### chunk number 26: 
###################################################
##ajuste dos dados da amostra
com2s.fit.ln <- lnormt.fit(com2.s)
## comparando as estimativas dos parametros
## Comunidade
coef(com2.fit.ln)
## Amostra
coef(com2s.fit.ln)


###################################################
### chunk number 27: 
###################################################
as.numeric(coef(com2s.fit.ln)[1])-log(0.1)


###################################################
### chunk number 28: 
###################################################
w.plot(com2.s)
lines(rad.lnormt(com2.s),col="blue")
lines(rad.lnormt(com2.s,sigma=coef(com2.fit.ln)[2]), col="red")
legend(x="topright",
       legend=c("Amostra",expression(paste(sigma, " comunidade"))),
       lty=1,col=c("blue","red"))