How to simulate effectiveness of treatment in R?

Question

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?

Answer 1

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.

Answer 2

(It isn't clear if you want neither, only a particular one, or both marginal totals fixed. @Glen_b has provided a simulation algorithm that is based on having both marginal totals fixed; below, I provide algorithms for all three possibilities.)

The assumption of independence means that the cell probabilities are equal to the product of the probability an observation will occur within a given row times the probability that the observation will occur within the given column.

Neither margin is fixed:

Using the values from your contingency table, the following code will simulate the null. Note, however, that the exact number of, say, "yes" observations in each iteration will not necessarily equal 117. Nonetheless, the probability that an observation will be in the Placebo Gum row AND being in the "yes" column is equal to the product of the row probability times the column probability, which is the definition of independence. (Note further that to get a simple, single simulated table, just set B = 1.)

N     = 533               # this is the total number of observations in your table
r.ns  = c(178, 179, 176)  # these are the row counts
c.ns  = c(416, 117)       # these are the column counts
r.ps  = r.ns/N            # these are the row probabilities
c.ps  = c.ns/N            # these are the column probabilities
probs = r.ps%o%c.ps       # these are the cell probabilities under independence
probs
#           [,1]       [,2]
# [1,] 0.2606507 0.07330801
# [2,] 0.2621150 0.07371986
# [3,] 0.2577221 0.07248433
probs.v = as.vector(probs)      # notice that the probabilities read column-wise
probs.v
# [1] 0.26065071 0.26211504 0.25772205 0.07330801 0.07371986 0.07248433
cuts = c(0, cumsum(probs.v))    # notice that I add a 0 on the front
cuts
# [1] 0.0000000 0.2606507 0.5227658 0.7804878 0.8537958 0.9275157 1.0000000

set.seed(4847)                  # this makes the example exactly reproducible
B      = 10000                  # number of iterations in simulation
vals   = runif(N*B)             # generate random values / probabilities
 # cut the random uniform values into cell categories: 
cats   = cut(vals, breaks=cuts, labels=c("11","21","31","12","22","32"))
 # this reforms the single N*B vector into a matrix of N obs by B iterations:
cats   = matrix(cats, nrow=N, ncol=B, byrow=F)
 # here we get the number of observations w/i each cell for each iteration:
counts = apply(cats, 2, function(x){ as.vector(table(x)) })

From here, if you only made a single table (by having set B = 1), and just wanted to see it, you could use:

matrix(counts, nrow=3, ncol=2, byrow=T)  # if B had been 1
#      [,1] [,2]
# [1,]  137   36
# [2,]  125   47
# [3,]  148   40

To perform a full simulation of the null, make sure B was some large number and use:

 # some clean up of the workspace: 
rm(N, r.ns, c.ns, r.ps, c.ps, vals, probs, probs.v, cuts, cats) 
p.vals = vector(length=B)       # this will store the outputs
for(i in 1:B){
   # put the counts into the form that chisq.test() needs:
  mat       = matrix(counts[,i], nrow=3, ncol=2, byrow=T) 
  p.vals[i] = chisq.test(mat)$p.value  # here we store the p values
}
mean(p.vals<.05)  # we have 5% type I errors, as appropriate:
# [1] 0.0475

Only the row totals are fixed:

In this case, you will always have, say, exactly 179 observations in the Placebo Gum treatment. Nonetheless, the probability of being in the "yes" column will always be the same:

prob.yes = 117/533  # this is the probability of 'yes' under all 3 treatments

set.seed(192)       # this makes the example exactly reproducible
PG = rbinom(n=178, size=1, prob=prob.yes)  # these each generate a vector of 'yes'es &
XG = rbinom(n=179, size=1, prob=prob.yes)  #  'no's with the fixed row totals & the
XL = rbinom(n=176, size=1, prob=prob.yes)  #  constant probability of success

raw.observations = rbind(cbind("PG", PG),  # here I just make the dataset
                         cbind("XG", XG),
                         cbind("XL", XL) )
table(raw.observations[,1], raw.observations[,2])
#      0   1
# PG 142  36
# XG 133  46
# XL 140  36

Both margins are fixed:

In this case, you will always have, say, exactly 179 observations in the Placebo Gum treatment, and, say, exactly 117 observations in the "yes" column. Nonetheless, the probability of being in the Placebo Gum row AND being in the "yes" column is equal to the product of the row probability times the column probability:

X = rbind(cbind(rep("PG",129), rep("no", 129)),   # this just re-creates your table
          cbind(rep("XG",150), rep("no", 150)),
          cbind(rep("XL",137), rep("no", 137)),
          cbind(rep("PG", 49), rep("yes", 49)),
          cbind(rep("XG", 29), rep("yes", 29)),
          cbind(rep("XL", 39), rep("yes", 39)) )
table(X[,1],X[,2])
#     no yes
# PG 129  49
# XG 150  29
# XL 137  39

set.seed(6875)                 # this makes the simulation exactly reproducible
 # the sample() call is the key element:
X.perm = cbind(X[,1], sample(X[,2], nrow(X), replace=F))
table(X.perm[,1], X.perm[,2])
#     no yes
# PG 140  38
# XG 140  39
# XL 136  40