# Use file Dispersal Pairs for Synchrony analysis.xlsx

x1<-read.table(file="clipboard",header=T)
attach(x1)

# Load LME4, GLMMadaptive
library(lme4)
library(GLMMadaptive)

fit1<-glmer(Dispersal~Marked+RecapRate+Treat*EXP+(1|Mead),family=poisson)
summary(fit1)
plot(fit1)

fit2<-glmer.nb(Dispersal~Marked+RecapRate+Treat*EXP+(1|Mead))
summary(fit2)
plot(fit2)

anova(fit1,fit2)

fit2b<-glmer.nb(Dispersal~Marked+RecapRate+Treat+EXP+(1|Mead))
summary(fit2b)
plot(fit2b)

anova(fit2b,fit2)

# Calculate marginal means and SE

p2<-predict(fit2b, type="response")
x2<-data.frame(x1,p2)
tapply(p2,list(Treat,EXP),mean)
tapply(p2,list(Treat,EXP),function(x) sqrt(var(p2)/length(p2)))


#Alex ideas

### Separate analysis binary and gradient
library(ggplot2) ## load for figures

## Create a binary response
x1[,"DispersalYN"] <- ifelse(Dispersal > 0, 1, 0)

## A Dispersal event occurred
fit3 <-glmer(DispersalYN~Marked+RecapRate+Treat*EXP+(1|Mead), family="binomial", data=x1)
summary(fit3)
plot(fit3)
car::Anova(fit3, type=3)

## The number of dispersal events
fit4 <-glmer.nb(Dispersal~Marked+RecapRate+Treat*EXP+(1|Mead), data=subset(x1, DispersalYN != 0))
summary(fit4)
plot(fit4)
car::Anova(fit4, type=3)

## Create some plots of analyses
## Dispersal rate based on marked or recapture rate
markEff <- effects::effect( "Marked",fit3, xlevels =  list(Marked= seq(min(Marked), max(Marked),5)), se=T)
markEff <- data.frame(markEff)

ggplot(markEff, aes(x= Marked, y= fit)) + geom_line(lwd=2) + geom_ribbon(aes(ymin= lower, ymax=upper), alpha=0.3) + theme_bw()

recapEff <- effects::effect( "RecapRate",fit3, xlevels =  list(RecapRate= seq(min(RecapRate), max(RecapRate),0.01)), se=T)
recapEff <- data.frame(recapEff)

ggplot(recapEff, aes(x= RecapRate, y= fit)) + geom_line(lwd=2) + geom_ribbon(aes(ymin= lower, ymax=upper), alpha=0.3) + theme_bw()

## Treatment and control on dispersal
ggplot(subset(x1, DispersalYN != 0), aes(x=Treat, y=Dispersal)) + geom_boxplot() + facet_grid(~EXP)

recapEff <- effects::effect( "RecapRate",fit4, xlevels =  list(RecapRate= seq(min(RecapRate), max(RecapRate),0.01)), se=T)
recapEff <- data.frame(recapEff)


### Zero-inflated
library(glmmTMB)

## default model for binary component
fit5 <- glmmTMB(Dispersal~Marked+RecapRate+Treat*EXP+(1|Mead), family=nbinom2, data=x1, ziformula=~1)
summary(fit5)
plot(fitted(fit5), residuals(fit5))
car::Anova(fit5, type=3)

fit6 <- glmmTMB(Dispersal~Marked+RecapRate+Treat*EXP+(1|Mead), family=poisson, data=x1, ziformula=~1)
summary(fit6)
plot(fitted(fit6), residuals(fit6))
car::Anova(fit6, type=3)

anova(fit5, fit6)

fit7<-update(fit6, family=nbinom1)
summary(fit7)
plot(fitted(fit7), residuals(fit7))
anova(fit6,fit7)
anova(fit5,fit6)
AICtab(fit5,fit6,fit7)

fit8 <- glmmTMB(Dispersal~Marked+RecapRate+Treat*EXP+(1|Mead), family=nbinom2, data=x1, ziformula=~Marked+RecapRate)
summary(fit8)
plot(fitted(fit8), residuals(fit8))
car::Anova(fit8, type=3)

anova(fit8,fit5)
fit9<-update(fit8, family=nbinom1)
fit10<-update(fit8, family=poisson)
AICtab(fit8,fit9,fit10,fit5,fit2b)

fit8b <- glmmTMB(Dispersal~Marked+RecapRate+Treat+EXP+(1|Mead), family=nbinom2, data=x1, ziformula=~Marked+RecapRate)
summary(fit8b)
plot(fitted(fit8b), residuals(fit8b))
car::Anova(fit8b, type=3)

anova(fit8,fit8b)

p8b<-predict(fit8b, type="response")
x8b<-data.frame(x1,p8b)
tapply(p8b,list(Treat,EXP),mean)
tapply(p8b,list(Treat,EXP),function(x) sqrt(var(p8b)/length(p8b)))

ggplot(x1, aes(x=EXP, y=Dispersal)) + geom_boxplot() + facet_grid(~Treat)

ggplot(x1, aes(x=EXP, y=p8b)) + geom_boxplot() + facet_grid(~Treat)