Vignette to reproduce the graph in ‘How Social Learning Shapes the Efficacy of Preventative Health Behaviors in an Outbreak’

Main Paper

Fig2

This Figure explore the parameters of the sigmoid function sig:

\[ p \sim \frac{1}{1+e^{-\kappa\times(x-\nu)}} \]

which dictates the probability of individuals to switch behaviors given the percentage of infected people around them.

par(mar=c(4,4,1,1))
x=seq(0,1,.01)
inp=seq(0,1,.05)
clrs=colorRampPalette(c("purple4","yellow"))(length(inp))
plot(x,sig(x),type="n",ylim=c(0,1),xlim=c(0,1),xlab="% infected",ylab=expression(P[switch]),main="" )
for(i in 1:length(inp)) lines(x,sig(x,b=inp[i],a=4),ylim=c(0,1),xlim=c(0,1),col=clrs[i],lwd=2)
leg=seq(1,length(inp),length.out=5)
legend("bottomright",legend=rev(sapply(inp[leg],function(d)as.expression(bquote(.(d))))),title=expression(nu),col=rev(clrs[leg]),lwd=2,cex=.8,bg="white")

tstp=25
stp=rev(seq(-3,3,length.out=tstp))
clrs=colorRampPalette(c("purple4","yellow"))(tstp)
plot(x,sig(x),type="n",ylim=c(0,1),xlim=c(0,1),xlab="% infected",ylab=expression(P[switch]),main="" )
for(i in rev(1:length(stp))) lines(x,sig(x,a=10^stp[i]),ylim=c(0,1),xlim=c(0,1),col=clrs[i],lwd=2)
leg=seq(1,tstp,length.out=5)
legend("bottomright",legend=sapply(round(stp[leg]),function(d)as.expression(bquote(.(10^d)))),col=clrs[leg],lwd=2,title=expression(kappa),cex=.8,bg="white")

illustration of switch parameters

Fig3.tiff

A figure to illustrate different way to flatten the curve. We ran several simulations in two different setups to show the differences between a situation in which there is no social distancing and one in which there is.

Static simulations

The simulations have been run using:

library(parallel)
xsize=ysize=100
pop=generatePopulation(N=500,xsize=xsize,ysize=ysize,recovery=c(8,14)*25,speed=c(1,.2),behavior=rep(B,500))
cl <- makeForkCluster(18,outfile="")
neutralNSD=parSapply(cl,1:500,function(i){
                     print(i);
                     simu=slsirSimu(tstep=1500,p=c(1,1),di=2,i0=1,xsize=xsize,ysize=ysize,
                                    pop=pop,
                                    inf=9,
                                    sat=1000,
                                    p_i=1,
                                    strategy="best",
                                    ts=T,ap=F,visu=F
                                    )$timeseries[,2]

})
save(file="neutralNSD.bin",neutralNSD)

pop[,"behavior"]=G
neutralSD=parSapply(cl,1:500,function(i){
                    print(i);
                    simu=slsirSimu(tstep=1500,p=c(.1,.1),di=2,i0=1,xsize=xsize,ysize=ysize,
                                   pop=pop,
                                   inf=9,
                                   sat=1000,
                                   p_i=1,
                                   strategy="best",
                                   ts=T,ap=F,visu=F
                                   )$timeseries[,2]

})
save(file="neutralSD.bin",neutralSD)


neutralSD2=parSapply(cl,1:500,function(i){
                     print(i);
                     simu=slsirSimu(tstep=1500,p=c(.1,.1),di=2,i0=1,xsize=xsize,ysize=ysize,
                                    pop=pop,
                                    inf=9,
                                    sat=1000,
                                    p_i=1,
                                    strategy="random",
                                    ts=T,ap=F,visu=F
                                    )$timeseries[,2]

})
save(file="neutralSD2.bin",neutralSD2)
stopCluster(cl)

plot figure

par(mar=c(2,2,2,2))
plot(1,1,type="n",ylim=c(0,500),xlim=c(1,1500),ylab="",xlab="",xaxt="n",yaxt="n")
legend("topright",legend=c("No social distancing","Social distancing"),fill=c(myred,mygreen),cex=.95)
mtext("time",1,1)
mtext("number of infected people",2,1)

load("../data/neutralSD.bin")
load("../data/neutralNSD.bin")
load("../data/neutralSD2.bin")
#we use precomputed data to compute a mean curve, those data have been generated using the scripts in "exec/graphPaper.R"

meanNSD=apply(neutralNSD,1,mean)
meanSD=apply(neutralSD,1,mean)

polygon(0:1500,meanSD,col=adjustcolor(mygreen,.5),border=NA)
polygon(0:1500,meanNSD,col=adjustcolor(myred,.5),border=NA)
lines(meanSD,col=adjustcolor(mygreen,.9),lwd=2)
lines(meanNSD,col=adjustcolor(myred,.9),lwd=2)

abline(h=max(meanNSD),lty=2,col=myred)
abline(v=which.max(meanNSD),lty=2,col=myred)
abline(h=max(meanSD),lty=2,col=mygreen)
abline(v=which.max(meanSD),lty=2,col=mygreen)

text(bquote(paste("max infected (",I[max],") no SD")),x=1500,y=max(meanNSD)+10,pos=2,cex=.8)
mtext(bquote(paste("time max (",tau,") no SD")),at=which.max(meanNSD),adj=1,cex=.8)
text(bquote(paste("max infected (",I[max],") SD")),x=1500,y=max(meanSD)+10,pos=2,cex=.8)
mtext(bquote(paste("time max (",tau,") SD")),at=which.max(meanSD),adj=0,cex=.8)

abline(v=which.max(meanNSD),lty=2,col=myred)
abline(h=max(meanSD),lty=2,col=mygreen)
abline(v=which.max(meanSD),lty=2,col=mygreen)

Fig4.tiff

In this figure we present the overall results of the main simulations, with some some analysis of \(\delta\), the score of each simulation.

Run all simulations

To run all simulations one need to execute the script exec/parallelRuns.R.

Rscript exec/parallelRuns.R 20 10000 midCurveAllBadBestSLS

This will run 10000 simulations on 20 parallel cores, with priors and parameters defined in the script (cf the chunk below for the detail about the cod in the script).

We run this script multiple time on a cluster, until having enough simulation. This represent

old <- Sys.time() 

args=commandArgs(trailingOnly = TRUE) #pass number of slave by command line, should be #node - 1 as one node is used by the master 

ns=args[1]#first argument is the number of slave
nsm=args[2]#second argument the number of simulation 
mainfold=args[3] #third argument = name of the folder where to store the results

if(is.na(mainfold) | mainfold=="") mainfold=Sys.info()['nodename']

fi=0
fold=paste0(mainfold,fi)
while(file.exists(fold)){
    fi=fi+1
    fold=paste0(mainfold,fi)
}

print(paste0("ABC results will be stored in folder: ",fold))
dir.create(fold)


source("./R/slsirSimu.R")
library(parallel)
xsize=ysize=100

nsim=nsm

### Prior definition
inf=runif(nsim,0,1)
inf_r=runif(nsim,0,1)
sat=10^runif(nsim,-1,3)
sat_r=10^runif(nsim,-1,3)
pind=runif(nsim)
sl_rad=sample(100,nsim,replace=T)

allparameter=cbind.data.frame(inf=inf,sat=sat,inf_r=inf_r,sat_r=sat_r,pind=pind,sl_rad=sl_rad)


pg=0
behave=rep(c(G,B),500*c(pg,1-pg))

cl <- makeForkCluster(ns,outfile="")
allsummary=parSapply(cl,1:nsim,function(j)
                     {
                         print(paste("sim #",j,"/",nsim));
                         pop=generatePopulation(N=500,xsize=xsize,ysize=ysize,recovery=c(8,14)*25,speed=c(1,.2),behavior=behave)
                         simu=abmSIR(1500,p=c(1,.1),di=2,i0=1,recovery=c(8,14)*25,speed=c(1,.2),xsize=xsize,ysize=ysize,
                                     pop=pop,
                                     inf=allparameter$inf[j],
                                     sat=allparameter$sat[j],
                                     inf_r=allparameter$inf_r[j],
                                     sat_r=allparameter$sat_r[j],
                                     p_i=allparameter$pind[j],
                                     sl_rad=allparameter$sl_rad[j],
                                     strategy="best",
                                     ts=T,ap=F,visu=F,bt=0
                                     )
                         #save(file=file.path(fold,paste0("simu_",j,".bin")),simu)
                         max_infect=max(simu$timeseries[,2])
                         max_infect150=max(simu$timeseries[1:150,2])
                         max_infect250=max(simu$timeseries[1:250,2])
                         c(
                           max_infect=max_infect,
                           max_infect150=max_infect150,
                           max_infect250=max_infect250,
                           time_max=which.max(simu$timeseries[,2]),
                           time_max2=which.max(simu$timeseries[,2]>=(max_infect/2)),
                           time_max4=which.max(simu$timeseries[,2]>=(max_infect/4)),
                           time_max150=which.max(simu$timeseries[,2]>=(max_infect150)),
                           time_max250=which.max(simu$timeseries[,2]>=(max_infect250)),
                           final_size=sum(simu$timeseries[1500,2:3]),
                           size150=sum(simu$timeseries[150,2:3]),
                           size250=sum(simu$timeseries[250,2:3])
                           )
                     }
)


stopCluster(cl)

allresults=cbind(allparameter,t(allsummary))
save(file=file.path(fold,"allresults.bin"),allresults)
new <- Sys.time() - old # calculate difference
print(new) # print in nice format

To get the results of the simulation that have been run from the script you will need to group all file together, being sure that parameters and results stay do to so you can do:

name='../data/midCurveAllBadBestSLS/' #a folder where all results are stored
aa=lapply(unlist(lapply(list.dirs(name),list.files,pattern="allresults.bin",full.names=T)),function(f){load(f);return(allresults)})
allresults=do.call("rbind",aa)
if(length(which(allresults$max_infect <4))>0)
    allresults=allresults[-which(allresults$max_infect <4),] #we remove the results where no outbreak happen
allresults$distances=(1-bdistance(allresults$time_max) + bdistance(allresults$max_infect))/2

The resulting object is available here, if you copy it in the folder data one can load it using:

allresults=readRDS("../data/mainmodel.rds")

For this Figure we also want to compare how simulations with different score relates to the case with no social learning. For that we re-run simulations with no social learning. These simulations also take time but less space and are thus available in data/

Re-run best/worst results

library(parallel)
cl <- makeForkCluster(16,outfile="")

best=allresults[order(allresults$distances),][1:1000,]
worst=allresults[order(allresults$distances,decreasing=T),][1:1000,]

repetBest=parSapply(cl,1:500,function(i){
                    j=sample(nrow(best),1);print(paste(i,j));
                    simu=slsirSimu(500,1500,p=c(1,.1),di=2,i0=1,recovery=c(8,14)*25,speed=c(1,.2),xsize=xsize,ysize=ysize,
                                inf=best$inf[j],
                                sat=best$sat[j],
                                inf_r=best$inf_r[j],
                                sat_r=best$sat_r[j],
                                p_i=best$pind[j],
                                sl_rad=best$sl_rad[j],
                                strategy="best",
                                ts=T,ap=F,visu=F
                                )$timeseries[,2]
})
save(file="repetBest.bin",repetBest)

repetWorst=parSapply(cl,1:500,function(i){
                    j=sample(nrow(worst),1);print(paste(i,j));
                    simu=slsirSimu(500,1500,p=c(1,.1),di=2,i0=1,recovery=c(8,14)*25,speed=c(1,.2),xsize=xsize,ysize=ysize,
                                inf=worst$inf[j],
                                sat=worst$sat[j],
                                inf_r=worst$inf_r[j],
                                sat_r=worst$sat_r[j],
                                p_i=worst$pind[j],
                                sl_rad=worst$sl_rad[j],
                                strategy="best",
                                ts=T,ap=F,visu=F
                                )$timeseries[,2]
})
save(file="repetWorst.bin",repetWorst)

stopCluster(cl)

load("../data/repetBest.bin") 
load("../data/repetWorst.bin") 
meanWorst=apply(repetWorst,1,mean) #
meanBest=apply(repetBest,1,mean) #

generate Figure

We can now plot the graph of the paper

### left panel
par(mar=c(4,4,2,2))
plot(1,1,type="n",ylim=c(0,500),xlim=c(1,1500),ylab="",xlab="",xaxt="n",yaxt="n")
abline(v=150,lwd=2,col="gray")
par(xpd=NA)
text(150,-20,"t=150",col="1",pos=1)
par(xpd=F)
legend("topright",legend=c("No Learning (all NA)","No Learning (all A)","Worst Simulations","Best Simulations"),col=c(myred,mygreen),lty=c(1,1,2,2),lwd=2,cex=.95)
axis(1)
axis(2)
mtext("time",1,3)
mtext("number of infected people",2,3)


lines(meanSD,col=adjustcolor(mygreen,.9),lwd=2)
lines(meanNSD,col=adjustcolor(myred,.9),lwd=2)
lines(meanBest,col=adjustcolor(mygreen,.9),lwd=2,lty=2)
lines(meanWorst,col=adjustcolor(myred,.9),lwd=2,lty=2)

### right panel
subresults=allresults[sample(nrow(allresults),20000),]
rank=rank(subresults$distances)
orderscol=rank
color_class=rev(brewer.pal(5,"PuBu"))
colorbest=color_class[1]
orderscol=rep(adjustcolor(color_class[5],.4),length(orderscol))
nvl=c(1000,2500,5000,10000)
for(i in (length(nvl):1)) orderscol[rank<=nvl[i]]=adjustcolor(color_class[i],.4)
#plot(subresults$time_max,subresults$max_infect,bg=orderscol,main=expression(I[max]*" wrt "*tau),pch=21,xlab=expression(tau),ylab=expression(I[max]),lwd=.1)
plot(subresults$time_max,subresults$max_infect,bg=orderscol,main="",pch=21,xlab=expression(tau),ylab=expression(I[max]),lwd=.1)
legend("topright",legend=c(paste("<",nvl),"all"),pt.bg=color_class,pch=21,col=adjustcolor(1,.5),pt.lwd=.1,title="Rank")

Fig 4: mean trajectories and scores for simulations

Fig5

For this figure we need the posterior distribution of our model given by the distance to the data as measured by our \(\delta\)

Get the posterior

n=1000
posterior=list(
               best=allresults[order(allresults$distances),][1:n,],
               worst=allresults[order(allresults$distances,decreasing=T),][1:n,],
               time_max=allresults[order(allresults$time_max,decreasing=T),][1:n,],
               time_max150=allresults[order(allresults$time_max150,decreasing=T),][1:n,],
               time_max250=allresults[order(allresults$time_max250,decreasing=T),][1:n,],
               wtime_max150=allresults[order(allresults$time_max150,decreasing=F),][1:n,],
               wtime_max250=allresults[order(allresults$time_max250,decreasing=F),][1:n,],
               max_infect=allresults[order(allresults$max_infect,decreasing=F),][1:n,],
               max_infect150=allresults[order(allresults$max_infect150,decreasing=F),][1:n,],
               max_infect250=allresults[order(allresults$max_infect250,decreasing=F),][1:n,],
               wmax_infect150=allresults[order(allresults$max_infect150,decreasing=T),][1:n,],
               wmax_infect250=allresults[order(allresults$max_infect250,decreasing=T),][1:n,]
               )

library(hdrcde)

duocolrs=adjustcolor(c(myred,colorbest),.8)
currange=range=list(dim1=c(0,1),dim2=c(-1,3))
dimlab=list(dim1=expression(nu),dim2=expression(kappa))

expnames=c(distribA="Worst 150",distribB="Best 150")
marginAndJoin(
              distribA=list(dim1=posterior$wmax_infect150$inf,dim2=posterior$wmax_infect150$sat),
              distribB=list(dim1=posterior$max_infect150$inf,dim2=posterior$max_infect150$sat),
              dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
              expnames=expnames, range=currange,cols=duocolrs,log="y",
              main=""
              )

expnames=c(distribA="Worst 150",distribB="Best 150")
marginAndJoin(
              distribA=list(dim1=posterior$wmax_infect150$inf_r,dim2=posterior$wmax_infect150$sat_r),
              distribB=list(dim1=posterior$max_infect150$inf_r,dim2=posterior$max_infect150$sat_r),
              dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
              expnames=expnames, range=currange,cols=duocolrs,log="y",
              main=""
              )

expnames=c(distribA="Worst",distribB="Best")
marginAndJoin(
              distribA=list(dim1=posterior$worst$inf,dim2=posterior$worst$sat),
              distribB=list(dim1=posterior$best$inf,dim2=posterior$best$sat),
              dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
              expnames=expnames, range=currange,cols=duocolrs,log="y",
              main=""
              )


expnames=c(distribA="Worst",distribB="Best")
marginAndJoin(
              distribA=list(dim1=posterior$worst$inf_r,dim2=posterior$worst$sat_r),
              distribB=list(dim1=posterior$best$inf_r,dim2=posterior$best$sat_r),
              dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
              expnames=expnames, range=currange,cols=duocolrs,log="y",
              main=""
              )

Top row: distribution of parameters yielding the best results when the simulation reach 150 steps. Bottom row: best results when the simulation reaches 1500. Left column: parameters to swtich from NA to A, right column to switch from NA to A

Supplementary Material

FigS1.tiff

    posterior=listallposterior[["midCurveAllBadBestSLS"]]
    colors=c(P="NA",A=adjustcolor(colorbest,1),B=adjustcolor(colorbest,1))
    plot2dens(prior=(allresults$inf),A=(posterior$best$inf),B=(posterior$max_infect150$inf),cols=colors,from=0,to=1,main="posterior inflection point",xlab=expression(nu),xaxt="n")
    plot2dens(prior=(allresults$inf_r),A=(posterior$best$inf_r),B=(posterior$max_infect150$inf_r),cols=colors,from=0,to=1,main="posterior revert inflection point",xlab=expression(nu[r]),xaxt="n")
    legend("topleft",legend=c("prior","final time","150 timesteps"),fill=c(0,colorbest,colorbest),density=c(NA,NA,20))
    plot2dens(prior=1-allresults$pind,A=1-posterior$best$pind,B=1-posterior$max_infect250$pind,cols=colors,from=0,to=1,main="posterior proba social learning",xlab="proba social learning")

    plot2dens(prior=log10(allresults$sat),A=log10(posterior$best$sat),B=log10(posterior$max_infect150$sat),cols=colors,from=-1,to=3,xlab=expression(kappa),main="",xaxt="n",log=T)
    plot2dens(prior=log10(allresults$sat_r),A=log10(posterior$best$sat_r),B=log10(posterior$time_max150$sat_r),cols=colors,from=-1,to=3,main="",xlab=expression(kappa[r]),log=T)

    plot2dens(prior=allresults$sl_rad,A=posterior$best$sl_rad,B=posterior$max_infect150$sl_rad,cols=colors,from=1,to=100,main="",xlab="radius social learning")

All posterior distributions for the 6 parameters used in our simulations. The white area represents the distribution of the parameters of all simulations (the priors), the blue area represents the distribution for the top 1000 simulations rank using δ while the shaded area represent the parameters of the 1000 simulations with the lowest number of infected people after 150 timesteps.Left column represents the parameters value to switch from Non-Adherent to Adherent, middle column from Adherent to Non-Adherent, and right column the paramaters that defines the probability and the radius of social learning.

FigS2.tiff

for(d in nvl){
    best=allresults[rank(allresults$distances)<d,]
    par(mar=c(4,4,1,1))
    hdr.boxplot.2d(best$inf_r,log10(best$sat_r),prob=seq(20,100,10),shadecols=adjustcolor(colorbest,.9),xlim=c(0,1),ylim=c(-1,3),xlab=expression(nu[r]),ylab=expression(kappa[r]))
    hdr.boxplot.2d(best$inf,log10(best$sat),prob=seq(20,100,10),shadecols=adjustcolor(colorbest,.9),xlim=c(0,1),ylim=c(-1,3),xlab=expression(nu[r]),ylab=expression(kappa))
    hdr.boxplot.2d(best$inf,best$pind,prob=seq(20,100,10),shadecols=adjustcolor(colorbest,.9),xlim=c(0,1),ylim=c(0,1),xlab=expression(nu[r]),ylab="p")
    hdr.boxplot.2d(log10(best$sat),best$pind,prob=seq(20,100,10),shadecols=adjustcolor(colorbest,.9),xlim=c(-1,3),ylim=c(0,1),xlab=expression(kappa),ylab="p")
    hdr.boxplot.2d(best$sl_rad,best$pind,prob=seq(20,100,10),shadecols=adjustcolor(colorbest,.9),xlim=c(1,100),ylim=c(0,1),xlab="r",ylab="p")
    hdr.boxplot.2d(best$sl_rad,best$inf,prob=seq(20,100,10),shadecols=adjustcolor(colorbest,.9),xlim=c(1,100),ylim=c(0,1),xlab="r",ylab=expression(nu))
}

Joint posteriors of representative pair of parameters of the model for different level of what we consider as the “best” simulations. Each distribution represent the parameter distribution of simulations for which the metric δ is ranked below different level, from top to bottom: 10, 000; 5, 000; 2, 500; 1, 000. The 2d areas represent the High Density Regions, ie the smallest regions within which falls a certain percentage of the distribution of the parameters our selected simulation. The lighter areas represent the area within which all the simulations are distributed and darker areas represent regions for smaller HDR

Figure in S1_text.pdf

In S1_text we explore multiples variants of the model.

Load data from all model:

library(RColorBrewer)
modelnames=c("midCurveAllBadBestSLS","midCurve15Good","midCurve15GoodBestSLS","midCurve100GoodBestSLS","midCurveAllBad","burnin100BestSLS","testPA0.2","testPA0.6")
names(modelnames)=modelnames
modelcol=brewer.pal(length(modelnames),name="Dark2")
names(modelcol)=modelnames
modelsubnames=sapply(modelnames,function(mod)paste0(strsplit(mod,"")[[1]][-(1:8)],collapse=""))
modelsubnames=LETTERS[1:8]

setwd("../data/")

allexpe=lapply(modelnames,function(name){
               aa=lapply(unlist(lapply(list.dirs(name),list.files,pattern="allresults.bin",full.names=T)),function(f){load(f);return(allresults)})
               allresults=do.call("rbind",aa)
               if(sum(which(allresults$max_infect <4)>2))allresults=allresults[-which(allresults$max_infect <4),]
               allresults$distances=(1-bdistance(allresults$time_max) + bdistance(allresults$max_infect))/2
               allresults$distances2=(1-bdistance(allresults$time_max150) + bdistance(allresults$max_infect150))/2
               allresults$distances3=(1-bdistance(allresults$time_max250) + bdistance(allresults$max_infect250))/2
               return(allresults)})
nsims=sapply(allexpe,nrow)

Due to limited computational power, variants has been run an uneven number of time:

for(i in seq_along(modelsubnames)){
    cat(paste("* Variant",modelsubnames[i],":", nsims[i],"simulations\n\r "))
}

• Variant A : 826892 simulations

• Variant B : 153499 simulations

• Variant C : 158400 simulations

• Variant D : 251000 simulations

• Variant E : 130199 simulations

• Variant F : 309998 simulations

• Variant G : 84999 simulations

• Variant H : 84998 simulations

We can compute and compare the score of these different variants:

allscores=c("distances","distances2","distances3","time_max","time_max150","time_max250","max_infect","max_infect150","max_infect250")
names(allscores)=allscores
tablesscores=sapply(allscores,function(as){
                    allscores=lapply(allexpe,"[[",as)
                    allnames=rep(names(allscores),lengths(allscores))
                    table(factor(allnames[order(unlist(allscores),decreasing=F)][1:(10000)],levels=names(allexpe)))
})
normalized=apply(tablesscores,2,"/",nsims[rownames(tablesscores)])
dis=normalized[,1]

We can calculate the Bayes Factors for all models:

##calculate bayes factor, by row, uie mat_rat[i,j] = dis[i]/dis[j]
mat_rat=sapply(1:length(dis),function(i)sapply(1:length(dis),function(j)dis[i]/dis[j]))  
colnames(mat_rat)=names(dis)
rownames(mat_rat)=names(dis)
knitr::kable(mat_rat)

	midCurveAllBadBestSLS	midCurve15Good	midCurve15GoodBestSLS	midCurve100GoodBestSLS	midCurveAllBad	burnin100BestSLS	testPA0.2	testPA0.6
midCurveAllBadBestSLS	1.0000000	1.4168466	1.4199611	1.626095	1.4886475	1.0199993	1.0592723	0.2659396
midCurve15Good	0.7057927	1.0000000	1.0021982	1.147686	1.0506766	0.7199080	0.7476267	0.1876983
midCurve15GoodBestSLS	0.7042447	0.9978067	1.0000000	1.145168	1.0483721	0.7183290	0.7459869	0.1872866
midCurve100GoodBestSLS	0.6149704	0.8713187	0.8732340	1.000000	0.9154741	0.6272693	0.6514211	0.1635450
midCurveAllBad	0.6717507	0.9517677	0.9538598	1.092330	1.0000000	0.6851852	0.7115669	0.1786451
burnin100BestSLS	0.9803929	1.3890663	1.3921197	1.594211	1.4594594	1.0000000	1.0385030	0.2607253
testPA0.2	0.9440443	1.3375660	1.3405062	1.535105	1.4053492	0.9629245	1.0000000	0.2510588
testPA0.6	3.7602517	5.3276999	5.3394111	6.114525	5.5976894	3.8354539	3.9831305	1.0000000

Then we get the posterior distribution for all models:

n=1000
listallposterior=lapply(allexpe,function(e)
                        {
                            list(
                                 best=e[order(e$distances),][1:n,],
                                 worst=e[order(e$distances,decreasing=T),][1:n,],
                                 time_max=e[order(e$time_max,decreasing=T),][1:n,],
                                 time_max150=e[order(e$time_max150,decreasing=T),][1:n,],
                                 time_max250=e[order(e$time_max250,decreasing=T),][1:n,],
                                 max_infect=e[order(e$max_infect,decreasing=F),][1:n,],
                                 max_infect150=e[order(e$max_infect150,decreasing=F),][1:n,],
                                 max_infect250=e[order(e$max_infect250,decreasing=F),][1:n,],
                                 wtime_max=e[order(e$time_max,decreasing=F),][1:n,],
                                 wtime_max150=e[order(e$time_max150,decreasing=F),][1:n,],
                                 wtime_max250=e[order(e$time_max250,decreasing=F),][1:n,],
                                 wmax_infect=e[order(e$max_infect,decreasing=T),][1:n,],
                                 wmax_infect150=e[order(e$max_infect150,decreasing=T),][1:n,],
                                 wmax_infect250=e[order(e$max_infect250,decreasing=T),][1:n,]
                                 )
                        }
)

Distribution of the different metrics used to approximate the posteriors for all 4 variants:

par(oma=c(0,2,0,0),mar=c(4,1,1,0))
for(v in c(1,4,7)){
    var=allscores[v]
        xlimits=range(sapply(allexpe,"[[",var))
        alld=lapply(allexpe,function(i,var)density(i[[var]]),var=var)
        ylimits=range(sapply(alld,"[[","y"))
        t=""
        if(v==1)t=expression(delta)
        if(v==4)t=expression(tau)
        if(v==7)t=expression(I[max])
        plot(1,1,ylim=ylimits,xlim=xlimits,type="n",main=t,yaxt="n",xlab=t,ylab="density")
        lapply(names(alld),function(i)lines(alld[[i]],col=modelcol[i],lwd=2))

}
legend("topleft",legend=modelsubnames,col=modelcol,lwd=3,lty=1)
mtext("density",2,0,outer=T)

Distribution of the metrics described in the section Results for all simulations for all model (approx 100k simulations per model).

2D compound posteriors

duocolrs=adjustcolor(c(myred,colorbest),.8)
currange=range=list(dim1=c(0,1),dim2=c(-1,3))
dimlab=list(dim1=expression(nu),dim2=expression(kappa))

for(i in seq_along(allexpe)){

    cat(paste0("#### Model Variant ",modelsubnames[i],"\n\r"), sep="")

    posterior=listallposterior[[i]]
    worst=allexpe[[i]]

    expnames=c(distribA="Worst",distribB="Best")
    marginAndJoin(
                  distribA=list(dim1=posterior$worst$inf,dim2=posterior$worst$sat),
                  distribB=list(dim1=posterior$best$inf,dim2=posterior$best$sat),
                  dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )
    mtext(paste("Model variant",modelsubnames[i]),3,0)

    expnames=c(distribA="Worst 150",distribB="Best 150")
    marginAndJoin(
                  distribA=list(dim1=posterior$wmax_infect150$inf,dim2=posterior$wmax_infect150$sat),
                  distribB=list(dim1=posterior$max_infect150$inf,dim2=posterior$max_infect150$sat),
                  dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )
    mtext(paste("Model variant",modelsubnames[i]),3,0)

    ### Compound marginal with 2d revert learning

    expnames=c(distribA="Worst",distribB="Best")
    marginAndJoin(
                  distribA=list(dim1=posterior$worst$inf_r,dim2=posterior$worst$sat_r),
                  distribB=list(dim1=posterior$best$inf_r,dim2=posterior$best$sat_r),
                  dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )
    mtext(paste("Model variant",modelsubnames[i]),3,0)

    expnames=c(distribA="Worst 150",distribB="Best 150")
    marginAndJoin(
                  distribA=list(dim1=posterior$wmax_infect150$inf_r,dim2=posterior$wmax_infect150$sat_r),
                  distribB=list(dim1=posterior$max_infect150$inf_r,dim2=posterior$max_infect150$sat_r),
                  dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )
    mtext(paste("Model variant",modelsubnames[i]),3,0)

    ################ best vs prior

    ### Compound marginal with 2d  learning

    expnames=c(distribA="Prior",distribB="Best")
    marginAndJoin(
                  distribA=list(dim1=worst$inf,dim2=worst$sat),
                  distribB=list(dim1=posterior$best$inf,dim2=posterior$best$sat),
                  dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )

    expnames=c(distribA="Prior",distribB="Best 150")
    marginAndJoin(
                  distribA=list(dim1=worst$inf,dim2=worst$sat),
                  distribB=list(dim1=posterior$max_infect150$inf,dim2=posterior$max_infect150$sat),
                  dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )
    mtext(paste("Model variant",modelsubnames[i]),3,0)

    ### Compound marginal with 2d revert learning

    expnames=c(distribA="Prior",distribB="Best")
    marginAndJoin(
                  distribA=list(dim1=worst$inf_r,dim2=worst$sat_r),
                  distribB=list(dim1=posterior$best$inf_r,dim2=posterior$best$sat_r),
                  dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )
    mtext(paste("Model variant",modelsubnames[i]),3,0)

    expnames=c(distribA="Prior",distribB="Best 150")
    marginAndJoin(
                  distribA=list(dim1=worst$inf_r,dim2=worst$sat_r),
                  distribB=list(dim1=posterior$max_infect150$inf_r,dim2=posterior$max_infect150$sat_r),
                  dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )
    mtext(paste("Model variant",modelsubnames[i]),3,0)
    cat(paste0("\n\r----\n\r"),sep="")
}

Model Variant A

---

Model Variant B

---

Model Variant C

---

Model Variant D

---

Model Variant E

---

Model Variant F

---

Model Variant G

---

Model Variant H

---

posterior for each parameter, each models

for(i in seq_along(allexpe)){
    cat(paste0("#### Model Variant ",modelsubnames[i],"\n\r"), sep="")
    posterior=listallposterior[[i]]
    mar=c(4,3,2,1)

    colors=c(P="NA",A=adjustcolor(colorbest,1),B=adjustcolor(colorbest,1))

    #par(mar=mar)
    plot2dens(prior=1-allresults$pind,A=1-posterior$best$pind,B=1-posterior$max_infect250$pind,cols=colors,from=0,to=1,main="posterior proba social learning",xlab="proba social learning")
    plot2dens(prior=(allresults$inf_r),A=(posterior$best$inf_r),B=(posterior$max_infect150$inf_r),cols=colors,from=0,to=1,main="posterior revert inflection point",xlab=expression(nu[r]),xaxt="n")
    legend("topleft",legend=c("prior","final time","150 timesteps"),fill=c(0,colorbest,colorbest),density=c(NA,NA,20))
    plot2dens(prior=(allresults$inf),A=(posterior$best$inf),B=(posterior$max_infect150$inf),cols=colors,from=0,to=1,main="posterior inflection point",xlab=expression(nu),xaxt="n")

    plot2dens(prior=allresults$sl_rad,A=posterior$best$sl_rad,B=posterior$max_infect150$sl_rad,cols=colors,from=1,to=100,main="",xlab="radius social learning")

    plot2dens(prior=log10(allresults$sat_r),A=log10(posterior$best$sat_r),B=log10(posterior$time_max150$sat_r),cols=colors,from=-1,to=3,main="",xlab=expression(kappa[r]),log=T)

    plot2dens(prior=log10(allresults$sat),A=log10(posterior$best$sat),B=log10(posterior$max_infect150$sat),cols=colors,from=-1,to=3,xlab=expression(kappa),main="",xaxt="n",log=T)

    cat(paste0("\n\r----\n\r"),sep="")
    

}

Model Variant A

---

Model Variant B

---

Model Variant C

---

Model Variant D

---

Model Variant E

---

Model Variant F

---

Model Variant G

---

Model Variant H

---