#Load the code for running NBDA from https://lalandlab.st-andrews.ac.uk/freeware/ #function for calculating estimated proportion of events occurring by social learning oadaPropSolveByST<-function(par,oadata,innovations=1,additive=T,asocialVar=NULL){ if(is.character(oadata)) oadata<- combineOaCoxData(oadata); s<-par[1] beta<-par[-1] if(is.null(asocialVar)|!additive){ data2<-oadata@coxdata[oadata@coxdata$status==1,6] number<-sum((s * data2)/(1+ s * data2)) total<-length(data2)-innovations }else{ data2<-oadata@coxdata[oadata@coxdata$status==1,c(6,6+asocialVar)] number<-sum((s * data2[,1])/(exp(as.matrix(data2[,-1])%*%beta)+ s * data2[,1])) total<-dim(data2)[1]-innovations } return(number/total) } #Set number of simulations noSims<-1000 #Set number of individuals noInd<-100 #Set s values to be considered sVect<-seq(0,4,0.5) #Set noise to be considered NoiseVect<-c(0,0.05,0.1,0.15,0.2) #asocial learning set to exp(asocialLP) = 1 asocialLP<-rep(0,noInd) asoc<-cbind(rep(0,noInd)) #Set base rate of asocial learning baseRate<-1/100 #set up records to record outputs of models sRecord<-aicRecord<-pRecord<-sRecordW<-aicRecordW<-pRecordW<-TsRecord<-TaicRecord<-TpRecord<-TsRecordW<-TaicRecordW<-TpRecordW<-profLikForS <-profLikForSW <-TprofLikForS <-TprofLikForSW <-estPropSolvesST<-truePropSolvesST<-TestPropSolvesST<-array(NA,dim=c(noSims,length(sVect),length(NoiseVect))) for(l in 1:length(NoiseVect)){ Noise<-NoiseVect[l] for(k in 1:length(sVect)){ s<-sVect[k] for(j in 1:noSims){ #generate a social network socialNet<-matrix(0, ncol=noInd, nrow=noInd) for(i in 0:9){socialNet[i*10+(1:10),i*10+(1:10)]<-runif(100,0.5,1)} #simulate diffusion z<-rep(0,noInd) orderAcq<-timeAcq <-rep(NA,noInd) runningTime<-0 for(i in 1:noInd){ rate<-baseRate*(exp(asocialLP)+s*z%*%t(socialNet))*(1-z) times<-rexp(noInd,rate) times[is.nan(times)]<-Inf orderAcq [i]<-which(times==min(times))[1] runningTime<-runningTime+min(times) timeAcq[i]<-runningTime z[which(times==min(times))[1]]<-1 } #set recorded social net according to level of noise required socialNet<-socialNet+rnorm(noInd*noInd,0,sd=Noise) socialNet[socialNet<0]<-0 socialNet[socialNet<0]<-1 #Fit OADA and TADA oaDataObject<-oaData(assMatrix= socialNet, taskid="1", groupid="1", asoc=asoc, orderAcq= orderAcq) taDataObject<-taData(timeAcq,max(timeAcq)+1,oadata= oaDataObject) model<-addFit(oaDataObject,interval=c(0,9999)) Tmodel<-tadaFit(taDataObject) #Record outputs from OADA #Record s sRecord[j,k,l]<-model@optimisation$minimum #Record AIC aicRecord[j,k,l]<-model@aicc #Record p value pRecord[j,k,l]<-model@ LRTsocTransPV #Record the profile likelihood for the real value of s (used to determine if it is within CI) profLikForS[j,k,l]<-profileLikelihoodAddOADA(s,1, model, oaDataObject)-model@optimisation$objective #Estimated proportion of solves by social transmission estPropSolvesST[j,k,l]<-oadaPropSolveByST(par=model@optimisation$minimum,oadata= oaDataObject,innovations=1,additive=T,asocialVar=NULL) #True proportion of solves by social transmission truePropSolvesST[j,k,l]<-oadaPropSolveByST(par=s,oadata= oaDataObject,innovations=1,additive=T,asocialVar=NULL) #Record the equivalent things for TADA TsRecord[j,k,l]<-Tmodel@optimisation$par[1] TaicRecord[j,k,l]<-Tmodel@aicc TpRecord[j,k,l]<-Tmodel@ LRTsocTransPV TprofLikForS[j,k,l]<-profileLikelihoodTADA(s,1, Tmodel, taDataObject)-Tmodel@optimisation$objective TestPropSolvesST[j,k,l]<-oadaPropSolveByST(par=Tmodel@optimisation$par[1],oadata= oaDataObject,innovations=1,additive=T,asocialVar=NULL) #Save it all in an object save(sRecord,aicRecord,pRecord,sRecordW,aicRecordW,pRecordW,TsRecord,TaicRecord,TpRecord,TsRecordW, TaicRecordW,TpRecordW,profLikForS ,profLikForSW ,TprofLikForS , file="NBDANoiseinNetRandomss100.dat") } } }