data(finch)
## Simulate graphs and record the networks statistics ##
## specified by 'coincidences(x)' ##
sim<-simulate_sis(finch~coincidences(0:17), nsim=10000)
observed.stats<-summary(finch~coincidences(0:17))
sampled.stats<-sim %n% "samplestatistics"
## Calculate 95\% confidence intervals for the network statistics ##
library(Hmisc)
p<-exp(sim %n% "log.prob")
p<-p/sum(p)
maxs<-apply(sampled.stats,2,wtd.quantile,
weights=p,probs=0.975,normwt=TRUE)
mins<-apply(sampled.stats,2,wtd.quantile,
weights=p,probs=0.025,normwt=TRUE)
means<-apply(sampled.stats,2,wtd.mean,weights=p)
## Plot the means and CIs for the null distributions of the ##
## statistics. Also plot the observed statistics ##
plot(0:17, means, type="b", ylim=c(0,24.5), lwd=3, lty=3,
xlab="Number of Islands", ylab="Pairs of Finches")
for(i in 1:18)
{
points(rep(i-1,2), c(maxs[i],mins[i]), type="l", lwd=2)
}
points(0:17, observed.stats, type="b", pch=4, lwd=3)
## Calculate the p-value for 'coincidences(0)' ##
r0<-(p%*%sweep(sampled.stats,2,observed.stats,"<"))[1,]
r1<-(p%*%sweep(sampled.stats,2,observed.stats,">"))[1,]
round(apply(cbind(r0,r1),1,min),digits=8)[1]
Run the code above in your browser using DataLab