library(rjags) library(coda) library(R2WinBUGS) #------# # Data #------# dtot=read.table("Dissection_data.txt",header=T,sep="\t") int=read.table("Ingestion_data.txt",header=T,sep="\t") #---------------# # Time interval #---------------# # Time slot dtot$RT_int=cut(dtot$RT_h, breaks=c(0,1.5,3,6,9,12,15,18,21,24,30,36,42,48,54),right=T,include.lowest=F) # Continuous time hours=c(1,3,6,9,12,15,18,21,24,30,36,42,48,54) for(i in 1:14){ dtot[which(dtot$RT_int==levels(dtot$RT_int)[i]),"RT"]=hours[i] } #-----------------------------------------------# # Without CVC, + weight of the dissecetd sample #-----------------------------------------------# d=subset(dtot,Psp!="CVC") d$Psp=factor(d$Psp) d$faeces=paste("F",d$Faeces,sep="") #-----------------------------------------------------------# # Seeds counted in the dissected sample [animal,time,plant] #-----------------------------------------------------------# ilist=list() for (i in c("CV","JE","PM","PV","RF","TP")){ # For a plant tp=subset(d,Psp==i)[,c("RT","Rep","All_seeds")] # All the faeces collected during a same time slot tp2=aggregate(tp$All_seeds,by=list(tp$RT,tp$Rep),FUN="sum") colnames(tp2)=colnames(tp) # Individual animal * time, in order m1=reshape(tp2[,c("RT","Rep","All_seeds")],direction="wide",idvar="Rep",timevar="RT") m2=m1[,-1] rownames(m2)=m1[,1] colnames(m2)=substring(colnames(m2),first=11,last=nchar(colnames(m2))) m2=m2[rownames(m2),as.character(c(1,3,6,9,12,15,18,21,24,30,36,42,48,54))] m2=cbind(0,m2) # Add to the list ilist[[i]]=m2 } Nf=array(unlist(ilist),dim=c(17,15,6),dimnames=list(rownames(m2),colnames(m2),c("CV","JE","PM","PV","RF","TP"))) #-----------# # Seedlings #-----------# glist=list() for (i in c("CV","JE","PM","PV","RF","TP")){ # For a plant tp=subset(d,Psp==i)[,c("RT","Rep","Germ_seeds")] # All the faeces collected during a same time slot tp2=aggregate(tp$Germ_seeds,by=list(tp$RT,tp$Rep),FUN="sum") colnames(tp2)=colnames(tp) # Individual animal * time, in order m1=reshape(tp2[,c("RT","Rep","Germ_seeds")],direction="wide",idvar="Rep",timevar="RT") m2=m1[,-1] rownames(m2)=m1[,1] colnames(m2)=substring(colnames(m2),first=12,last=nchar(colnames(m2))) m2=m2[rownames(m2),as.character(c(1,3,6,9,12,15,18,21,24,30,36,42,48,54))] m2=cbind(0,m2) # Add to the list glist[[i]]=m2 } Gf=array(unlist(glist),dim=c(17,15,6),dimnames=list(rownames(m2),colnames(m2),c("CV","JE","PM","PV","RF","TP"))) #------------------------------------------------# # Matrix with the weight of the dissected sample #------------------------------------------------# # Weights per time slot and individual dp=unique(d[,c("Rep","RT","Faeces","SFW","FW")]) dp1=aggregate(dp[,c("SFW","FW")],by=list(dp$Rep,dp$RT),FUN="sum") colnames(dp1)=c("Rep","RT","SFW","FW") # Weight of the dissected sample dp2=reshape(dp1[,c("Rep","RT","SFW")],timevar="RT",idvar="Rep",direction="wide") dp2b=dp2[,-1] rownames(dp2b)=dp2[,1] colnames(dp2b)=substring(colnames(dp2b),first=5,last=nchar(colnames(dp2b))) dp2b=cbind(0,dp2b) # Faeces weight dp3=reshape(dp1[,c("Rep","RT","FW")],timevar="RT",idvar="Rep",direction="wide") dp3b=dp3[,-1] rownames(dp3b)=dp3[,1] colnames(dp3b)=substring(colnames(dp3b),first=4,last=nchar(colnames(dp3b))) dp3b=cbind(0,dp3b) #-----------# # Ingestion #-----------# int=subset(int,Ind%in%unique(d$Rep)) int1=int[,-1] rownames(int1)=int[,1] colnames(int1)=c("CV","JE","PM","PV","RF","TP") #---------------------------# # Putting the data in order #---------------------------# Y=Nf[dimnames(Nf)[[1]],dimnames(Nf)[[2]],dimnames(Nf)[[3]]] Yg=Gf[dimnames(Nf)[[1]],dimnames(Nf)[[2]],dimnames(Nf)[[3]]] poids=dp2b[dimnames(Nf)[[1]],dimnames(Nf)[[2]]] poidstot=dp3b[dimnames(Nf)[[1]],dimnames(Nf)[[2]]] Ntot=int1[dimnames(Nf)[[1]],dimnames(Nf)[[3]]] temps=c(0,hours) nanimal=dim(Y)[1] nplante=dim(Y)[3] ntemps=length(temps) ntemps1=ntemps-1 # Seed size t1=unique(d[,c("Psp","SL")]) rownames(t1)=t1[,1] t1=t1[dimnames(Y)[[3]],] taille=as.vector(scale(t1[,2])) #---------------# # Germination #---------------# # Matrix with the weight of the germinated sample (4g for roe deer, 8g for red deer and wild boar) poidsgerm=matrix(0,nrow=17,ncol=15) for(i in 2:ncol(poidsgerm)){poidsgerm[,i]=c(rep(4,5),rep(8,12))} rownames(poidsgerm)=rownames(poids) colnames(poidsgerm)=colnames(poids) #-----------------------------# # Preparing the data for Jags #-----------------------------# poidsRfeces <- poids/poidstot poidsRfeces[is.na(poidsRfeces)] <- 0.5 poidstot[is.na(poidstot)] <- 0 defec<-ifelse(poidstot>0,1,0) poidsRgerm <- poidsgerm/poidstot poidsRgerm<-as.matrix(poidsRgerm) poidsRgerm[is.na(poidsRgerm)|!is.finite(poidsRgerm)] <- 0.5 toto <- strsplit(as.character(dimnames(Y)[[1]]),split="_") toto <- matrix(unlist(toto),ncol=2,byrow=T) Spani <- as.numeric(as.factor(toto[,1])) datax=list(Y=ifelse(is.na(Y),0,Y), # replaced NA by zeros G=ifelse(is.na(Yg),0,Yg), poidstot=as.matrix((poidstot)), defec=as.matrix(defec), poidsRfeces=as.matrix((poidsRfeces)), poidsRgerm=as.matrix((poidsRgerm)), temps=temps, difftemps=diff(temps), meant=mean(temps), meant2=mean(temps^2), nanimal=nanimal, nplante=nplante, ntemps=unlist(ntemps), Ntot=as.matrix(Ntot), Spani=Spani) #-------# # Inits #-------# Nt=array(NA,dim=c(17,6,15)) Nt[,,1] <- as.matrix(Ntot) Ybis <- Y Ybis[is.na(Ybis)] <- 0 Ybis[,1,] <- NA Ygbis <- Yg Ygbis[is.na(Ygbis)] <- 0 Ygbis[,1,] <- NA Nf <- Ybis+Ygbis Nfc <- Ybis Nfg<-Ygbis aex<-matrix(-10,ncol=6,nrow=3) bex<-matrix(0,ncol=6,nrow=3) cex<-matrix(0,ncol=6,nrow=3) ranefex=array(0,dim=c(11,6,15)) ranefex[,,1]<-NA inits=list(aex=aex, bex=bex, cex=cex, agerm=aex, bgerm=bex, Nt=Nt, Nf=Nf, Nfc=Nfc, Nfg=Nfg, sigexp=rep(0.1,3), sigeps=rep(0.1,3), siggps=rep(0.1,3) ) #-------# # Model #-------# debmodini=Sys.time() model=jags.model(file="model.r",data=datax,inits=inits,n.chains=3,n.adapt=15000,quiet=FALSE) finmodini=Sys.time() ; print(finmodini-debmodini) debmodup=Sys.time() update(model,n.iter=1000000) finmodup=Sys.time() ; print(finmodup-debmodup) debmodm=Sys.time() chain2jag=coda.samples(model,c("adef","bdef","aex","bex","cex","ranefex","pexeff","agerm","bgerm","pg","Nf","Y.rep","G.rep"),n.iter=1000000,thin=500) finmodm=Sys.time() ; print(finmodm-debmodm) save(chain2jag,file="CODASmodel.R")