How to simulate effectiveness of treatment in R? Let's say I want to write a simulation for the table below to decide if Xylitol treatment and ear infections are independent. How would I go about doing this?

 A: If you want to simulate one random cell (under independence) with fixed margins, that's effectively hypergeometric sampling, which we can apply recursively, so one approach is
pick one cell; 
  given the margins that cell has a hypergeometric distribution, so
  simulate from that hypegeometric 
once you have that value, that affects possible values of other cells, which
can be generated in turn, each conditional on all previous values

In the case of $3\times 2$ tables such as yours (or $k\times 2$ tables more generally), you need only simulate two ($k-1$) values, and the rest are determined. If you look at the (1,1) cell you can treat the situation as $2\times 2$ (by combining the remaining row categories) and so generate the (1,1) cell; then (1,2) is determined. After subtraction of the first row from the column totals you're then left with a $2\times 2$ (more generally $(k-1)\times 2$) table for the lower rows which is then done in the same fashion.
[Note: gung suggests a simpler-to-understand and (in some cases), perhaps faster approach to simulation with fixed margins in the comments, and gives some code in his answer.]
In R, you can just use r2dtable; it uses Patefield's algorithm[1].
[1]: Patefield, W. M. (1981),
"Algorithm AS159. An efficient method of generating r x c tables with given row and column totals,"
Applied Statistics 30, 91–97. 
