¿Cómo simular la eficacia del tratamiento en R?

Question

Digamos que quiero escribir una simulación para la tabla de abajo para decidir si el tratamiento con xilitol y las infecciones de oído son independientes. ¿Cómo lo haría?

Answer 1

(No queda claro si no quiere ninguno de los dos, sólo uno de ellos, o ambos totales marginales fijos. @Glen_b ha proporcionado un algoritmo de simulación que se basa en tener ambos totales marginales fijos; a continuación, proporciono algoritmos para las tres posibilidades).

La suposición de independencia significa que las probabilidades de las celdas son iguales al producto de la probabilidad de que una observación ocurra dentro de una fila dada por la probabilidad de que la observación ocurra dentro de la columna dada.

Ningún margen es fijo:

Utilizando los valores de su tabla de contingencia, el siguiente código simulará el nulo. Tenga en cuenta, sin embargo, que el número exacto de, por ejemplo, "yes" observaciones en cada iteración no serán necesariamente iguales 117 . No obstante, la probabilidad de que una observación esté en el Placebo Gum fila Y estar en la "yes" es igual al producto de la probabilidad de la fila por la probabilidad de la columna, que es la definición de independencia. (Obsérvese además que para obtener una tabla simple, de una sola simulación, basta con establecer 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)) })

A partir de aquí, si sólo has hecho una única tabla (por haber puesto B = 1 ), y sólo quería verlo, podría utilizarlo:

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

Para realizar una simulación completa del nulo, asegúrese de que B era algún número grande y el uso:

 # 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

Sólo se fijan los totales de las filas:

En este caso, siempre tendrás, digamos, exactamente 179 observaciones en el Placebo Gum tratamiento. No obstante, la probabilidad de estar en el "yes" columna será siempre la misma:

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

Ambos márgenes son fijos:

En este caso, siempre tendrás, digamos, exactamente 179 observaciones en el Placebo Gum tratamiento, y, digamos, exactamente 117 observaciones en el "yes" columna. No obstante, la probabilidad de estar en la Placebo Gum fila Y estar en la "yes" es igual al producto de la probabilidad de la fila por la probabilidad de la columna:

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

Answer 2

Si quieres simular una celda aleatoria (bajo independencia) con márgenes fijos, eso es efectivamente un muestreo hipergeométrico, que podemos aplicar recursivamente, por lo que una aproximación es

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

En el caso de $3\times 2$ tablas como la suya (o $k\times 2$ de forma más general), sólo hay que simular dos ( $k-1$ ), y el resto están determinados. Si se observa la celda (1,1) se puede tratar la situación como $2\times 2$ (combinando las categorías de las filas restantes) y así generar la celda (1,1); entonces se determina (1,2). Después de restar la primera fila de los totales de las columnas, queda un $2\times 2$ (más generalmente $(k-1)\times 2$ ) para las filas inferiores, lo que se hace de la misma manera.

[ Nota gung sugiere un enfoque más sencillo de entender y (en algunos casos), quizás más rápido, para la simulación con márgenes fijos en los comentarios, y da algo de código en su respuesta].

En R, puede utilizar simplemente r2dtable utiliza el algoritmo de Patefield[1].

[1]: Patefield, W. M. (1981),
"Algoritmo AS159. Un método eficiente para generar tablas r x c con totales de fila y columna dados,"
Estadística aplicada 30 , 91-97.