How to create a dataset with conditional probability?

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 diseaseD 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.

• Common courtesy is to denote that you've cross-posted at another SE site. stackoverflow.com/questions/7291935/… Commented Sep 3, 2011 at 8:07
• I flagged your question on SO for migration. Please, don't cross-post!
– chl
Commented Sep 4, 2011 at 11:46

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.

• Awesome solution, it works great! Now I can create a dataset showing what the Bayes formula can calculate. Thanks a lot! Commented Sep 3, 2011 at 15:54
• Told you someone would come along with something more elegant ;) Commented Sep 3, 2011 at 18:51
• @henry I would be really happy if you can give a look at my new question here: stats.stackexchange.com/questions/15202/…. It is a generalization of this question, with 2 symptoms. Commented Sep 5, 2011 at 7:32

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

• Perfect, nice and elegant trick! I added some new information in my answer, to formalize better what I'm looking for. Commented Sep 3, 2011 at 15:36

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.

• Thanks for the answer; the SuperEbola example is really educative and useful! The rest of your answer remains quite unclear, to me, especially when you say "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". How to compute this baseline risk? Commented Sep 3, 2011 at 9:33
• Honestly, its a pain to do. If I were you, I'd change my example slightly - rather than asserting that the overall risk in the population is 0.15, I'd say the baseline risk in the non-diseased is, say, 0.15 or 0.10, then determine the increase in risk I want in the diseased and let the overall risk fall where it may, rather than trying to set it. It's considerably easier to code, though you will possibly not have numbers that are quite as clean at the end. Commented Sep 3, 2011 at 18:47

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.

• Yes, the true value is 0.18, sorry for wrong typing. The second part of your answer is correct, but the problem is to create a dataset (in R) that really has 9 people with disease and symptom. The "sample" function correctly creates 50 and 30 "yes" for, respectively, symptom and disease; but it does not ensure that 9 people (out of 30) are also in the group of "yes-disease". Commented Sep 3, 2011 at 9:31
• Again afraid you might need someone more familiar with R than myself to help you on the use of this sample function. However, you could always generate a much larger population and then randomly pick 10000 samples from that.
– Charles Keepax
Commented Sep 3, 2011 at 12:10