shy.lw<-read.table("shy.l-w.txt",header=T)
summary(shy.lw)

length(shy.lw[,1])
length(unique(shy.lw$tag))
min(table(shy.lw$tag))
max(table(shy.lw$tag))
mean(table(shy.lw$tag))
sd(table(shy.lw$tag))


####Lab-Wild relationship####

plot(log(shy.lw$lshy),log(shy.lw$wshy))

cor.test(shy.lw$lshy,shy.lw$wshy,method="spearman")

shy.lw$burrow<-as.factor(shy.lw$burrow)

hist(shy.lw$wshy)

require(MCMCglmm)


#shyness lab-wild relationship

priorX1 = list(R = list(V = 1, n = 0.2),
               G = list(G1 = list(V = 1, n = 0.2),G2=list(V=1,n=0.2)))

shy.l_w<-MCMCglmm(wshy~lshy*gap,random=~tag+burrow,rcov=~units,family="poisson",
                  prior=priorX1,data=shy.lw,nitt=100000,thin=25,burnin=1000,verbose=F)
summary(shy.l_w)
plot(shy.l_w)

posterior.mode(shy.l_w$Sol)
posterior.mode(shy.l_w$VCV)

w.b.rpt<-shy.l_w$VCV[,1]/
  (shy.l_w$VCV[,1]+shy.l_w$VCV[,2]+shy.l_w$VCV[,3]+log(1/ exp(shy.l_w$Sol[,1])+1))
posterior.mode(w.b.rpt)
HPDinterval(w.b.rpt)

#activity lab - wild relationship

ae.lw<-read.table("ae.l-w.txt",header=T)
summary(ae.lw)

length(ae.lw[,1])
length(unique(ae.lw$tag))
min(table(ae.lw$tag))
max(table(ae.lw$tag))
mean(table(ae.lw$tag))
sd(table(ae.lw$tag))

hist(ae.lw$wactive)

plot(ae.lw$lactive,ae.lw$wactive)
cor.test(ae.lw$lactive,ae.lw$wactive,method="spearman")

ae.lw$lactive2<-round(ae.lw$lactive*1000)

priorX2 = list(R = list(V = 1, n = 0.2),
               G = list(G1 = list(V = 1, n = 0.2)))


act.lw<-MCMCglmm(wactive~lactive*gap+minutes,random=~tag,rcov=~units,family="poisson",
                 prior=priorX2,data=ae.lw,nitt=100000,thin=25,burnin=1000,verbose=F)
summary(act.lw)
plot(act.lw)

posterior.mode(act.lw$Sol)
posterior.mode(act.lw$VCV)

w.a.rpt<-act.lw$VCV[,1]/
  (act.lw$VCV[,1]+act.lw$VCV[,2]+log(1/ exp(act.lw$Sol[,1])+1))
posterior.mode(w.a.rpt)
HPDinterval(w.a.rpt)

plot(ae.lw$lexplore,ae.lw$wexplore)
cor.test(ae.lw$lexplore,ae.lw$wexplore,method="spearman")

#exploration lab-wild relationship

exp.lw<-MCMCglmm(wexplore~lexplore*gap+minutes,random=~tag,rcov=~units,family="poisson",
                 prior=priorX2,data=ae.lw,nitt=100000,thin=25,burnin=1000,verbose=F)
summary(exp.lw)
plot(exp.lw)

posterior.mode(exp.lw$Sol)
posterior.mode(exp.lw$VCV)

w.e.rpt<-exp.lw$VCV[,1]/
  (exp.lw$VCV[,1]+exp.lw$VCV[,2]+log(1/ exp(exp.lw$Sol[,1])+1))
posterior.mode(w.e.rpt)
HPDinterval(w.e.rpt)

####Lab repeatability####



bae.13<-read.table("bae.13.txt",header=T)
summary(bae.1314)

bae.13<-subset(bae.13, age!="NA" )
summary(bae.13)

length(bae.13[,1][bae.13$out==1])
length(unique(bae.13$tag[bae.13$out==1]))
min(table(bae.13$tag[bae.13$out==1]))
max(table(bae.13$tag[bae.13$out==1]))
mean(table(bae.13$tag[bae.13$out==1]))
sd(table(bae.13$tag[bae.13$out==1]))

#shyness
shy.l2<-MCMCglmm(round(shy)~sex+scale(age)+scale(temp),random=~tagyear,rcov=~units,family="poisson",
                prior=priorX2,data=bae.13,nitt=100000,thin=25,burnin=1000,verbose=F)
summary(shy.l2)
plot(shy.l2)

posterior.mode(shy.l2$Sol)
posterior.mode(shy.l2$VCV)

l.s.rpt2<-shy.l2$VCV[,1]/
  (shy.l2$VCV[,1]+shy.l2$VCV[,2]+log(1/ exp(shy.l2$Sol[,1])+1))
posterior.mode(l.s.rpt2)
HPDinterval(l.s.rpt2)

#activity
bae.13$lactive<-round(bae.13$rate*1000)
act.l2<-MCMCglmm(lactive~sex+scale(age)+scale(temp),random=~tagyear,rcov=~units,family="poisson",
                 prior=priorX2,data=bae.13,nitt=100000,thin=25,burnin=1000,verbose=F)
summary(act.l2)
plot(act.l2)

posterior.mode(act.l2$Sol)
posterior.mode(act.l2$VCV)

l.a.rpt2<-act.l2$VCV[,1]/
  (act.l2$VCV[,1]+act.l2$VCV[,2]+log(1/ exp(act.l2$Sol[,1])+1))
posterior.mode(l.a.rpt2)
HPDinterval(l.a.rpt2)


#exploration
exp.l2<-MCMCglmm(explore~sex+scale(age)+scale(temp),random=~tagyear,rcov=~units,family="poisson",
                 prior=priorX2,data=bae.13,nitt=100000,thin=25,burnin=1000,verbose=F)
summary(exp.l2)
plot(exp.l2)

posterior.mode(exp.l2$Sol)
posterior.mode(exp.l2$VCV)

l.e.rpt2<-exp.l2$VCV[,1]/
  (exp.l2$VCV[,1]+exp.l2$VCV[,2]+log(1/ exp(exp.l2$Sol[,1])+1))
posterior.mode(l.e.rpt2)
HPDinterval(l.e.rpt2)


####Plots####
require(ggplot2)

shylw.plot <- ggplot(shy.lw, aes(lshy,log(wshy+1)))
shylw.plot + stat_smooth(method = "lm",col="black") + geom_point()+theme(panel.background = element_blank())+
  ylab("Log of wild shyness") +xlab("Lab shyness (seconds)")+
  theme(axis.line = element_line (colour = "black"),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.text.y = element_text (colour = "black", size = 12),
        axis.ticks = element_line (colour = "black"),
        axis.title=element_text(size=12,face="bold"))

actlw.plot <- ggplot(ae.lw, aes(lactive2,wactive))
actlw.plot + stat_smooth(method = "lm",col="black") + geom_point()+theme(panel.background = element_blank())+
  ylab("Wild activity (number of 'leaves' per day)") +xlab("Lab activity (number of wires crossed per 1000 seconds)")+
  theme(axis.line = element_line (colour = "black"),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.text.y = element_text (colour = "black", size = 12),
        axis.ticks = element_line (colour = "black"),
        axis.title=element_text(size=12,face="bold"))

explw.plot <- ggplot(ae.lw, aes(lexplore,wexplore))
explw.plot + stat_smooth(method = "lm",col="black") + 
  geom_point(position = position_jitter(w = 0.1, h = 0.1))+theme(panel.background = element_blank())+
  ylab("Wild exploration (unique burrows per day)") +xlab("Lab exploration (unique tripwires in first minute)")+
  theme(axis.line = element_line (colour = "black"),
        axis.text.x = element_text(colour = "black", size = 14),
        axis.text.y = element_text (colour = "black", size = 14),
        axis.ticks = element_line (colour = "black"),
        axis.title=element_text(size=14,face="bold"))



####Samples sizes in the wild####
baew<-read.table("bae.w2.txt",header=T)
summary(baew)

baew<-subset(baew,age!="NA")
length(baew[,1])
length(unique(baew$tag))
length(unique(baew$burrow))
min(table(baew$tag))
max(table(baew$tag))
mean(table(baew$tag))
sd(table(baew$tag))


aew<-read.table("ae.w.txt",header=T)
summary(aew)

length(aew[,1])
length(unique(aew$tag))
min(table(aew$tag))
max(table(aew$tag))
mean(table(aew$tag))
sd(table(aew$tag))

####Bivariate method####

bbivar<-read.table("bbivar.txt",header=T)
summary(bbivar)

bbivarxn<-subset(bbivar,age!="NA")
summary(bbivarxn)


priorX7 = list(R = list(V = diag(2), n = 0, fix=1),G = list(G1 = list(V = diag(2), n = 2)))

#shyness
lwbbivar<-MCMCglmm(fixed=cbind(round(lshy),wbold)~trait+trait:(sex + scale(age)),
                   random=~us(trait):tag,rcov=~us(trait):units,family=c("poisson","poisson"),
                   prior=priorX7,data=bbivarxn,nitt=100000,thin=25,burnin=1000,verbose=F)
summary(lwbbivar)
plot(lwbbivar)

posterior.mode(lwbbivar$Sol)
posterior.mode(lwbbivar$VCV)

corr.b<-lwbbivar$VCV[,2] / (lwbbivar$VCV[,1]*lwbbivar$VCV[,4])^0.5
posterior.mode(corr.b)
HPDinterval(corr.b)

#without fixed effects, using all data
lwbbivar2<-MCMCglmm(fixed=cbind(round(lshy),wbold)~trait,
                   random=~us(trait):tag,rcov=~us(trait):units,family=c("poisson","poisson"),
                   prior=priorX7,data=bbivar,nitt=100000,thin=25,burnin=1000,verbose=F)


aebivar<-read.table("aebivar.txt",header=T)
summary(aebivar)

aebivar$lactive2<-round(aebivar$lactive*1000)

aebivarxn<-subset(aebivar,age!="NA")
summary(aebivarxn)

hist(aebivar$lactive2)
hist(abivar$wactive)

#activity
lwabivar<-MCMCglmm(fixed=cbind(lactive2,wactive)~trait+trait:(sex + scale(age)),
                   random=~us(trait):tag,rcov=~us(trait):units,family=c("poisson","poisson"),
                   prior=priorX7,data=aebivarxn,nitt=100000,thin=25,burnin=1000,verbose=F)
summary(lwabivar)
plot(lwabivar)

posterior.mode(lwabivar$Sol)
posterior.mode(lwabivar$VCV)

corr.a<-lwabivar$VCV[,2] / (lwabivar$VCV[,1]*lwabivar$VCV[,4])^0.5
posterior.mode(corr.a)
HPDinterval(corr.a)

#without sex and age as fixed effects, using all data
lwabivar2<-MCMCglmm(fixed=cbind(lactive2,wactive)~trait,
                   random=~us(trait):tag,rcov=~us(trait):units,family=c("poisson","poisson"),
                   prior=priorX7,data=aebivar,nitt=100000,thin=25,burnin=1000,verbose=F)

#exploration
lwebivar<-MCMCglmm(fixed=cbind(lexplore,wexplore)~trait+trait:(sex + scale(age)),
                   random=~us(trait):tag,rcov=~us(trait):units,family=c("poisson","poisson"),
                   prior=priorX7,data=aebivarxn,nitt=100000,thin=25,burnin=1000,verbose=F)
summary(lwebivar)
plot(lwebivar)

posterior.mode(lwebivar$Sol)
posterior.mode(lwebivar$VCV)

corr.e<-lwebivar$VCV[,2] / (lwebivar$VCV[,1]*lwebivar$VCV[,4])^0.5
posterior.mode(corr.e)
HPDinterval(corr.e)

#without sex and age as fixed effects, using all data
lwebivar2<-MCMCglmm(fixed=cbind(lexplore,wexplore)~trait,
                   random=~us(trait):tag,rcov=~us(trait):units,family=c("poisson","poisson"),
                   prior=priorX7,data=aebivar,nitt=100000,thin=25,burnin=1000,verbose=F)

####Wild analyses####

bae.w<-read.table("bae.w2.txt",header=T)
summary(bae.w)

fburrow<-as.factor(bae.w$burrow)
length(table(fburrow)) #see how many unique burrows we have
#turn time into a number
pbrush <- as.POSIXct(bae.w$tbrush, format = "%H:%M:%S") 
sbrush<-cut.POSIXt(pbrush,breaks="min",labels=F)
#turn date into a number
pdate <- as.POSIXct(bae.w$date, format = "%d/%m/%Y") 
sdate<-cut.POSIXt(pdate,breaks="day",labels=F)
fpair<-as.factor(bae.w$pair)
#stick them all together
bae.w1<-cbind(bae.w,fburrow,sbrush,sdate,fpair)

summary(bae.w1)
bae.w2<-subset(bae.w1,age!="NA")

length(table(bae.w1$tag))

mean(table(bae.w2$tag))
sd(table(bae.w2$tag))

priorX1 = list(R = list(V = diag(1), n = 2),G = list(G1 = list(V = diag(1), n = 2),G2=list(V=diag(1),n=2)))

wbold1<-MCMCglmm(wbold~sex+scale(age)+fpair+scale(temp)+weather+scale(sdate)+scale(sbrush)
                 +sex:scale(age)+sex:weather+sex:fpair+scale(age):weather+scale(age):fpair,
                 random=~tag+fburrow,rcov=~units,family="poisson",prior=priorX1,
                 data=bae.w2,nitt=50000,burnin=1000,thin=25,verbose=F)
summary(wbold1)
plot(wbold1)

posterior.mode(wbold1$VCV)
posterior.mode(wbold1$Sol)

#repeatability of tag
ICC.wb.t <- wbold1$VCV[,1]/ (wbold1$VCV[,1]+ wbold1$VCV[,2]+wbold1$VCV[,3]+log(1/exp(wbold1$Sol[,1])+1))
posterior.mode(ICC.wb.t)
HPDinterval(ICC.wb.t)

#repeatability of burrow
ICC.wb.b <- wbold1$VCV[,2]/ (wbold1$VCV[,1]+ wbold1$VCV[,2]+wbold1$VCV[,3]+log(1/exp(wbold1$Sol[,1])+1))
posterior.mode(ICC.wb.b)
HPDinterval(ICC.wb.b)


#remove burrow and repeat to see if tag repeatabiity goes up

priorX2 = list(R = list(V = diag(1), n = 2),G = list(G1 = list(V = diag(1), n = 2)))
wbold2<-MCMCglmm(wbold~sex+scale(age)+fpair+scale(temp)+weather+scale(sdate)+scale(sbrush)
                 +sex:scale(age)+sex:weather+sex:fpair+scale(age):weather+scale(age):fpair,
                 random=~tag,rcov=~units,family="poisson",prior=priorX2,
                 data=bae.w2,nitt=50000,burnin=1000,thin=25,verbose=F)
summary(wbold2)
plot(wbold2)
posterior.mode(wbold2$VCV)
#repeatability of tag
ICC.wb.t2 <- wbold2$VCV[,1]/ (wbold2$VCV[,1]+ wbold2$VCV[,2]+log(1/exp(wbold2$Sol[,1])+1))
posterior.mode(ICC.wb.t2)
HPDinterval(ICC.wb.t2)



ae.w<-read.table("ae.w.txt", header=T)
summary(ae.w)

ae.w<-subset(ae.w,age>=0)



#turn date into a number
pdate <- as.POSIXct(ae.w$date, format = "%d/%m/%Y") 
sdate<-cut.POSIXt(pdate,breaks="day",labels=F)

#stick them all together
ae.w1<-cbind(ae.w,sdate)
length(ae.w1$tag)
length(unique(ae.w1$tag))

mean(table(ae.w1$tag))
sd(table(ae.w1$tag))


#only have tag as a random effect as burrow likely to vary over the day
priorX2 = list(R = list(V = diag(1), n = 2),G = list(G1 = list(V = diag(1), n = 2)))

wact1<-MCMCglmm(wactive~scale(minutes)+sex+scale(age)+scale(avtemp)+scale(avsun)+scale(sumrain)+scale(sdate)
                +sex:scale(age)+sex:scale(avsun)+sex:scale(sumrain)+scale(age):scale(avsun)+scale(age):scale(sumrain),
                 random=~tag,rcov=~units,family="poisson",prior=priorX2,
                 data=ae.w1,nitt=50000,burnin=1000,thin=25,verbose=F)
summary(wact1)
plot(wact1)

posterior.mode(wact1$VCV)
posterior.mode(wact1$Sol)

#repeatability of tag
ICC.wa.t <- wact1$VCV[,1]/ (wact1$VCV[,1]+ wact1$VCV[,2]+log(1/exp(wact1$Sol[,1])+1))
posterior.mode(ICC.wa.t)
HPDinterval(ICC.wa.t)

plot(ae.w1$age[ae.w1$sex=="M"],ae.w1$wactive[ae.w1$sex=="M"])
plot(ae.w1$age[ae.w1$sex=="F"],ae.w1$wactive[ae.w1$sex=="F"])

hist(ae.w1$wexplore)
priorX3 = list(R = list(V = diag(2), n = 2),G = list(G1 = list(V = diag(1), n = 2)))

wexp1<-MCMCglmm(wexplore~scale(minutes)+sex+scale(age)+scale(avtemp)+scale(avsun)+scale(sumrain)+scale(sdate)
                +sex:scale(age)+sex:scale(avsun)+sex:scale(sumrain)+scale(age):scale(avsun)+scale(age):scale(sumrain),
                random=~tag,rcov=~us(trait):units,family="zapoisson",prior=priorX3,
                data=ae.w1,nitt=50000,burnin=1000,thin=25,verbose=F)
summary(wexp1)
plot(wexp1)

posterior.mode(wexp1$VCV)
posterior.mode(wexp1$Sol)

#repeatability of tag
ICC.we.t <- wexp1$VCV[,1]/ (wexp1$VCV[,1]+ wexp1$VCV[,2]+log(1/exp(wexp1$Sol[,1])+1))
posterior.mode(ICC.we.t)
HPDinterval(ICC.we.t)