#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<-c(0,0.1,0.2,0.3,0.4,0.5,1,1.5,2,3,4) #Set variation in Bj to be considered BNoiseVect<-c(0,0.1,0.25,0.5,1) #asocial learning set to exp(asocialLP) = 1 asocialLP<-rep(0,noInd) asoc<-cbind(rep(0,noInd)) 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 <-TestPropSolvesST<- truePropSolvesST <-estPropSolvesSTW <-TestPropSolvesSTW <-array(NA,dim=c(noSims,length(sVect),length(BNoiseVect))) for(l in 1:length(BNoiseVect)){ BNoise<-BNoiseVect[l] for(k in 1:length(sVect)){ s<-sVect[k] for(j in 1:noSims){ #generate a social network socialNet<-matrix(runif(noInd+noInd,0,0.3), ncol=noInd, nrow=noInd) for(i in 0:9){socialNet[i*10+(1:10),i*10+(1:10)]<-runif(100,0.5,1)} #Add in variation in performance rate BVect<-exp(rnorm(noInd,log(1),sd=BNoise)) z<-rep(0,noInd) orderAcq<-timeAcq <-rep(NA,noInd) runningTime<-0 #simulate diffusion for(i in 1:noInd){ rate<-baseRate*(exp(asocialLP)+s*z%*%(t(socialNet)*BVect))*(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 } #Fit OADA and TADA models weighted and unweighted oaDataObject<-oaData(assMatrix= socialNet, taskid="1", groupid="1", asoc=asoc, orderAcq= orderAcq) oaDataObjectW<-oaData(assMatrix= socialNet, taskid="1", groupid="1", asoc=asoc, orderAcq= orderAcq,weights= BVect) taDataObject<-taData(timeAcq,max(timeAcq)+1,oadata= oaDataObject) taDataObjectW<-taData(timeAcq,max(timeAcq)+1,oadata= oaDataObjectW) model<-addFit(oaDataObject,interval=c(0,9999)) modelW<-addFit(oaDataObjectW,interval=c(0,9999)) Tmodel<-tadaFit(taDataObject) TmodelW<-tadaFit(taDataObjectW) #Record outputs from OADA and TADA (see "Simulations association net noise.R" for explanation) #W suffix indicates output from a weighted model sRecord[j,k,l]<-model@optimisation$minimum aicRecord[j,k,l]<-model@aicc pRecord[j,k,l]<-model@ LRTsocTransPV profLikForS[j,k,l]<-profileLikelihoodAddOADA(s,1, model, oaDataObject)-model@optimisation$objective sRecordW[j,k,l]<-modelW@optimisation$minimum aicRecordW[j,k,l]<-modelW@aicc pRecordW[j,k,l]<-modelW@ LRTsocTransPV profLikForSW[j,k,l]<-profileLikelihoodAddOADA(s,1, modelW, oaDataObjectW)-modelW@optimisation$objective 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 TsRecordW[j,k,l]<-TmodelW@optimisation$par[1] TaicRecordW[j,k,l]<-TmodelW@aicc TpRecordW[j,k,l]<-TmodelW@ LRTsocTransPV TprofLikForSW[j,k,l]<-profileLikelihoodTADA(s,1, TmodelW, taDataObjectW)-TmodelW@optimisation$objective estPropSolvesST[j,k,l]<-oadaPropSolveByST(par=model@optimisation$minimum,oadata= oaDataObject,innovations=1,additive=T,asocialVar=NULL) estPropSolvesSTW[j,k,l]<-oadaPropSolveByST(par=modelW@optimisation$minimum,oadata= oaDataObjectW,innovations=1,additive=T,asocialVar=NULL) truePropSolvesST[j,k,l]<-oadaPropSolveByST(par=s,oadata= oaDataObject,innovations=1,additive=T,asocialVar=NULL) TestPropSolvesST[j,k,l]<-oadaPropSolveByST(par=Tmodel@optimisation$par[1],oadata= oaDataObject,innovations=1,additive=T,asocialVar=NULL) TestPropSolvesSTW[j,k,l]<-oadaPropSolveByST(par=TmodelW@optimisation$par[1],oadata= oaDataObjectW,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 ,TprofLikForSW, estPropSolvesST , TestPropSolvesST , truePropSolvesST , estPropSolvesSTW , TestPropSolvesSTW ,file="NBDANoiseinBAverageRate1.dat")