#R Code for Figure 1 & associated analyses

  library(tidyverse)
  library(binom)

  unlogit <- function(x) exp(x) / (1 + exp(x))
  
  dat <- read.csv("Table 2.csv")
  
  dat <- dat %>%
         mutate(n = suc + fail,
                p = suc / n)
 
# calculate confidence interval 
  out <- binom.confint(dat$suc, dat$n, method = "wilson")

  m1 <- glm(cbind(suc, fail) ~ log10(pp), family = binomial, data = dat)
  summary(m1)
  
  yl <- c(0, 1)
  plot(p ~ log10(pp), data = dat, pch = 19, cex = 1.5, bty = "l", ylim = yl,
       xlab = "log10(propagule pressure)", ylab = "Probability of establishment")
    arrows(log10(dat$pp), out$lower, log10(dat$pp), out$upper, length = 0)
    
  xx <- seq(min(log10(dat$pp)), max(log10(dat$pp)), 0.01)
  yy <- unlogit(coef(m1)[1] + coef(m1)[2] * xx)
  
  lines(yy ~ xx, lwd = 2)
  
# convert to 0,1 data to show that this produces the same results as aggregating above
  pp.suc <- rep(dat$pp, dat$suc)
  suc <- rep(rep(1, length(dat$suc)), dat$suc)
  pp.fail <- rep(dat$pp, dat$fail)
  fail <- rep(rep(0, length(dat$fail)), dat$fail)
  
  pa <- data.frame(outcome = c(suc, fail), pp = c(pp.suc, pp.fail))
  dim(pa)
  sum(dat$suc) + sum(dat$fail)
  
  m2 <- glm(outcome ~ log10(pp), family = binomial, data = pa)
  summary(m2)
  summary(m1)
  
  
  
  
  #R Code for Figure 2 & associated analyses
  
  library(arm)
  library(tidyverse)
  library(TeachingDemos)
  library(lme4)
  library(MuMIn)
  
  rm(list=ls())
  
  bird <- read.csv("Moultondata.csv")
  
  x <- bird$lnum
  y <- bird$fate
  
  
  x.lim <- 7
  xx <- 0:x.lim
  
  yy <- ifelse(y==1, 1.05, -0.05)
  plot(jitter(yy, 0.15) ~ x, bty="l", pch=3, cex=0.5, col="grey", 
       xlim=c(0.3, x.lim), ylab="Establishment Probability", xlab="Log10 Propagule Pressure", cex.lab=1.4)
  abline(h=0, col="grey")
  abline(h=1, col="grey")
  
  out.av <- cut_width(x, 1)
  ymean <- tapply(y, out.av, mean)
  xmean <- tapply(x, out.av, mean)
  points(ymean ~ xmean, pch=19, cex=1.3)
  
  # fit simple logistic 
  
  m1 <- glm(bird$fate~bird$lnum, family=binomial(link="logit"))
  summary(m1)
  
  new.lnum = seq(0,6.6,0.1)
  lin <- -1.09852 + (0.32804 *new.lnum)
  p <- exp(lin) / (1+exp(lin))
  lines(p~new.lnum)
  
  
  # fit logistic with random effects
  
  m8 <-glmer(fate ~ lnum  + (1|Location_of_introduction) + (1|Order_.HM4./Family_.HM4.), family = binomial, data = bird)
  summary(m8)
  r.squaredGLMM(m8)
  
  new.lnum2 = seq(0,6.6,0.1)
  lin2 <- -1.6581 + (0.69594 *new.lnum2)
  p <- exp(lin2) / (1+exp(lin2))
  lines(p~new.lnum2, lty = 2)
  
  #fit line for Galliformes only 
  
  bird <- subset(bird, SolOrder == 'Galliformes')
  
  m5 <- glm(bird$fate~bird$lnum, family=binomial(link="logit"))
  summary(m5)
  
  gall.lnum = seq(0,6.6,0.1)
  lin <- -2.2110 + (0.5811 *gall.lnum)
  p <- exp(lin) / (1+exp(lin))
  lines(p~gall.lnum, lty = 3)
  
  
  #fit model for CR species only
  
  CR <- subset(bird, Number < 50)
  
  summary(glmer(fate ~ lnum  + (1|Location_of_introduction) + (1|Order_.HM4./Family_.HM4.), family = binomial, data = CR))