#Authors: N.G. Yoccoz and N. Lecomte
# simulations with TurnoverSim
#a simulation tool for stable isotope analyses to test whether experiments were able to achieve the #requirements for measuring isotopic incorportation into tissues
#available at ‘www.arctic-predators.uit.no’. 

nind=200
# Number of individuals
ndays=200
 # Number of days
nt = nind*ndays
# Number of total data points
days=rep(1:ndays,times=nind)
str(days)
# days of sampling
id = factor(rep(1:nind,each = ndays))  # id of individuals
n.sim=200
mat.fix0=mat.fix1=mat.fix2=mat.fix3=mat.fix4=mat.fix5=array(0,dim=c(n.sim,3))
sd.res=0.2
mean.int= -24.2            # mean value of the yinit for each individual
mean.a= -4.2               # mean value of the slope
mean.asymp= -19.3          # mean value of the asymptotic value
days.sel0=c(1,10,20,30,40,50,70,90,110,150,200)
days.sel1=c(1,10,20,30,40,50,70,90,110,150)
days.sel2=c(1,10,20,30,40,50,70,90,110)
days.sel3=c(1,10,20,30,40,50,70,90)
days.sel4=c(1,10,20,30,40,50,70)
days.sel5=c(1,10,20,30,40,50)

for (i in 1:n.sim) {
intercept= rep(rnorm(nind, mean=mean.int, sd=0.2),each=ndays)
# simulated value of the yinit for each individual
a=rep(rnorm(nind, mean=mean.a, sd=0.01),each=ndays)
# simulated value of the slope
asymp=rep(rnorm(nind, mean=mean.asymp, sd=0.02),each=ndays)
# simulated value of the asymptotic value
d13c = asymp+(intercept-asymp)*exp(-exp(a)*(1:ndays)) + rnorm(n=nt, 0 , sd=sd.res)
aas=as.data.frame(cbind(days,id,d13c))
mat.fix0[i,]=fixed.effects(nlme(d13c ~ SSasymp(days,Asym,R0,lrc), fixed=Asym+R0+lrc~1, random=R0~1|id, start=c(Asym=-15,R0=-25,lrc=-4),data=aas[which(aas$days %in% days.sel0),]))
mat.fix1[i,]=fixed.effects(nlme(d13c ~ SSasymp(days,Asym,R0,lrc), fixed=Asym+R0+lrc~1, random=R0~1|id, start=c(Asym=-15,R0=-25,lrc=-4),data=aas[which(aas$days %in% days.sel1),]))
mat.fix2[i,]=fixed.effects(nlme(d13c ~ SSasymp(days,Asym,R0,lrc), fixed=Asym+R0+lrc~1, random=R0~1|id, start=c(Asym=-15,R0=-25,lrc=-4),data=aas[which(aas$days %in% days.sel2),]))
mat.fix3[i,]=fixed.effects(nlme(d13c ~ SSasymp(days,Asym,R0,lrc), fixed=Asym+R0+lrc~1, random=R0~1|id, start=c(Asym=-15,R0=-25,lrc=-4),data=aas[which(aas$days %in% days.sel3),]))
mat.fix4[i,]=fixed.effects(nlme(d13c ~ SSasymp(days,Asym,R0,lrc), fixed=Asym+R0+lrc~1, random=R0~1|id, start=c(Asym=-15,R0=-25,lrc=-4),data=aas[which(aas$days %in% days.sel4),]))
mat.fix5[i,]=fixed.effects(nlme(d13c ~ SSasymp(days,Asym,R0,lrc), fixed=Asym+R0+lrc~1, random=R0~1|id, start=c(Asym=-15,R0=-25,lrc=-4),data=aas[which(aas$days %in% days.sel5),]))
}

apply(mat.fix1,MAR=2,mean)
apply(mat.fix1,MAR=2,sd)
apply(mat.fix2,MAR=2,mean)
apply(mat.fix2,MAR=2,sd)
apply(mat.fix3,MAR=2,mean)
apply(mat.fix3,MAR=2,sd)

par(mfrow=c(3,1))
boxplot(mat.fix0[,1],mat.fix1[,1],mat.fix2[,1],mat.fix3[,1],mat.fix4[,1],mat.fix5[,1],names=c("200","150","110","90","70","50"),las=1)
abline(h=mean.asymp); title("Asymptote")
boxplot(mat.fix0[,2],mat.fix1[,2],mat.fix2[,2],mat.fix3[,2],mat.fix4[,2],mat.fix5[,2],names=c("200","150","110","90","70","50"),las=1)
abline(h=mean.int); title("Intercept")
boxplot(log(2)/exp(mat.fix0[,3]),log(2)/exp(mat.fix1[,3]),log(2)/exp(mat.fix2[,3]),log(2)/exp(mat.fix3[,3]),log(2)/exp(mat.fix4[,3]),log(2)/exp(mat.fix5[,3]),
names=c("200","150","110","90","70","50"),las=1)
abline(h=log(2)/exp(mean.a)); title("Half-life")
#end of the simulations. Further updates of the scripts would be available at www.arctic-predators.uit.no