---
title: "Effects of environmental and individual
variation on time patterns of extinction in
small experimental populations - extinction analysis" 
contributors:
  - name: Souleyman Bakker (corresponding author)
    affiliation: Institut d’Écologie et des Sciences de l’Environnement de Paris     (iEES-Paris), Sorbonne Université, CNRS, INRAe, IRD, Université Paris Créteil, Université Paris cité, 75005, Paris, France
    email: souleyman1@hotmail.fr
  - name: Tom J.M. Van Dooren
    affiliation: Institut d’Écologie et des Sciences de l’Environnement de Paris (iEES-Paris), Sorbonne Université, CNRS, INRAe, IRD, Université Paris Créteil, Université Paris cité, 75005, Paris, France
    email: tvdooren@gmail.com
  - name: Thomas Tully
    affiliation: Institut d’Écologie et des Sciences de l’Environnement de Paris (iEES-Paris), Sorbonne Université, CNRS, INRAe, IRD, Université Paris Créteil, Université Paris cité, 75005, Paris, France
    email: thomas.tully@sorbonne-universite.fr
date: "`r format(Sys.Date(), '%B %d, %Y')`"
output: pdf_document
---


# Libraries

```{r setup, include=FALSE}

require("knitr")
knitr::opts_chunk$set(echo = TRUE)
#opts_knit$set(root.dir = )

```

```{r,ECHO=FALSE}
# Data handling libraries
library(lubridate)
library(tidyr)
library(dplyr) 

# Graphic libraries
library(ggplot2)
library(cowplot) 
library(gridExtra)

# Survival analysis libraries
library(muhaz) 
library(MASS) 
library(bshazard) 
library(survminer)
library(survival) 
library(ggsurvfit) 
library(interval) 
library(BiocManager) 
library(mgcv) 
library(HelpersMG)
library(Hmisc)
library(multcomp)
library(gee)
library(lme4)

```

# Data formatting

```{r Import and format data,echo=FALSE}
## Import and data entry for the third experiment 

#setwd()
data3<- read.table(file="Donnees_janvier_2018_D.txt" ,header=T, sep="\t", dec=",", stringsAsFactors = T)
tail(data3)

dim(data3) # 290 rows 104 variables

length(unique(data3$boite)) # 290

# Date formatting
data3$DateTime_demarrage=dmy_hm(data3$Date_demarrage)

# Observations after October removed
tempdata<-data3[,-c(93:104)]
dim(tempdata)
tempdata[,9:92][tempdata[,9:92]>0]<-1 # replace all larger than zero by one

# variables counting observation numbers
npops<-nrow(tempdata) # number of lines
min1<-rep(0,npops) # first one observed per population
max1<-rep(0,npops) # last instance of "1" per population
maxNA<-rep(0,npops) # last instance of NA per population
for(i in 1:npops){
	ones<-which(tempdata[i,9:92]==1)
  	min1[i]<-min(ones) # observations in columns 9 to 92
  	max1[i]<-max(ones)
  	maxNA[i]<-max(which(is.na(tempdata[i,9:92])))
}
# check the result
max1-min1
(maxNA-max1)>0 # there are some populations which were overlooked for a while
tempdata$boite[(maxNA-max1)>0]

# min1 is sometimes Inf, max1 sometimes -Inf when there are no ones
# fill intermittent zeroes for further calculations
for(i in 1:nrow(tempdata)){
  if(min1[i]<Inf)profile<-tempdata[i, (8+min1[i]):(8+max1[i])]
profile<-ifelse(profile==0,1,profile)
tempdata[i, (8+min1[i]):(8+max1[i])]<-profile
# should respect the NA's 
}

min0<-rep(0,npops) # first zero: extinction time
for(i in 1:npops){ # determine min0
min0[i]=min(which(tempdata[i,9:92]==0))
}

tempdata[(min0<min1)&min1<Inf,] # no leading zeroes remaining before the population was initiated
tempdata$min1<-min1
tempdata$max1<-max1
tempdata$maxNA<-maxNA
tempdata$min0<-min0 

#Date formatting
datetime<-dmy_hm(gsub("^D(.*)_a.*$","\\1",names(tempdata[,9:92]))) 
beforeafter<-as.factor(gsub("^.*_.*_(.*)$","\\1", names(tempdata[,9:92])))

intervaldays<-diff(datetime)/(60*60*24) 
plot(intervaldays,type="l",ylab="interval lengths in days",xlab="Interval order") # interval lengths in days

tempdata$extdatetime<-datetime[min0]
tempdata$extobscat<-beforeafter[min0]
tempdata$startdatetime<-dmy_hm(gsub("^D(.*)_a.*$","\\1",tempdata$Date_demarrage))
tempdata$lifespan<-ifelse(min0==Inf,
interval(tempdata$startdatetime,datetime[84])%>%as.numeric('days'),interval(tempdata$startdatetime,datetime[min0])%>%as.numeric('days')) # crude way, not interval-censored, at right limit interval
tempdata$censor<-ifelse(min0==Inf,0,1)
startdt<-ifelse(max1==-Inf,0,interval(tempdata$startdatetime,datetime[max1])%>%as.numeric('days')) # for interval-censored analysis, start of interval where ext occurs
tempdata$tstart<-ifelse(tempdata$censor==0,tempdata$lifespan,startdt)
tempdata$tend<-ifelse(tempdata$censor==0,Inf,tempdata$lifespan)

datasurv<-tempdata # data frame for survival analysis

# LONG FORMAT
#############

# The next line turns this into long format

data3l<-datasurv %>% dplyr::select(-c(mort,Date_deniere_observation)) %>% gather("DateTimeAA","Presence", "D12.01.2018_15h00_avant":"D24.10.2018_15h00_apres", factor_key=T, na.rm=T) 

dim(data3l) 
# 22606    20

# new variables
data3l$AvantApres=as.factor(gsub("^.*_.*_(.*)$","\\1", data3l$DateTimeAA))
data3l$DateTime=dmy_hm(gsub("^D(.*)_a.*$","\\1", data3l$DateTimeAA)) 
data3l$TotalDateTime=data3l$DateTime-data3l$startdatetime # this is duration since the start of experiment 

data3lap<-filter(data3l,AvantApres=="apres")## We only take the observations after culling for several of the analyses so we delete before lines

data3lap$atrisk<-ifelse(!is.na(data3lap$extdatetime)&(data3lap$DateTime>data3lap$extdatetime),NA,1)
head(data3lap)

data3lbp<-subset(data3lap,!is.na(atrisk)) # only observations at risk
# sort by box
data3lbp<-arrange(data3lbp,boite)
# week index

subset(data3lbp,boite=="U_A_05") # example

data3lbp$week<-sequence(table(data3lbp$boite))
data3lbp$t<-data3lbp$week
data3lbp$resp<-ifelse(data3lbp$Presence==0,1,0) # one = extinction

tempdatagam<-data.frame() ##temporary dataframe
for (i in levels(data3lbp$boite)){ 
  subd<-filter(data3lbp,boite==i) 
  subd$interval <- c(subd$DateTime[1]- subd$startdatetime[1],subd$TotalDateTime[-1]-subd$TotalDateTime[-nrow(subd)])  # in days
  subd$midpoint <- as.numeric(subd$TotalDateTime-subd$interval/2)/(60*60*24*7) # 
 tempdatagam<-bind_rows(tempdatagam,subd)
}
data3lbp<-tempdatagam
## data3lbp is a df where a time step equal another culling, the interval length between two cullings is took in account in the analysis 

data3lbp$TimeWeeks<- data3lbp$TotalDateTime
units(data3lbp$TimeWeeks)<-"weeks"
data3lbp$IntervalWeeks<- data3lbp$interval
units(data3lbp$IntervalWeeks)<-"weeks"
data3lbp<-subset(data3lbp,as.numeric(data3lbp$IntervalWeeks)<3) # only observations at risk and excluding five boxes that were omitted for several weeks


data3lbp10<-subset(data3lbp,t<11) # data subset ten cullings or less

```

# Semi-parametric or Non-parametric/parametric comparison of multiple survival curves
```{r Survival curves comparisons}

#Proportionnal hazard assumption test

coxm<-coxph(Surv(lifespan,censor)~factor(predation)*stade_intro, data=datasurv) 
cox.zph(coxm)

# Assumptions are violated for culling effect, interaction culling and founder stage effects and for global data

# Interval censoring comparisons of survival curves

datasurv_interval <- Surv(datasurv$tstart,datasurv$tend,type="interval2")
global_test<-ictest(datasurv_interval~interaction(datasurv$traitement,datasurv$stade_intro),scores="logrank1") #Global test : all the survival curves distribution are not equal
global_test

# We set the culling intensity parameter to compare all founder ages 
# Test : For a specific culling treatment, survival curves distribution are equal
ictest(Surv(tstart,tend, type="interval2")~stade_intro,data=datasurv,subset=traitement=="C") #ns
ictest(Surv(tstart,tend, type="interval2")~stade_intro,data=datasurv,subset=traitement=="U") # 0.001
ictest(Surv(tstart,tend, type="interval2")~stade_intro,data=datasurv,subset=traitement=="D") #***
ictest(Surv(tstart,tend, type="interval2")~stade_intro,data=datasurv,subset=traitement=="T") #***
# survival distributions not equal except for control treatment

# We fix the founder age to compare all culling treatments per founder age
# Test : For a specific founder stage, survival curves distribution are equal bewteen culling treatments 
ictest(Surv(tstart,tend, type="interval2")~traitement,data=datasurv,subset=stade_intro=="A") #***
ictest(Surv(tstart,tend, type="interval2")~traitement,data=datasurv,subset=stade_intro=="J") #***
ictest(Surv(tstart,tend, type="interval2")~traitement,data=datasurv,subset=stade_intro=="N") #***
# survival distributions not equal between culling treatment for a specific founder stage

#For all culling treatments, there was no difference between small adults and large adults. There is a difference between large adults and newborns.

```
# Improving the BP predictions and estimate delta

## Parameters for model and import model functions
```{r Model parameters}
source("BP_Functions_Oikos.R") #Import functions with delta and pi.
# Reminder: we will use the function "Parrays" and the values pext in there.
# predictions using Parrays, warnings concerning ultextinct irrelevant.
# Dont forget to import also parms!
invlogit<-function(x)exp(x)/(1+exp(x)) # for variables which have to stay within ]0,1[
logit<-function(x)log(x/(1-x))


# Model parameters : s1 -> s6 : survival of stage 1 to stage 2, stage 2 to stage 3 etc...
# f1 -> f6 : fecundity for stage 1 to 6, averaged over all individuals
# 1 - gamma = culling intensity
# delta : proportion of missed clutches
# p1 -> p6 : probability that an individual produces a clutch.
# T : upper time limit of forward simulation
# initn : number of initial individual corresponding to the initial number of branching

paramsfixed16<-list(s1 = 0.994,  s2 = 0.997,  s3 = 0.977,  s4 = 0.943,  s5 = 0.913,  s6 = 0.887,
  f1 = 0,  f2 = 0.174,  f3 = 12.311,  f4 = 52.505,  f5 = 82.520,  f6 = 96.619,  gamma = 1, delta=0,p1=1,p2=1,p3=1,p4=1,p5=1,p6=1,
  T = 16,  initn = 1)

paramsfixed41<-list(s1 = 0.994,  s2 = 0.997,  s3 = 0.977,  s4 = 0.943,  s5 = 0.913,  s6 = 0.887,
  f1 = 0,  f2 = 0.174,  f3 = 12.311,  f4 = 52.505,  f5 = 82.520,  f6 = 96.619,  gamma = 1, delta=0,p1=1,p2=1,p3=1,p4=1,p5=1,p6=1,
  T = 41,  initn = 1)
```

## Functions for predictions, optimisations and model comparisons
```{r Model functions, ECHO = FALSE}

parms<-list(s1 = 0.994,  s2 = 0.997,  s3 = 0.977, s4 = 0.943, s5 = 0.913, s6 = 0.887,  f1=0,  f2 = 0, f3 = 12.311, f4 = 52.505, f5 = 82.52, f6 = 96.619, gamma=1,  delta=0, p1 = 1, p2 = 1, p3 = 1, p4 = 1, p5 = 1, p6 = 1, T = 16,  initn = 1)


# predictext generates predicted values from a parameter vector which can take on any real value
# predict with six arguments always

predictext <- function(vec){  # for delta , s2, s4 f2 p4 and p6
  params<- parms
  params$delta<-invlogit(vec[1])
  params$s2<-invlogit(vec[2])
  params$s4<-invlogit(vec[3])
  params$f2<-exp(vec[4])
  params$p4<-invlogit(vec[5])
  params$p6<-invlogit(vec[6])
 
  # we need to predict extinctions for 3 initial stages and four predation levels

  params1<-params;params1$gamma=1;
  params75<-params;params75$gamma=.75;
  params5<-params;params5$gamma=.5;
  params25<-params;params25$gamma=.25;

  extgamma1<-Parrays(params1)$pext[,c(2,4,6)]
  extgamma075<-Parrays(params75)$pext[,c(2,4,6)]
  extgamma050<-Parrays(params5)$pext[,c(2,4,6)]
  extgamma025<-Parrays(params25)$pext[,c(2,4,6)]
  
T<-params$T
  pred<-c(extgamma1,extgamma075,extgamma050,extgamma025)
  stade_intro<-factor(rep(rep(c("N","J","A"),each=T+1,4)))
  predation<-factor(rep(c(0,0.25,0.5,0.75),each=(T+1)*3))
  timevec<-c(0:T)
  t<-rep(timevec,12)
    predictions<-data.frame(pred,stade_intro,predation,t=t)
  return(predictions)} 


predictsurv <- function(vec){ # # for delta , s2, s4 f2 p4 and p6
  params<- parms 
params$delta<-invlogit(vec[1])
  params$s2<-invlogit(vec[2])
  params$s4<-invlogit(vec[3])
  params$f2<-exp(vec[4])
  params$p4<-invlogit(vec[5])
  params$p6<-invlogit(vec[6])
 
  # we need to predict survival for 3 initial stages and four predation levels
  params1<-params;params1$gamma=1;
  params75<-params;params75$gamma=.75;
  params5<-params;params5$gamma=.5;
  params25<-params;params25$gamma=.25;
  extgamma1<-Parrays(params1)$psurv[,c(2,4,6)]
  extgamma075<-Parrays(params75)$psurv[,c(2,4,6)]
  extgamma050<-Parrays(params5)$psurv[,c(2,4,6)]
  extgamma025<-Parrays(params25)$psurv[,c(2,4,6)]
  
  T<-params$T
  pred<-c(extgamma1,extgamma075,extgamma050,extgamma025)
  stade_intro<-factor(rep(rep(c("N","J","A"),each=T+1,4)))
  predation<-factor(rep(c(0,0.25,0.5,0.75),each=(T+1)*3))
  timevec<-c(0:T)
  t<-rep(timevec,12)
    predictions<-data.frame(pred,stade_intro,predation,t=t)
  return(predictions)} 

loglik<-function(vec){
  # compute the log likelihood of the model compared to empirical data
  predictions<-predictext(vec) # make the model predictions for "params"
  
  comb<-merge(predictions,data3lbp10)
   # print(dim(comb))
sum(dbinom(x=comb$resp[!is.na(comb$resp)],size=1,prob=comb$pred[!is.na(comb$resp)],log=T))
  
}

deviance<-function(vec,data){
  # compute the deviance of the model compared to empirical data
  predictions<-predictext(vec) # make the model predictions for "params"
  
  comb<-merge(predictions,data)
  
  -2*sum(dbinom(x=comb$resp[!is.na(comb$resp)],size=1,prob=comb$pred[!is.na(comb$resp)],log=T)) 

  #Indicates the fit of the ML model. The lower the fit, the larger the sum of the residual squares (distance between observed and predicted data). The aim is to reduce this value to obtain a better fit.
}

# Calculate confidence intervals (e.g., using z-scores for 95% CI)
z <- qnorm(0.975)  # 0.975 for 95% CI (two-tailed)

# versatile optimizer
dsfpoptim<-function(names){
initvec<-c(-3,rep(1,5))
tochange<-match(names,c("delta","s2","s4","f2","p4","p6"))
initint<-initvec[tochange]
logliki<-function(vecint){
vec<-unlist(parms)[c(14,2,4,8,18,20)] # 18 and 20 are always one
vec<-c(logit(vec[1:3]),log(vec)[4],17,17)
vec[tochange]<-vecint
loglik(vec)
}

out<- optim(par=initint, fn=logliki, hessian=TRUE,method = "L-BFGS-B",control=list(fnscale=-1)) 
standard_errors<-SEfromHessian(-out$hessian, silent = F)
confidence_intervals <- out$par + outer(standard_errors, c(-z, z))
parameters_CI <- data.frame(out$par,confidence_intervals) # not back-transformed!
colnames(parameters_CI) = c("Parameter","Lower 95%", "Upper 95") 
rownames(parameters_CI )<-names
list(out,standard_errors,parameters_CI)
}
# Optimising multiple parameters
Fitdsfp<-dsfpoptim(c("delta","s2","s4","f2","p4","p6"))

Fitdsf<-dsfpoptim(c("delta","s2","s4","f2"))
Fitdsp<-dsfpoptim(c("delta","s2","s4","p4","p6"))
Fitdfp<-dsfpoptim(c("delta","f2","p4","p6"))
Fitfsp<-dsfpoptim(c("s2","s4","f2","p4","p6"))

Fitds<-dsfpoptim(c("delta","s2"))
Fitdf<-dsfpoptim(c("delta","f2"))
Fitdp<-dsfpoptim(c("delta","p4","p6"))
Fitsf<-dsfpoptim(c("s2","s4","f2"))
Fitsp<-dsfpoptim(c("s2","s4","p4","p6"))
Fitfp<-dsfpoptim(c("f2","p4","p6"))

Fitd<-dsfpoptim(c("delta"))
Fits<-dsfpoptim(c("s2","s4"))
Fitf<-dsfpoptim(c("f2","p4","p6"))
Fitp<-dsfpoptim(c("p4","p6"))

vecnull<- c(-Inf,logit(c(parms$s2,parms$s4)),-Inf,17,17)
Fitnull<-loglik(c(-Inf,logit(c(parms$s2,parms$s4)),-Inf,17,17))

# model comparisons with AIC
AICccalc<-function(obj){
# -2*obj[[1]]$val+2*nrow(obj[[3]]) # AIC
(-2 * obj[[1]]$val) + (2 * nrow(obj[[3]])) * (1 + ((nrow(obj[[3]]) + 1) / (1799 - nrow(obj[[3]]) - 1))) # 1799 observations in this experiment
}

AICcnull<-Fitnull*(-2)
AICccalc(Fitdsfp) # <- best AIC
AICccalc(Fitdsf)
AICccalc(Fitdsp)
AICccalc(Fitdfp)
AICccalc(Fitfsp)
AICccalc(Fitds)
AICccalc(Fitdf)
AICccalc(Fitdp)
AICccalc(Fitsf)
AICccalc(Fitsp)
AICccalc(Fitfp)
AICccalc(Fitd)
AICccalc(Fits)
AICccalc(Fitf)
AICccalc(Fitp)

rescaleconf<-Fitdsfp[[3]]
rescaleconf[c(1,2,3,5,6),]<-invlogit(rescaleconf[c(1,2,3,5,6),])
rescaleconf[c(4),]<-exp(rescaleconf[c(4),])
rescaleconf
```

## Growthrates to identify subcritical and critical branching in function of culling parameters  
```{r Growthrates}

#Parameters with optimised delta, s1,s2, f2, p4, p6
paramsfixed16_dsfp<-list(s1 = 0.4658531,  s2 = 0.3023434,  s3 = 0.977,  s4 = 0.9365130,  s5 = 0.913,  s6 = 0.887,
  f1 = 0,  f2 = 0.4519951,  f3 = 12.311,  f4 = 52.505,  f5 = 82.520,  f6 = 96.619,  gamma = 1, delta= 0.6128728,p1=1,p2=1,p3=1,p4= 0.9306119 ,p5=1,p6= 0.7932153,
  T = 16,  initn = 1)

# Parameters subsets for each culling intensity

params1_16_dsfp<-paramsfixed16_dsfp;params1_16_dsfp$gamma=1;
params75_16_dsfp<-paramsfixed16_dsfp;params75_16_dsfp$gamma=.75;
params5_16_dsfp<-paramsfixed16_dsfp;params5_16_dsfp$gamma=.5;
params25_16_dsfp<-paramsfixed16_dsfp;params25_16_dsfp$gamma=.25;

#Growth rates for each culling levels. It correspond to log(lambda)
growth_gamma1_dsf=log(growthrate(meanmatrix(params1_16_dsfp)))
growth_gamma75_dsf=log(growthrate(meanmatrix(params75_16_dsfp)))
growth_gamma50_dsf=log(growthrate(meanmatrix(params5_16_dsfp)))
growth_gamma25_dsf=log(growthrate(meanmatrix(params25_16_dsfp)))

#Compute the gamma value at which lambda falls below 1

for(j in seq(1,0,-0.01)){ #Gamma from 1 to 0 with 0.01 increment 
  paramsfixed16_dsfp_test <- paramsfixed16_dsfp;paramsfixed16_dsfp_test$gamma=j
  growth_gammatest <- growthrate(meanmatrix(paramsfixed16_dsfp_test))
  if(growth_gammatest<1){
    return(j)
  }
  else{}
}
# for gamma = .4 the process become supercritical 


```



# Statistical modelling of extinction probabilities over time using data3lbp
#Plot extinction probabilities per culling event, placing points at number of days since founding. Note the widths of confidence intervals, which depend on sample size.

```{r}

data3lbpagg<-aggregate(cbind(atrisk,resp) ~ week + predation + stade_intro+TimeWeeks+IntervalWeeks+midpoint, data = data3lbp,FUN= function(x) sum(x))
data3lbpagg$extp<- data3lbpagg$resp/data3lbpagg$atrisk


binci<-binconf(data3lbpagg$resp, data3lbpagg$atrisk)# Wilson method
data3lbpagg$lowCI<- binci[,2]
data3lbpagg$highCI<- binci[,3]
head(data3lbpagg)
data3lbpagg$highCI<- binci[,3]
data3lbpagg$numweeks<-as.numeric(data3lbpagg$TimeWeeks)

pplot <- ggplot(data3lbpagg, aes(numweeks, extp, colour = stade_intro))
pplot + geom_point()+ geom_linerange(aes(ymin = lowCI, ymax = highCI)) + facet_wrap(~predation, ncol = 4) # width increase as the sample size is decreasing 


```

## Investigation of time effects in general. 
# Categorical variable for two periods: 
# probabilities of extinction predicted by BP are non-zero, or they are zero
# first analysis: counts after extinction, combining culling and autonomous extinction.
````{r Time effects}
# we could try a gam here, but can't integrate interactions with intervals, so a glm is used 
# create period variable
# We need parms with T=41 for the predicted values for the graphs 
parms <- list(s1 = 0.994,  s2 = 0.997,  s3 = 0.977,  s4 = 0.943,  s5 = 0.913,  s6 = 0.887,
  f1 = 0,  f2 = 0.,  f3 = 12.311,  f4 = 52.505,  f5 = 82.520,  f6 = 96.619,  gamma = 1, delta=0,p1=1,p2=1,p3=1,p4=1,p5=1,p6=1,
  T = 41,  initn = 1)


predictedforglm<- predictext( Fitdsfp[[1]]$par)
predictedforglm2<- predictext( c(Fitd[[1]]$par,unlist(parms)[c(2,4,8,18,20)]))
names(predictedforglm)[c(4)]<-c("week")
data3glmagg<-merge(predictedforglm , data3lbpagg)

# if rate is 10^-9 or smaller, probability to observe an extinction in experiment is
1-((1-10^-9)^290)^41
1-((1-10^-6)^290)^41 # less than 5% probability
# Using cutoff of extinction probability would make that some populations are never small. Others are never considered established or at QSD
# a ten or twenty week limit seems more reasonable in the end, to have appreciable amounts of data points in both periods
data3glmagg$period<-ifelse(data3glmagg$week<20 ,"Small","Estab") #
data3glmagg$numinterval<-as.numeric(data3glmagg$IntervalWeeks)



# Modelling survival and extinction probabilities after an interval + culling event.
# model for survival probability
Gmmext0<-glm(I(atrisk-resp)~offset(log(atrisk))+ (period*predation+ stade_intro*predation+period*stade_intro)*(log(numweeks)+midpoint*numinterval),family="poisson",data=data3glmagg)
Gmmext1<-stepAIC(Gmmext0)
summary(Gmmext1)
drop1(Gmmext1,test="C")
Gmmext2<-update(Gmmext1,.~.-stade_intro:midpoint- period:midpoint -period-midpoint)
drop1(Gmmext2,test="C")
summary(Gmmext2)

# extinction probability per population
data3glmagg$predsgmm2<-1-exp(predict(Gmmext2)- (log(data3glmagg$atrisk)))

pplot <- ggplot(data3lbpagg, aes(numweeks, extp, colour = stade_intro))
pplot + geom_point()+ geom_linerange(aes(ymin = lowCI, ymax = highCI)) +
geom_line(aes(numweeks, predsgmm2),data=data3glmagg,lwd=1)+facet_wrap(~predation, ncol = 4)
# very rough and not very useful
# separate modelling before and during culling uses more data. Can deal with all variables and interactions better and with constraints.

```

## Format for gam/glm analysis separate periods
```{r Formatting for gam and glm}
# We add the number of populations at risk, the numbers of populations alive, the length of each interval and its midpoint. 
data3lgam<-data3l

dim(data3lgam) # 22606    23

data3lgam$atrisk<-ifelse(!is.na(data3lgam$extdatetime)&(data3lgam$DateTime>data3lgam$extdatetime),NA,1)
data3lgam<-subset(data3lgam,!is.na(atrisk)) # only observations at risk
# sort by box
data3lgam<-arrange(data3lgam,boite)
data3lgam$AvantApres<-relevel(data3lgam$AvantApres,"avant")
data3lgam$TimeWeeks<- data3lgam$TotalDateTime
units(data3lgam$TimeWeeks)<-"weeks"
data3lgam$resps<-ifelse(data3lgam$Presence==1,1,0) # one = survival, analyse probability of survival
data3lgam$respe<-ifelse(data3lgam$Presence==0,1,0) # one = extinction, analyse probability of extinction


tempdatagam<-data.frame() ##temporary dataframe
for (i in levels(data3lgam$boite)){ 
  subd<-filter(data3lgam,boite==i) 
subd<-subd[order(subd$TotalDateTime,subd$AvantApres),]
  subd$interval <- as.numeric(subd$TotalDateTime-c(0,subd$TotalDateTime[-nrow(subd)]))/(60*60*24)  # in days
  subd$midpoint <- as.numeric(subd$TotalDateTime)/(60*60*24) -subd$interval/2 # 
 tempdatagam<-bind_rows(tempdatagam,subd)
}
data3lgam<-tempdatagam
data3lgam<-subset(data3lgam, midpoint>0) # remove founding events

data3lgam$IntervalWeeks<- data3lgam$interval/7
data3lgam<-subset(data3lgam,as.numeric(data3lgam$IntervalWeeks)<3) # only observations at risk and excluding five boxes that were omitted for several weeks

# add periods, first combine with data3glmagg
data3lwagg<-data3glmagg[,c(1,2,5,14)]
data3lgamp<-merge(data3lgam, unique(data3lwagg))
data3lgamp$period<-factor(data3lgamp$period)

dim(data3lgamp) # 9537   32

forgambefore<-subset(data3lgamp,subset=AvantApres=="avant"&!is.na(interval)) #before culling
forgamafter<-subset(data3lgamp,subset=AvantApres=="apres"&!is.na(interval)) #after culling 

table(forgambefore$respe)
table(forgamafter$respe)

# group data for Nsurv = prob Nrisk, analysis after
by_gamgroup<-forgamafter%>%group_by(stade_intro,period,traitement,midpoint)
forgamafterN<-by_gamgroup%>%summarise(across(c("atrisk", "resps"),function(x) sum(x)))

#inspect
data.frame(subset(forgamafterN,stade_intro=="J"&traitement=="U"))

# Looking for the boxes going through small then established

box_forarticle <- data3lgamp %>%
  group_by(boite) %>%
  filter(all(c("Estab", "Small") %in% period)) %>%
  ungroup() 

unique(box_forarticle$boite)

boxes_s_to_est <- subset(data3lgamp, boite %in% unique(box_forarticle$boite)) #tab with these boxes


data3lgamp$ID <- paste(data3lgamp$traitement,data3lgamp$replicat,sep="_")
unique(data3lgamp$ID)

gee_1 <- gee(respe~ stade_intro*traitement,id=ID,family=binomial,data=data3lgamp,corstr = "independence",na.action = na.omit)

```

## GAMs and Glms before culling, estimates log survival probabilities

```{r glm pre culling extinctions}
Gamavant0n<-gam(respe~offset(log(interval))+stade_intro*traitement*period+s(midpoint, by=interaction(stade_intro:traitement:period)),family=binomial(link="cloglog"),data=forgambefore,control = list(maxit=200),method= "ML")
plot(Gamavant0n,scale=F)
summary(Gamavant0n) # only one spline significant, N in U


glmavant0<-glm(respe~offset(log(interval))+stade_intro*traitement*period*midpoint,family=binomial(link="cloglog"),data=forgambefore,control = list(maxit=200)) #Signfiicativ full interaction 
glmavant1<-stepAIC(glmavant0)
updateglmavant<-glmavant1
drop1(updateglmavant,test="Chisq") ## three way significant

summary(updateglmavant)      
#The tab i used for interculling extinction modeling

#GLMM with box ID as a random factor effect 
forgambefore$ID <- paste(forgambefore$traitement,forgambefore$replicat,sep="_")

forgambefore_full_random <- glmer(respe~offset(log(interval))+stade_intro*traitement*period*midpoint+(1|ID),family=binomial,data=forgambefore)

gee_2 <- gee(respe~ offset(log(interval))+stade_intro*traitement*period*midpoint,id=ID,family=binomial,data=data3lgamp,corstr = "independence",na.action = na.omit)
```

## Glms moment of culling
```{r glm post culling extinctions}
Gamafter0p<-gam(resps~stade_intro*traitement*period+s(midpoint, by=interaction(stade_intro:traitement:period)),family=poisson,data=forgamafter,control = list(maxit=200),method= "ML")
plot(Gamafter0p,scale=F) # counting process
Gamafter0b<-gam(respe~stade_intro*traitement*period+s(midpoint, by=interaction(stade_intro:traitement:period)),family=binomial,data=forgamafter,control = list(maxit=200),method= "ML")
plot(Gamafter0b,scale=F)

# classical glm, estimates probability of extinction
glmafter_full <- glm(respe~ stade_intro*traitement*period*midpoint,family=binomial,data=forgamafter)
summary(glmafter_full)


updateglmafter<-stepAIC(glmafter_full)
drop1(updateglmafter,test="Chisq") #significant third order interaction
# and treatment effect
summary(updateglmafter)
updateglmafter

#GLMM to test random effect of the treatments 
forgamafter$ID <- paste(forgamafter$traitement,forgamafter$replicat,sep="_")

glmafter_full_random <- glmer(respe~ stade_intro*traitement*period*midpoint+(1|ID),family=binomial,data=forgamafter)

summary(glmafter_full_random)

# Tukey comparisons 
confint(glht(updateglmafter, linfct = mcp("traitement" = "Tukey")))


# Data for figure, predicted extinction probabilities at each culling event
# predicted values, probabilities of extinction
avpred<-predict(updateglmavant,se.fit=T) # extinction 
forgambefore$predictav<-avpred$fit
forgambefore$predictavse<-avpred$se.fit
forgraphbefore<-unique(forgambefore[,c(1,2,3,4,30,32,33)])
appred<- predict(updateglmafter,se.fit=T) # extinction
forgamafter$predictaf<-appred$fit
forgamafter$predictafse<-appred$se.fit
forgraphafter<-unique(forgamafter[,c(1,2,3,4,32,33)])

datagraphglm<-merge(data3lbpagg,forgraphbefore)
datagraphglm<-merge(datagraphglm,forgraphafter)

# calculate predicted extinction probability over interval & culling
datagraphglm$predext<-1- (1-invlogit(datagraphglm$predictav))*(1- invlogit(datagraphglm$predictaf))
datagraphglm$predextll<-1- (1-invlogit(datagraphglm$predictav- 2*datagraphglm$predictavse))*(1- invlogit(datagraphglm$predictaf-2*datagraphglm$predictafse))

datagraphglm$predextul<-1- (1-invlogit(datagraphglm$predictav+ 2*datagraphglm$predictavse))*(1- invlogit(datagraphglm$predictaf+2*datagraphglm$predictafse))
# combine with data in figure

```
