library(circular) library(insol) library(pbapply) library(mgcv) library(activity) # The activity package has been updated to include the necesary functions # Update to the latest version (1.3). # 1. Underlying distribution functions ---- # von mises mixture distribution pdf # prm is named vector of parameters: # k1, k2... mixture conentrations # m1, m2... mixture positions in radians # p2... mixture proportions (p1 is implicitly 1-sum(p2,...)) dvmsmix <- function(x,prm){ k <- prm[grep("k",names(prm))] m <- prm[grep("m",names(prm))] p <- prm[grep("p",names(prm))] if(length(k)==1) p <- 1 else p <- c(1-sum(p), p) res <- rep(0, length(x)) for(i in 1:length(k)) res <- res + p[i]*dvonmises(circular(x),circular(m[i]),k[i]) res } #random numbers sampled from von mises mixture distribution, prm as above rvmsmix <- function(n,prm){ k <- prm[grep("k",names(prm))] m <- prm[grep("m",names(prm))] p <- prm[grep("p",names(prm))] if(length(k)==1) p <- 1 else p <- c(1-sum(p), p) nn <- as.vector(rmultinom(1,n,p)) res <- vector(length=0) for(i in 1:length(nn)){ if(nn[i]>0) res <- c(res,as.numeric(rvonmises(nn[i],circular(m[i]),k[i]))) } res } #create parameter vector for von mises mixture distribution pars <- function(anchR, anchS, type=c("crep","diur","cath"), k=10, nn=20, poff=0, p=0.8){ type <- match.arg(type) if(type=="crep") prm <- c(k1=k,k2=k, m1=anchR,m2=anchS, p2=0.5) else if(type=="diur"){ relp <- seq(-p,p,len=nn)^2+p/2 prm <- c(rep(k, nn), seq(anchR, anchS, len=nn), (relp/sum(relp))[-1]) names(prm) <- c(paste0("k",1:nn), paste0("m",1:nn), paste0("p",2:nn)) } else if(type=="cath"){ sqq <- seq(2*pi/nn, 2*pi, len=nn) prm <- c(k, 2*(cos(sqq)+1.1), anchS, sqq, rep(p/nn, nn)) names(prm) <- c(paste0("k",1:(nn+1)), paste0("m",1:(nn+1)), paste0("p",2:(nn+1))) } prm } #generate random values from VMmix distribution given sun rise and set vectors rtimes <- function(rise, set, type=c("crep", "diur", "cath"), k=10, nn=20, poff=0, p=0.8){ type <- match.arg(type) f <- function(i){ prm <- pars(rise[i], set[i], type, k=k, nn=nn, poff=poff, p=p) rvmsmix(1, prm) } sapply(1:length(rise), f) } fit.actmods <- function(dd, lat, cmo, type, k, nn, poff, p, cday=21){ tseq <- seq(0,2*pi,len=512) jd0 <- JDymd(2015,cmo,21)-round(dd/2) #julian day of start jdseq <- jd0+0:(dd-1) #sequence of julian days across study stimes <- daylength(lat,0,jdseq,0) #sun times on study days Rd <- stimes[,"sunrise"]*pi/12 #daily sunrise times Sd <- stimes[,"sunset"]*pi/12 #daily sunset times jd <- sample(jdseq, n, rep=TRUE)-jd0+1 #random julian days of events Tr <- Rd[jd] #sunrise clock times Ts <- Sd[jd] #sunset clock times Tc <- rtimes(Tr, Ts, type, k, nn, poff, p) #event clock times Tn <- transtime(Tc, Ts, type="sing") #solar time, nouvellet anchor Te <- transtime(Tc, cbind(Tr,Ts), type="eq") #solar time, equinoctial anchor Ta <- transtime(Tc, cbind(Tr,Ts)) #solar time, average anchor true <- dvmsmix(tseq, pars(mean(Rd), mean(Sd), type, k, nn, poff, p)) true.eq <- dvmsmix(tseq, pars(pi/2, pi*3/2, type, k, nn, poff, p)) c(clock=fitact(Tc,sample="n")@act, single=fitact(Tn,sample="n")@act, equinox=fitact(Te,sample="n")@act, average=fitact(Ta,sample="n")@act, true=1/(2*pi*max(true)), true.eq=1/(2*pi*max(true.eq))) } fig1 <- function(sq, true, true.eq, type, k, nn, poff, p, ylm, letter=FALSE, ylab=""){ true <- dvmsmix(sq, pars(mean(Rd), mean(Sd), type, k, nn, poff, p)) true.eq <- dvmsmix(sq, pars(pi/2, pi*3/2, type, k, nn, poff, p)) #Generate and transform some random event data, then fit models jd <- sample(jdseq, n, rep=TRUE)-jd0+1 #random julian days of events Tr <- Rd[jd] #sunrise clock times Ts <- Sd[jd] #sunset clock times Tc <- rtimes(Tr, Ts, type, k, nn, poff, p) #event clock times Ta <- transtime(Tc, cbind(Tr,Ts)) #solar time, average anchor Te <- transtime(Tc, cbind(Tr,Ts), type="eq") #solar time, equinoctial anchor Tn <- transtime(Tc, Ts, type="sing") #solar time, nouvellet anchor dat <- list(Tc,Tn,Te,Ta) mods <- lapply(dat, fitact, sample="n", adj=1.5) names <- c("Clock", "Single-anchored", "Equinoctial double-anchored", "Double-anchored") #Plot true underlying patterns plotted at regular intervals across the study, then models yxp <- 0.2 #top y axis tickmark position mgp <- c(2,0.6,0) #axis label positions nt <- 6 #number of trends to plot cols <- rainbow(nt, end=0.7)[nt:1] d <- round(seq(1,dd,len=nt)) #days at which to plot plot(range(sq*12/pi), c(0,0.5), xlab="", ylab="", las=1, type="n", xaxt="n", xaxp=c(0,24,4), ylim=c(0,ylm), yaxp=c(0,yxp,2), mgp=mgp) axis(1, 6*(0:4), labels=F) for(i in 1:nt){ prm <- pars(Rd[d[i]], Sd[d[i]], type, k, nn, poff, p) lines(sq*12/pi, dvmsmix(sq, prm)*pi/12, col=cols[i]) } if(letter) text(1,0.95*ylm,"a)") rads <- round((0:4)*pi/2, 2) for(i in 1:4){ if(i==2) ylb <- ylab else ylb <- "" if(i==4) xlb <-"Time of day" else xlb <- "" if(i==3) tr <- true.eq*pi/12 else tr <- true*pi/12 plot(mods[[i]], yunit="d", las=1, ylim=c(0,ylm), tline=list(col=1), dline=list(col="grey50"), xaxis=list(xaxt="n"), xlab=xlb, ylab=ylb, yaxp=c(0,yxp,2), mgp=mgp) if(i==4){ axis(1, 6*(0:4))#, padj=-0.8) } else axis(1, 6*(0:4), labels=F) lines(sq*12/pi, tr, col=3) if(i==3){ lines(rep(6,2), 0:1, col=2) lines(rep(18,2), 0:1, col=2) } else if(i>1){ lines(rep(mean(Sd)*12/pi,2), c(0,ylm), col=2) if(i==4) lines(rep(mean(Rd)*12/pi,2), c(0,ylm), col=2) } if(letter) text(1,0.95*ylm,paste(letters[i+1],")",sep="")) } true.act <- 1/(2*pi*max(true)) trueq.act <- 1/(2*pi*max(true.eq)) unlist(lapply(mods, function(m) m@act)) / c(true.act, true.act, trueq.act, true.act) } fig2dat <- function(sq, type, k, nn, poff, p){ res1 <- pbsapply(sq, fit.actmods, lat=60, cmo=3, type=type, k=k, nn=nn, poff=poff, p=p) res2 <- pbsapply(sq, fit.actmods, lat=60, cmo=6, type=type, k=k, nn=nn, poff=poff, p=p) res3 <- pbsapply(sq, fit.actmods, lat=40, cmo=3, type=type, k=k, nn=nn, poff=poff, p=p) res4 <- pbsapply(sq, fit.actmods, lat=40, cmo=6, type=type, k=k, nn=nn, poff=poff, p=p) res5 <- pbsapply(sq, fit.actmods, lat=20, cmo=3, type=type, k=k, nn=nn, poff=poff, p=p) res6 <- pbsapply(sq, fit.actmods, lat=20, cmo=6, type=type, k=k, nn=nn, poff=poff, p=p) mods <- list(c1 = gam(res1[1,]/res1[5,]~te(sq,k=3)), c2 = gam(res2[1,]/res2[5,]~te(sq,k=3)), c3 = gam(res3[1,]/res3[5,]~te(sq,k=3)), c4 = gam(res4[1,]/res4[5,]~te(sq,k=3)), c5 = gam(res5[1,]/res5[5,]~te(sq,k=3)), c6 = gam(res6[1,]/res6[5,]~te(sq,k=3)), s1 = gam(res1[2,]/res1[5,]~te(sq,k=3)), s2 = gam(res2[2,]/res2[5,]~te(sq,k=3)), s3 = gam(res3[2,]/res3[5,]~te(sq,k=3)), s4 = gam(res4[2,]/res4[5,]~te(sq,k=3)), s5 = gam(res5[2,]/res5[5,]~te(sq,k=3)), s6 = gam(res6[2,]/res6[5,]~te(sq,k=3)), e1 = gam(res1[3,]/res1[5,]~te(sq,k=3)), e2 = gam(res2[3,]/res2[5,]~te(sq,k=3)), e3 = gam(res3[3,]/res3[5,]~te(sq,k=3)), e4 = gam(res4[3,]/res4[5,]~te(sq,k=3)), e5 = gam(res5[3,]/res5[5,]~te(sq,k=3)), e6 = gam(res6[3,]/res6[5,]~te(sq,k=3)), a1 = gam(res1[4,]/res1[5,]~te(sq,k=3)), a2 = gam(res2[4,]/res2[5,]~te(sq,k=3)), a3 = gam(res3[4,]/res3[5,]~te(sq,k=3)), a4 = gam(res4[4,]/res4[5,]~te(sq,k=3)), a5 = gam(res5[4,]/res5[5,]~te(sq,k=3)), a6 = gam(res6[4,]/res6[5,]~te(sq,k=3)) ) gamtr <- lapply(mods, function(m) predict(m)) } fig2panel <- function(trends, ylim, xlab="", ylab="", main="", letter="", xlabels=TRUE){ plot(dsq, trends[[1]], ylim=ylim, type="n", las=1, ylab=ylab, xlab=xlab, xaxt="n", main=main) axis(1, labels=xlabels) for(j in 1:6) lines(dsq, trends[[j]], lty=(j+1) %% 2 + 1, col=(6-j) %/% 2 + 2) text(80,ylim[2], letter) } # 2. Set parameters ---- n <- 5000 #number of events to generate dd <- 90 #number of days for study lat <- 60 #latitude of study jd0 <- JDymd(2015,12,21) #julian day of start #Set up daily sun times jdseq <- jd0+0:(dd-1) #sequence of julian days across study stimes <- daylength(lat,0,jdseq,0) #sun times on study days Rd <- stimes[,"sunrise"]*pi/12 #daily sunrise times Sd <- stimes[,"sunset"]*pi/12 #daily sunset times tsq <- seq(0,2*pi,len=256) #time of day sequence dsq <- seq(60,360,60) #study length sequence # 2. a) Crepuscular ---- type1="crep" ylm1 <- 0.2 #y axis limit for plotting k1=10 #von mises concentration parameter (higher is sharper) nn1=NA #number of mixtures poff1=NA #time offfset of peaks from sun rise/set p1=NA #parameter controling depth of midday trough in diurnal / height of peak in cathermeral act.ests1 <- fit.actmods(90, 60, 2, type1, k1, nn1, poff1, p1, cday=4) trends1 <- fig2dat(dsq, type1, k1, nn1, poff1, p1) #2. b) Diurnal----- type2="diur" ylm2 <- 0.22 #y axis limit for plotting k2=100 nn2=50 poff2=0 p2=0.4 act.ests2 <- fit.actmods(90, 60, 2, type2, k2, nn2, poff2, p2, cday=4) trends2 <- fig2dat(dsq, type2, k2, nn2, poff2, p2) #2. c) Cathemeral---- type3="cath" ylm3 <- 0.12 #y axis limit for plotting k3=8 nn3=10 poff3=0 p3=0.8 act.ests3 <- fit.actmods(90, 60, 2, type3, k3, nn3, poff3, p3, cday=4) trends3 <- fig2dat(dsq, type3, k3, nn3, poff3, p3) # 3. Relative activity levels for FIG 1 ---- act.ests1["clock.act"] / act.ests1["true"] act.ests2["clock.act"] / act.ests2["true"] act.ests3["clock.act"] / act.ests3["true"] act.ests1["single.act"] / act.ests3["true"] act.ests2["single.act"] / act.ests2["true"] act.ests3["single.act"] / act.ests3["true"] act.ests1["equinox.act"] / act.ests1["true.eq"] act.ests2["equinox.act"] / act.ests2["true.eq"] act.ests3["equinox.act"] / act.ests3["true.eq"] act.ests1["average.act"] / act.ests1["true"] act.ests2["average.act"] / act.ests2["true"] act.ests3["average.act"] / act.ests3["true"] # 4. Figures ---- # 4. a) Figure 1. Activity pattern simulations---- layout(matrix(1:15, ncol=3)) par(mar=c(3,3,0,0)) fig1(tsq, true1, true.eq1, type1, k1, nn1, poff1, p1, ylm1, TRUE, "PDF") fig1(tsq, true2, true.eq2, type2, k2, nn2, poff2, p2, ylm2) fig1(tsq, true3, true.eq3, type3, k3, nn3, poff3, p3, ylm3) # 4. b) Figure 2. Relative Activity level simulations ---- layout(matrix(1:6, ncol=3)) par(mar=c(4,4,0.5,0.5), cex=0.8) ylim <- c(0.9, 3) fig2panel(trends1[grep("c", names(trends1))], ylim, ylab="Relative activity level", xlabels=FALSE, letter="a)") fig2panel(trends1[grep("a", names(trends1))], ylim, ylab="Relative activity level", letter="d)") fig2panel(trends2[grep("c", names(trends2))], ylim, xlabels=FALSE, letter="b)") fig2panel(trends2[grep("a", names(trends2))], ylim, xlab="Study length (days)", letter="e)") fig2panel(trends3[grep("c", names(trends3))], ylim, xlabels=FALSE, letter="c)") fig2panel(trends3[grep("a", names(trends3))], ylim, letter="f)") legend(90,2.8, c("60", "40", "20"), title="Latitude", lty=1, col=4:2, cex=0.8, bty="n") legend(90,2.2, c("equinox","solstice"), title="Study centre", lty=1:2, col=1, cex=0.8, bty="n") # 5. Illustrating day length changes for different study lengths ---- par(mfrow=c(4,2)) dd <- 91 root <- as.POSIXct("1970:03:21", format="%Y:%m:%d", tz="UTC") jd0 <- JDymd(2015,3,21)-round(dd/2) #julian day of start (centring study period on given date) jdseq <- jd0+0:(dd-1) #sequence of julian days across study Date <- root+(jdseq-median(jdseq))*24*3600 stimes <- daylength(lat,0,jdseq,0) #sun times on study days plot(Date, stimes[,"daylen"], ylab="Day length", type="l", main=paste("Monotonic", dd, "days")) root <- as.POSIXct("1970:06:21", format="%Y:%m:%d", tz="UTC") jd0 <- JDymd(2015,6,21)-round(dd/2) #julian day of start (centring study period on given date) jdseq <- jd0+0:(dd-1) #sequence of julian days across study Date <- root+(jdseq-median(jdseq))*24*3600 stimes <- daylength(lat,0,jdseq,0) #sun times on study days plot(Date, stimes[,"daylen"], ylab="Day length", type="l", main=paste("Symmetric", dd, "days")) dd <- 182 root <- as.POSIXct("1970:03:21", format="%Y:%m:%d", tz="UTC") jd0 <- JDymd(2015,3,21)-round(dd/2) #julian day of start (centring study period on given date) jdseq <- jd0+0:(dd-1) #sequence of julian days across study Date <- root+(jdseq-median(jdseq))*24*3600 stimes <- daylength(lat,0,jdseq,0) #sun times on study days plot(Date, stimes[,"daylen"], ylab="Day length", type="l", main=paste("Monotonic", dd, "days")) root <- as.POSIXct("1970:06:21", format="%Y:%m:%d", tz="UTC") jd0 <- JDymd(2015,6,21)-round(dd/2) #julian day of start (centring study period on given date) jdseq <- jd0+0:(dd-1) #sequence of julian days across study Date <- root+(jdseq-median(jdseq))*24*3600 stimes <- daylength(lat,0,jdseq,0) #sun times on study days plot(Date, stimes[,"daylen"], ylab="Day length", type="l", main=paste("Symmetric", dd, "days")) dd <- 273 root <- as.POSIXct("1970:03:21", format="%Y:%m:%d", tz="UTC") jd0 <- JDymd(2015,3,21)-round(dd/2) #julian day of start (centring study period on given date) jdseq <- jd0+0:(dd-1) #sequence of julian days across study Date <- root+(jdseq-median(jdseq))*24*3600 stimes <- daylength(lat,0,jdseq,0) #sun times on study days plot(Date, stimes[,"daylen"], ylab="Day length", type="l", main=paste("Monotonic", dd, "days")) root <- as.POSIXct("1970:06:21", format="%Y:%m:%d", tz="UTC") jd0 <- JDymd(2015,6,21)-round(dd/2) #julian day of start (centring study period on given date) jdseq <- jd0+0:(dd-1) #sequence of julian days across study Date <- root+(jdseq-median(jdseq))*24*3600 stimes <- daylength(lat,0,jdseq,0) #sun times on study days plot(Date, stimes[,"daylen"], ylab="Day length", type="l", main=paste("Symmetric", dd, "days")) dd <- 365 root <- as.POSIXct("1970:03:21", format="%Y:%m:%d", tz="UTC") jd0 <- JDymd(2015,3,21)-round(dd/2) #julian day of start (centring study period on given date) jdseq <- jd0+0:(dd-1) #sequence of julian days across study Date <- root+(jdseq-median(jdseq))*24*3600 stimes <- daylength(lat,0,jdseq,0) #sun times on study days plot(Date, stimes[,"daylen"], ylab="Day length", type="l", main=paste("Monotonic", dd, "days")) root <- as.POSIXct("1970:06:21", format="%Y:%m:%d", tz="UTC") jd0 <- JDymd(2015,6,21)-round(dd/2) #julian day of start (centring study period on given date) jdseq <- jd0+0:(dd-1) #sequence of julian days across study Date <- root+(jdseq-median(jdseq))*24*3600 stimes <- daylength(lat,0,jdseq,0) #sun times on study days plot(Date, stimes[,"daylen"], ylab="Day length", type="l", main=paste("Symmetric", dd, "days"))