How to create a dataset with conditional probability?

Question

Suppose that a certain disease ($D$) has a prevalence of $\dfrac3{1000}$. Also suppose that a certain symptom ($S$) has a prevalence (in the general population = people with that disease$ $D and people without that disease [probably with other disease, but it's not important]) of $\dfrac5{1000}$. In a previous research, it was discovered that the conditional probability $P(S|D) = 30\%$ (the probability to have the symptom $S$, given the disease $D$ is $30\%$).

First question: Could be $P(S|D)$ interpreted as equivalent to the prevalence of the symptom $S$ in the group of people having the disease $D$?

Second question: I want to create in R a dataset, which shows that:

$$P(D|S) = \frac{P(S|D)P(D)} {P(S)}$$ With my fictional data, we can compute $P(D|S)=0.18$, that is interpreted in this way: given a patient with the symptom $S$, the probability that he has the disease $D$ is $18\%$.

How to do this? If I use simply the sample function, my dataset is lacking of the information that $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)

So my question is: how to create a good dataset, including the conditional probability I desire?

EDIT: I posted the same question also on stackoverflow.com (https://stackoverflow.com/questions/7291935/how-to-create-a-dataset-with-conditional-probability), because, in my opinion, my question is inherited to the R language program, but also to statistical theory.

Answer 1

You know the following marginal probabilities

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

and that a/(a+b) = 0.3 so this becomes

                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

and indeed a/(a+c) = 0.18 as you stated.

So in R you could code something like

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")

though you should note that 1000 is a small sample when one of the events has a probability of 0.0009 of happening.

Answer 2

The table function returns a matrix-like object:

> 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

So the Pr(D|S="yes") =

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

I changed the problem because the first time I ran it with your parameters, I got:

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

And I thought a Pr(D|S="yes") of 0 was rather boring. If you are going to run this many times you should construct a function and use that function with the replicate function.

Here is a method of constructing a dataset that applies a different probability of disease in the symptomatic group that it 3 times higher than is used in the asymptomatic group:

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

Answer 3

I'd argue your question isn't really that heavily dependent on the R language, and more appropriate here, because - to be blunt - the generation of data like this is mostly a statistical task, rather than a programming one.

First Question: p(S|D) is the risk of having symptom S in a population with disease D. It can be directly comparable to the prevalence with certain caveats, like the symptom having no impact on disease duration. Consider the following example: One of the symptoms of SuperEbola is Instant Death, with p(Death | Super Ebola) = 0.99. Here, your prevalence of the symptom would actually be extremely low (indeed, 0.00) as no one you can sample with the disease has the symptom.

Second Question: I would back into this in a somewhat stepwise fashion. First, calculate the baseline risk of the symptom you'll need to get 0.15 in the whole population, taking into account that 0.03% of your population will be at higher rate. Then essentially generate two probabilities:

• Risk of disease = 0.003
• Risk of symptom = calculated baseline risk + relative increase due to disease * binary indicator of disease status

Then generate two uniform random numbers. If the first is less than 0.003, they've got the disease. That then gets fed into the risk calculation for the second, and if the random number for each individual is less than their risk, they've got the symptom.

This is sort of a plodding, inelegant way to do things, and its likely someone will come by with a far more efficient approach. But I find in simulation studies spelling each step out in the code, and keeping it as close to how I would see a data set in the real world is useful.

Answer 4

First question:

Yes of course that is almost the definition, although you will have some error associated with your sample size. i.e. This is only exactly correct at an infinite sample size.

Second question:

This is called Bayes Theorem, but I presume you already know that. Now given the information you have provided I get the probability of P(D|S) as 0.18 or 18%:

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

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

= 0.18

Now unfortunately, I am not too familiar with R so can't really help you out with an exact program. But surely the quantities of people that fall into each group are quite easy to calculate:

For your 10000 sample set you need:

1. 50 people with symptoms (population*P(S))
2. 9 people should have symptoms and the disease (50*P(D|S))
3. 21 people with the disease and no symptoms (population*P(D)=30 and we already have 9)

Which should make generating a suitable population fairly trivial.