code for rebuttal of reviewer comments, additional Bayes analysis

gistfile1.txt

n1 = 9196 # total number of lysis positive genomes
y1 = 830  # number of lyssis positive genemes among colicin B positive genomes
n2 = 303 # total number of lysis negative genomes
y2 = 20 # lysis negative genomes among colicin B positive genomes

# SIMULATION
I = 10000 # simulations
theta1 = rbeta(I, y1 + 1, (n1 - y1) + 1)  
theta2 = rbeta(I, y2 + 1, (n2 - y2) + 1)
diff = theta1 - theta2  # simulated diffs

# OUTPUT
quantiles = quantile(diff, c(0.005,0.025,0.5,0.975,0.995))
print(quantiles,digits = 2)

print("Probability higher % lysis positive than lysis negative on colicinB:")
print(mean(theta1 > theta2))