This script was used to process raw data of Bactrocera oleae behaviour which were collected using JWatcher 1.0 and reported in Daher, E.; Cinosi, N.; Chierici, E.; Rondoni, G.; Famiani, F.; Conti, E. Field and laboratory efficacy of low-impact commercial products in preventing olive fruit fly, Bactrocera oleae, infestation. Insects 2022, 13, 213. https://doi.org/10.3390/insects13020213.

Load required packages

rm(list=ls())
library(knitr)
library(ggplot2)
library(survival)

options(max.print=100000)

Read the data and format the dataset

list<-read.csv2("list.csv", header=T,dec=".",sep=";")
dat<- read.csv2("output.csv", header=T,dec=".",sep=";")

dat$time<-0
dat$time<-sapply(strsplit(dat$time_old,":"),
  function(x) {  x <- as.numeric(x)
  x[1]*3600+x[2]*60+x[3]+x[4]/100  } )
head(dat)

lin<-which(dat$key=="-") #while recording and in case of mistyping, do press the key "-" followed by the correct key. This code will replace the wrong key with the correct one. 
if(length(lin)>0){
dat$key[lin-1]<-dat$key[lin+1]
rem<-sort(c(lin,lin+1))
dat<-dat[-rem,]
}

Assign keyboard keys to behavioural events

Relative to the 12 arenas, the 12 keys: q,w,e,r,t,y,u,i,o,p,z,x were pressed at the beginning and each time the insect was present in the external area of the arena. The 12 keys: 1,2,3,4,5,6,7,8,9,0,a,s were pressed each time the insect entered into the central area containing the olive twigs.

list$ttv<-list$tot<-list$choice<-list$no_choice<-as.numeric(0)

### change uppercase with lowercase 
dat$key<-gsub("Q", "q", dat$key)
dat$key<-gsub("W", "w", dat$key)
dat$key<-gsub("E", "e", dat$key)
dat$key<-gsub("R", "r", dat$key)
dat$key<-gsub("T", "t", dat$key)
dat$key<-gsub("Y", "y", dat$key)
dat$key<-gsub("U", "u", dat$key)
dat$key<-gsub("I", "i", dat$key)
dat$key<-gsub("O", "o", dat$key)
dat$key<-gsub("P", "p", dat$key)
dat$key<-gsub("A", "a", dat$key)
dat$key<-gsub("Z", "z", dat$key)
dat$key<-gsub("S", "s", dat$key)
dat$key<-gsub("X", "x", dat$key)

r1<-c("q","1");r2<-c("w","2");r3<-c("e","3");r4<-c("r","4");r5<-c("t","5");r6<-c("y","6");r7<-c("u","7");r8<-c("i","8");r9<-c("o","9");r10<-c("p","0"); r11<-c("z","a");r12<-c("x","s")

rep<-list(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12)

Raw data processing

exper<-unique(dat$experiment)
dat$duration<-dat$ch<-dat$tt<-dat$ti<-as.numeric(0)
for (i in 1:length(exper)){
new1<-subset(dat,dat$experiment==exper[i])  
for (j in 1:length(rep)){
new2<-new1[which(new1$key %in% rep[[j]]),]
startt<-new2$time[1] + 0   
minimum<-120 ## the first two min are left for acclimation [ref 30,31 of Daher et al., 2022]
maximum<-600 ## then, consider ten min of behavioural observations
new2$ti<-new2$time-startt
if(length(new2$time) > 1){
last<-length(which(new2$ti<=minimum))
      new2$ti[which(new2$ti<=minimum)]<-120 
      new2$tt<-new2$ti-minimum
      
  for (k in 1:(length(new2$time))-1) {
  if(max(new2$tt)>=maximum){
    new2$tt[which(new2$tt>=maximum)]<-maximum
    new2$duration[k]<- new2$tt[k+1]-new2$tt[k] 
}
if(max(new2$tt)<maximum){
 new2$duration[k]<- new2$tt[k+1]-new2$tt[k] 
 new2$duration[length(new2[,2])]<-maximum-new2$tt[length(new2[,2])]
}
     
new2$ch[which(new2$key==rep[[j]][1])]<-"0"
new2$ch[which(new2$key==rep[[j]][2])]<-"1"
new3<-new2[which(new2$ti>minimum),]
ttv<-new3$tt[which(new3$ch %in% 1)][1]
ag<-aggregate(duration ~ ch, sum, data=new2)

if(length(ag$ch)==1 && ag$ch[1]=="0"){
list$no_choice[which(list$experiment==i&list$arena==j)]<-ag[which(ag$ch==0),2]
list$tot[which(list$experiment==i&list$arena==j)]<-ag[which(ag$ch==0),2]
list$ttv[which(list$experiment==i&list$arena==j)]<-ttv  #time to visit
}

if(length(ag$ch)==2){
list$no_choice[which(list$experiment==i&list$arena==j)]<-ag[which(ag$ch==0),2]
list$choice[which(list$experiment==i&list$arena==j)]<-ag[which(ag$ch==1),2]
list$tot[which(list$experiment==i&list$arena==j)]<-ag[which(ag$ch==0),2]+ag[which(ag$ch==1),2]
list$ttv[which(list$experiment==i&list$arena==j)]<-ttv #time to visit
}
if(length(ag$ch)==1 && ag$ch[1]=="1"){
list$no_choice[which(list$experiment==i&list$arena==j)]<-0
list$choice[which(list$experiment==i&list$arena==j)]<-ag[which(ag$ch==1),2]
list$tot[which(list$experiment==i&list$arena==j)]<-ag[which(ag$ch==1),2]
list$ttv[which(list$experiment==i&list$arena==j)]<-ttv #time to visit
}}}}}
list2<-list

Inspect the new dataset

kable(head(list2,14))
experiment arena ind tr no_choice choice tot ttv
1 1 bo001 propolis 167.56 432.44 600 17.59
1 2 bo002 rock powder 600.00 0.00 600 NA
1 3 bo003 propolis + rock powder 600.00 0.00 600 NA
1 4 bo004 copper oxychloride 600.00 0.00 600 NA
1 5 bo005 copper sulphate 600.00 0.00 600 NA
1 6 bo006 control 191.12 408.88 600 8.13
1 7 bo007 propolis 379.19 220.81 600 8.83
1 8 bo008 rock powder 600.00 0.00 600 NA
1 9 bo009 propolis + rock powder 219.67 380.33 600 9.55
1 10 bo010 copper oxychloride 600.00 0.00 600 NA
1 11 bo011 copper sulphate 600.00 0.00 600 NA
1 12 bo012 control 112.73 487.27 600 9.79
2 1 bo013 propolis 600.00 0.00 600 NA
2 2 bo014 rock powder 147.51 452.49 600 147.51

Figure 1: Box plot for total residence time

list2<-list
mylab<-c("control", "copper\noxychloride","copper\nsulphate", "propolis","propolis +\nrock powder","rock\npowder")
bp <- ggplot(list2, aes(x = tr, y = choice, fill = "grey20"))
bp <- bp+scale_x_discrete(labels=mylab)
bp <- bp + geom_boxplot(notch = F, outlier.shape=1, width=0.5)
bp <- bp + ylim(c(0,610))
my_colors <- c("#999999")
bp <- bp + scale_fill_manual(values=my_colors)
bp <- bp + labs(x=" ", y = "Residence time (s)")
bp <- bp + theme_bw()
bp <- bp+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
bp <- bp + theme(legend.position="none")
bp <- bp + theme(axis.text.x = element_text(colour="black",size=13,angle=0,hjust=.5,vjust=.5,face="plain"),
axis.text.y = element_text(colour = "black",size=11,angle=0,hjust=1,vjust=.4,face="plain"),
axis.title.y = element_text(colour="black",size=15,angle=90,hjust=.5,vjust=2.5,face="plain"),
plot.title = element_text(size=18, face = "bold", hjust = 0.5))
bp

Figure 2: Cumulative hazards for B. oleae time to visit

list2<-list
list2$status<-1
list2$status[which(list2$ttv %in% NA)]<-0
list2$ttv[which(list2$ttv %in% NA)]<-maximum
cm<-coxph(Surv(ttv,status)~I(tr),data=list2)
aa1<-survfit(cm,data.frame(tr=c("control","propolis","propolis + rock powder","rock powder","copper oxychloride","copper sulphate")),fun="cumhaz")
plot(aa1,fun="cumhaz", xlab="Time (s)",ylab="Visit probability",lwd=2,lty=c(1,1,1,1,1),las=1)