¿Cómo crear un conjunto de datos con probabilidad condicional?

Supongamos que una determinada enfermedad ($D$) tiene una prevalencia de $\dfrac3{1000}$. Supongamos también que cierto síntoma ($S$) tiene una prevalencia (en la población general = personas con esa enfermedadD y personas sin esa enfermedad [probablemente con otra enfermedad, pero no es importante]) de $\dfrac5{1000}$. En una investigación previa, se descubrió que la probabilidad condicional$P(S|D) = 30\%$ (la probabilidad de tener el síntoma $S$, dada la enfermedad $D$ es $30\%$)

Primera pregunta : podría ser$P(S|D)$ interpretado como equivalente a la prevalencia del síntoma $S$ en el grupo de personas que tienen la enfermedad $D$?

Segunda pregunta : quiero crear en R un conjunto de datos, que muestra que:

¿Como hacer esto? Si uso simplemente la samplefunción, mi conjunto de datos carece de la información que$P(S|D)=30\%$:

symptom <- sample(c("yes","no"), 1000, prob=c(0.005, 0.995), rep=T)
disease <- sample(c("yes","no"), 1000, prob=c(0.002, 0.998), rep=T)

Entonces mi pregunta es: ¿cómo crear un buen conjunto de datos, incluida la probabilidad condicional que deseo?

EDITAR : publiqué la misma pregunta también en stackoverflow.com ( /programming/7291935/how-to-create-a-dataset-with-conditional-probability ), porque, en mi opinión, mi pregunta se hereda del programa de lenguaje R, pero también de la teoría estadística.

Conoces las siguientes probabilidades marginales

                Symptom        Total
                Yes     No
Disease Yes      a       b     0.003
        No       c       d     0.997  
Total           0.005   0.995  1.000

y eso se a/(a+b) = 0.3convierte en

                Symptom        Total
                Yes     No
Disease Yes     0.0009  0.0021 0.003
        No      0.0041  0.9929 0.997  
Total           0.005   0.995  1.000

y de hecho a/(a+c) = 0.18como dijiste.

Entonces en R podrías codificar algo como

diseaserate <- 3/1000
symptomrate <- 5/1000
symptomgivendisease <- 0.3

status  <- sample(c("SYDY", "SNDY", "SYDN", "SNDN"), 1000, 
            prob=c(diseaserate * symptomgivendisease,
                   diseaserate * (1-symptomgivendisease),
                   symptomrate - diseaserate * symptomgivendisease,
                   1 - symptomrate - diseaserate * (1-symptomgivendisease)),
            rep=TRUE)
symptom <- status %in% c("SYDY","SYDN")
disease <- status %in% c("SYDY","SNDY")

aunque debe tener en cuenta que 1000 es una muestra pequeña cuando uno de los eventos tiene una probabilidad de que ocurra 0.0009.

La tablefunción devuelve un objeto tipo matriz:

> symptom <- sample(c("yes","no"), 100, prob=c(0.2, 0.8), rep=TRUE)
> disease <- sample(c("yes","no"), 100, prob=c(0.2, 0.8), rep=TRUE)
> dataset <- data.frame(symptom, disease)
> dst_S_D <-with(dataset, table(symptom, disease))
> dst_S_D
       disease
symptom no yes
    no  65  13
    yes 17   5

Entonces el Pr (D | S = "sí") =

> probD_Sy <- dst_S_D[2, 2]/sum(dst_S_D[2, ] )
> probD_Sy
[1] 0.2272727

Cambié el problema porque la primera vez que lo ejecuté con sus parámetros, obtuve:

> dst_S_D <-with(dataset, table(symptom, disease)); dst_S_D
       disease
symptom   no  yes
    no  9954   22
    yes   24    0

Y pensé que un Pr (D | S = "yes") de 0 era bastante aburrido. Si va a ejecutar esto muchas veces, debe construir una función y usar esa función con la replicatefunción.

Aquí hay un método para construir un conjunto de datos que aplique una probabilidad diferente de enfermedad en el grupo sintomático que sea 3 veces mayor que la utilizada en el grupo asintomático:

symptom <- sample(c("yes","no"), 10000, prob=c(0.02, 0.98), rep=TRUE)
dataset <- data.frame(symptom, disease=NA)
dataset$disease[dataset$symptom == "yes"] <- 
       sample(c("yes","no"), sum(dataset$symptom == "yes"), prob=c(0.15, 1-0.15), rep=TRUE)
dataset$disease[dataset$symptom == "no"] <- 
        sample(c("yes","no"), sum(dataset$symptom == "no"), prob=c(0.05, 1-0.05), rep=TRUE)
 dst_S_D <-with(dataset, table(symptom, disease)); dst_S_D
#       disease
symptom   no  yes
    no  9284  509
    yes  176   31

Yo diría que su pregunta no es realmente tan dependiente del lenguaje R, y más apropiado aquí, porque, para ser franco, la generación de datos como esta es principalmente una tarea estadística, en lugar de una programación.

Primera pregunta: p (S | D) es el riesgo de tener un síntoma S en una población con enfermedad D. Puede ser directamente comparable a la prevalencia con ciertas advertencias, como el síntoma que no tiene impacto en la duración de la enfermedad. Considere el siguiente ejemplo: Uno de los síntomas de SuperEbola es Instant Death, con p (Death | Super Ebola) = 0.99. Aquí, su prevalencia del síntoma en realidad sería extremadamente baja (de hecho, 0.00) ya que nadie a quien pueda tomar muestras con la enfermedad tiene el síntoma.

Segunda pregunta: Volvería a esto de una manera un tanto gradual. Primero, calcule el riesgo inicial del síntoma que necesitará para obtener 0.15 en toda la población, teniendo en cuenta que el 0.03% de su población tendrá una tasa más alta. Entonces esencialmente genera dos probabilidades:

• Riesgo de enfermedad = 0.003
• Riesgo de síntomas = riesgo inicial calculado + aumento relativo debido a la enfermedad * indicador binario del estado de la enfermedad

Luego genera dos números aleatorios uniformes. Si el primero es inferior a 0,003, tienen la enfermedad. Eso luego se introduce en el cálculo del riesgo para el segundo, y si el número aleatorio de cada individuo es menor que su riesgo, tienen el síntoma.

Esta es una forma poco elegante y poco elegante de hacer las cosas, y es probable que alguien llegue con un enfoque mucho más eficiente. Pero encuentro en los estudios de simulación que detallan cada paso en el código, y es útil mantenerlo tan cerca de cómo vería un conjunto de datos en el mundo real.

Primera pregunta:

Sí, por supuesto, esa es casi la definición, aunque tendrá algún error asociado con el tamaño de su muestra. es decir, esto es exactamente correcto en un tamaño de muestra infinito.

Segunda pregunta:

Esto se llama Teorema de Bayes , pero supongo que ya lo sabes. Ahora, dada la información que ha proporcionado, obtengo la probabilidad de P (D | S) como 0.18 o 18%:

P(S|D)P(D)
----------
   P(S)

  0.3*(3/1000)
= ------------
    (5/1000)

= 0.18

Ahora, desafortunadamente, no estoy muy familiarizado con R, así que realmente no puedo ayudarte con un programa exacto. Pero seguramente las cantidades de personas que caen en cada grupo son bastante fáciles de calcular:

Para su conjunto de muestra 10000 necesita:

1. 50 personas con síntomas (población * P (S))
2. 9 personas deben tener síntomas y la enfermedad (50 * P (D | S))
3. 21 personas con la enfermedad y sin síntomas (población * P (D) = 30 y ya tenemos 9)

Lo que debería hacer que generar una población adecuada sea bastante trivial.