# packages
library(visreg)

# function to perform model selection
m.AICc <- function(models,n){
  LL <- sapply(models,logLik,REML=T) 
  k <- sapply(lapply(models,logLik,REML=T),attr,"df")
  AIC <- -2*LL+2*k
  AICc <- AIC+((2*k*(k+1))/(n-k-1))
  d.AICc <- AICc-min(AICc)
  w.AICc <- (exp(-0.5*d.AICc))/sum(exp(-0.5*d.AICc))
  data.frame(n.par=k,log.lik=round(LL,2),AICc=round(AICc,2),delta.AICc=round(d.AICc,2),w.AICc=round(w.AICc,2),row.names=names(LL))[order(d.AICc),]
}

# input data
data <- read.table("Data_Table_S1.txt", header=T)

# models
m1 <- glm(Nest ~ Native_2km, family=binomial, data=data)
m2 <- glm(Nest ~ Eucalyptus_2km, family=binomial, data=data)
m3 <- glm(Nest ~ Anthropogenic_2km, family=binomial, data=data)
m4 <- glm(Nest ~ Agriculture_2km, family=binomial, data=data)
m5 <- glm(Nest ~ Anthropogenic_2km + Agriculture_2km, family=binomial, data=data)
m6 <- glm(Nest ~ Eucalyptus_2km + Agriculture_2km, family=binomial, data=data)

m7 <- glm(Nest ~ Native_3km, family=binomial, data=data)
m8 <- glm(Nest ~ Eucalyptus_3km, family=binomial, data=data)
m9 <- glm(Nest ~ Anthropogenic_3km, family=binomial, data=data)
m10 <- glm(Nest ~ Agriculture_3km, family=binomial, data=data)
m11 <- glm(Nest ~ Anthropogenic_3km + Agriculture_3km, family=binomial, data=data)
m12 <- glm(Nest ~ Eucalyptus_3km + Agriculture_3km, family=binomial, data=data)
m.null <- glm(Nest ~ 1, family=binomial, data=data)

# model selection
m.AICc(list(m1,m2,m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m.null),24)

# plots
par(mar=c(4.8,5.3,1,1))
visreg(m11, "Anthropogenic_3km", scale="response", rug=2, xlab="Anthropogenic habitat cover (ha)", ylab="Probability of occurrence", line=list(col="black"), cex.axis=1.4, cex.lab=1.7)
text("A", x=350, y=1, cex=1.4)
visreg(m11, "Agriculture_3km", scale="response", rug=2, xlab="Agriculture habitat cover (ha)", ylab="Probability of occurrence", line=list(col="black"), cex.axis=1.4, cex.lab=1.7)
text("B", x=2100, y=1, cex=1.4)