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?

• What do you mean by "a simulation for the table"? Do you mean simulate under the assumption of independence? Is this treating both margins as fixed? Or in some way treating only one margin (presumably row totals) as fixed? (and in that case, how do you do want it to work?) Jan 21 '15 at 3:08

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.

• +1 for the fixed margins solution. In that case, is there any reason you prefer this over a permutation simulation / test? Jan 21 '15 at 3:44
• I looked up r2dtable but don't quite understand how to translate the table into this function. The documentation says that the rows and columns must both add up to the same value, which I'm kind of confused about. (Sorry I'm quite new to R)
– Gerk
Jan 21 '15 at 4:11
• @gung I was simply responding to a question about simulating tables. You could certainly use this simulation approach to perform an exact test using some (as yet unspecified) statistic, which I presume is the intent of the OP. But for simulation, I'm not sure what distinction you're making; could you clarify -- in respect of a simulation test, what are you suggesting? As far as a permutation test goes, again, what test are you suggesting? If you're saying "how is this distinct from Fisher's exact test?" -- it may not be (aside randomness of simulation), depending on what test statistic is used Jan 21 '15 at 4:14
• Let X be a matrix where each observation is a row, & where the observed column category is in col1 & the observed row category is in col2. Then you could run, say, Fisher's exact test as fisher.test(table(X[,1],X[,2])). From there, you can have independence & constant marginals by simply shuffling one of the columns of X, eg X.perm = cbind(X[,1], sample(X[,2], nrow(X), replace=F)). Jan 21 '15 at 4:19
• @gung Ah yes, thanks -- now I see what you're getting at. That simulation shouldn't differ in distribution from mine, and should be (up to simulation error) the same as Fishers test if you're using likelihood under the null as the test statistic. I think your suggested approach is easier to understand, and likely a good bit faster with small counts and large tables. Jan 21 '15 at 4:27

(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