## Code for the analysis of model fitting and selection fo
## Pardini R, Bueno AdA, Gardner TA, Prado PI, Metzger JP, 2010 Beyond the Fragmentation Threshold Hypothesis: Regime Shifts in Biodiversity Across Fragmented Landscapes. PLoS ONE 5(10): e13666. doi:10.1371/journal.pone.0013666 

################################################################################
##DATA Structure:
## explanatory variables were in data frame 'dri', which has patches as rows,
## and variable 'area' for patch area (ha) and variable 'paisag'
## as a factor that identify each landscape.
## Richness and abundances are in the object "sem.didelphis"
## See datastruct.Rout
################################################################################

#########################################################################################
## 1. Models for species richness
## Response was modelled as a Poisson variate
#########################################################################################

## 1.1.
## Log-likelihood functions of the models that are non-standard combinations
## of different models for each landscape context.
## Models that are a single function for all landscapes
## were set directly through glm, see function seleciona.mod.c.id below

## Model F
LL.mod5.c.id=function(a1,a2,a3,b2){
	LL1=-sum(dpois(dep.var[dri$paisag=="tap"],lambda=a1,log=T))
	LL2=-sum(dpois(dep.var[dri$paisag=="cau"],lambda=a2+b2*log10(dri$area[dri$paisag=="cau"]),log=T))
	LL3=-sum(dpois(dep.var[dri$paisag=="rg"],lambda=a3,log=T))
	LL1+LL2+LL3
      }

## Model D with different slopes ofr each landscape
LL.mod6.c.id=function(a1,a2,a3,b1,b2,b3){
	LL1=-sum(dpois(dep.var[dri$paisag=="tap"],lambda=a1+b1*log10(dri$area[dri$paisag=="tap"]),log=T))
	LL2=-sum(dpois(dep.var[dri$paisag=="cau"],lambda=a2+b2*log10(dri$area[dri$paisag=="cau"]),log=T))
	LL3=-sum(dpois(dep.var[dri$paisag=="rg"],lambda=a3+b3*log10(dri$area[dri$paisag=="rg"]),log=T))
	LL1+LL2+LL3
      }

##Model E, with same slope for the landscapes with patch-area effect
LL.mod7.c.id=function(a1,a2,a3,b23,k){
	LL1=-sum(dpois(dep.var[dri$paisag=="tap"],lambda=a1,log=T))
	LL2=-sum(dpois(dep.var[dri$paisag=="cau"],lambda=a2+b23*log10(dri$area[dri$paisag=="cau"]),log=T))
	LL3=-sum(dpois(dep.var[dri$paisag=="rg"],lambda=a3+b23*log10(dri$area[dri$paisag=="rg"]),log=T))
	LL1+LL2+LL3
      }

##Model E, with different slopes for the landscapes with patch-area effect
LL.mod8.c.id=function(a1,a2,a3,b2,b3){
	LL1=-sum(dpois(dep.var[dri$paisag=="tap"],lambda=a1,log=T))
	LL2=-sum(dpois(dep.var[dri$paisag=="cau"],lambda=a2+b2*log10(dri$area[dri$paisag=="cau"]),log=T))
	LL3=-sum(dpois(dep.var[dri$paisag=="rg"],lambda=a3+b3*log10(dri$area[dri$paisag=="rg"]),log=T))
	LL1+LL2+LL3
      }

## 1.2.
## Function to fit the models by minimizing log-likelihood functions
## and return the Model selection Table fo AICc
## Models were fit by minimization of the log-likelihood functions above or by ordinary Poisson glms
seleciona.mod.c.id=function(){
	##coeficients starting values
        cf3=as.numeric(coef(glm(dep.var~dri$paisag/log10(dri$area)-1,family=poisson(link="identity"))))
        cf1 <- tapply(dep.var,dri$paisag,mean)
        ##model fitting with glm and mle2 optimization
        ## model A
	mod1= glm(dep.var~1,family=poisson(link="identity"))
        ## model B
        mod2= mle2(dep.var~dpois(lambda=a1+a2*log10(dri$area)),start=list(a1=cf3[2],a2=cf3[5]))
        ## model C
        mod3= glm(dep.var~dri$paisag,family=poisson(link="identity"))
        ## model D, same slopes
        mod4= mle2(dep.var~dpois(lambda=a1+a2*log10(dri$area)),parameters=list(a1~dri$paisag),start=list(a1=cf3[2],a2=cf3[5]))
        ## model F
        mod5= mle2(LL.mod5.c.id,start=list(a1=cf1[1],a2=cf3[2],a3=cf1[3],b2=cf3[5]))
        ## model D, different slopes
        mod6= mle2(LL.mod6.c.id,start=list(a1=cf3[1],a2=cf3[2],a3=cf3[3],b1=cf3[4],b2=cf3[5],b3=cf3[6]))
        ## model E, same slope
        mod7= mle2(LL.mod7.c.id,start=list(a1=cf1[1],a2=cf3[2],a3=cf3[3],b23=mean(cf3[5:6])))
        ## model E, different slopes
        mod8= mle2(LL.mod8.c.id,start=list(a1=cf1[2],a2=cf3[2],a3=cf3[3],b2=cf3[5],b3=cf3[6]))
        ##Results
	list(tabela=AICctab(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,sort=T,delta=T,weights=T,nobs=length(dep.var)),
		m1=mod1,m2=mod2,m3=mod3,m4=mod4,m5=mod5,m6=mod6,m7=mod7,m8=mod8)
      }

## 1.3.
## Model fiting and model selection table
## Loading required packages
library(bbmle)
libray(MASS)
## Assigning the desired response variable to the object dep.var
## and run the analysis
## Richness of endemics were mean number of species in two years
## so it was multiplied by a fixed offset = 2
## Endemic species richness
dep.var <- 2*sem.didelphis$rich.end
rich.end <- seleciona.mod.c.id()
## Non-endemic species richness
dep.var <- sem.didelphis$ab.n.end
ab.n.end <- seleciona.mod.c.id()



#########################################################################################
## 2. Models for abundances
## Response was modelled as a Poisson variate or Negative Binomial variate (see text)
#########################################################################################

## 2.1.
## Log-likelihood functions of the models that are non-standard combinations
## of different models for each landscape context.
## Models that are a single function for all landscapes
## were set directly through glm, see function seleciona.mod.id below

## Model F
LL.mod5.id=function(a1,a2,a3,b2,k){
	LL1=-sum(dnbinom(dep.var[dri$paisag=="tap"],mu=a1,size=k,log=T))
	LL2=-sum(dpois(dep.var[dri$paisag=="cau"],lambda=a2+b2*log10(dri$area[dri$paisag=="cau"]),log=T))
	LL3=-sum(dnbinom(dep.var[dri$paisag=="rg"],mu=a3,size=k,log=T))
	LL1+LL2+LL3
      }

## Model D with different slopes of each landscape
LL.mod6.id=function(a1,a2,a3,b1,b2,b3){
	LL1=-sum(dpois(dep.var[dri$paisag=="tap"],lambda=a1+b1*log10(dri$area[dri$paisag=="tap"]),log=T))
	LL2=-sum(dpois(dep.var[dri$paisag=="cau"],lambda=a2+b2*log10(dri$area[dri$paisag=="cau"]),log=T))
	LL3=-sum(dpois(dep.var[dri$paisag=="rg"],lambda=a3+b3*log10(dri$area[dri$paisag=="rg"]),log=T))
	LL1+LL2+LL3
      }

##Model E, with same slope for the landscapes with patch-area effect
LL.mod7.id=function(a1,a2,a3,b23,k){
	LL1=-sum(dnbinom(dep.var[dri$paisag=="tap"],mu=a1,size=k,log=T))
	LL2=-sum(dpois(dep.var[dri$paisag=="cau"],lambda=a2+b23*log10(dri$area[dri$paisag=="cau"]),log=T))
	LL3=-sum(dpois(dep.var[dri$paisag=="rg"],lambda=a3+b23*log10(dri$area[dri$paisag=="rg"]),log=T))
	LL1+LL2+LL3
      }

##Model E, with different slopes for the landscapes with patch-area effect
LL.mod8.id=function(a1,a2,a3,b2,b3,k){
	LL1=-sum(dnbinom(dep.var[dri$paisag=="tap"],mu=a1,size=k,log=T))
	LL2=-sum(dpois(dep.var[dri$paisag=="cau"],lambda=a2+b2*log10(dri$area[dri$paisag=="cau"]),log=T))
	LL3=-sum(dpois(dep.var[dri$paisag=="rg"],lambda=a3+b3*log10(dri$area[dri$paisag=="rg"]),log=T))
	LL1+LL2+LL3
      }

## 2.2.
## Function to fit the models by minimizing log-likelihood functions
## and return the Model selection Table fo AICc
## Models were fit by minimization of the log-likelihood functions above or
## with the formula syntax of the mle2 function
seleciona.mod.id=function(...){
	##coeficients starting values
        cf1=as.numeric(coef(fitdistr(dep.var[dri$paisag!="cau"],"negative binomial")))
        cf3=as.numeric(coef(glm(dep.var~dri$paisag/log10(dri$area)-1,family=poisson(link="identity"))))
        ##model fitting with glm and mle2 optimization
        ## model A
	mod1= mle2(dep.var~dnbinom(mu=media,size=k),parameters=list(media~1,k~1),start=list(media=cf1[2],k=cf1[1]))
        ## model B
        mod2= mle2(dep.var~dpois(lambda=a1+a2*log10(dri$area)),start=list(a1=cf3[2],a2=cf3[5]))
        ## model C
        mod3= mle2(dep.var~dnbinom(mu=media,size=k),parameters=list(media~dri$paisag,k~1),start=list(media=cf1[2],k=cf1[1]),...)
        ## model D, same slopes
        mod4= mle2(dep.var~dpois(lambda=a1+a2*log10(dri$area)),parameters=list(a1~dri$paisag),start=list(a1=cf3[2],a2=cf3[5]))
        ## model F
        mod5= mle2(LL.mod5.id,start=list(a1=cf1[2],a2=cf3[2],a3=cf1[2],b2=cf3[5],k=cf1[1]))
        ## model D, different slopes
        mod6= mle2(LL.mod6.id,start=list(a1=cf3[1],a2=cf3[2],a3=cf3[3],b1=cf3[4],b2=cf3[5],b3=cf3[6]))
        ## model E, same slope
        mod7= mle2(LL.mod7.id,start=list(a1=cf1[2],a2=cf3[2],a3=cf3[3],b23=mean(cf3[5:6]),k=cf1[1]))
        ## model E, different slopes
        mod8= mle2(LL.mod8.id,start=list(a1=cf1[2],a2=cf3[2],a3=cf3[3],b2=cf3[5],b3=cf3[6],k=cf1[1]))
        ##Results
	list(tabela=AICctab(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,sort=T,delta=T,weights=T,nobs=length(dep.var)),
		m1=mod1,m2=mod2,m3=mod3,m4=mod4,m5=mod5,m6=mod6,m7=mod7,m8=mod8)
      }

## 2.3.
## Model fiting and model selection table
## Assigning the desired response variable to the object dep.var
## and run the analysis
## Abundance of endemic species
dep.var <- sem.didelphis$ab.end
ab.end <- seleciona.mod.id()
## Abundance of non-endemic species
dep.var <- sem.didelphis$ab.n.end
ab.n.end <- seleciona.mod.id()