## 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()