2
$\begingroup$

I am trying to randomly simulate a population of $N$ individuals among which a predefined number $K$ of them have an outcome. The trick is that I want to assign different probabilities to the individuals such that some are more 'at risk' than the others.

To give an example, let's say that I have $N = 10$ individuals and I know that $K = 4$ of them have the outcome. Besides, let's say that individuals 6 to 10 are twice as likely than the others to have the outcome, which means that a vector of probability would be $(1, 1, 1, 1, 1, 2, 2, 2, 2, 2) / 15$. How do I simulate such populations?

I tried using the R function sample without replacement to draw the indices of individuals having the outcome. However, for some reason, the obtained proportions in the end are not quite right.

Another way to put the problem is that I would like a distribution such as the multinomial one, but for which the count can't be larger than one for each category.

Is there such a distribution?

$\endgroup$
  • 1
    $\begingroup$ I do not understand the model: either each of the $N$ individuals has a probability $p_i$ to have the outcome, in which case the $p_i$'s do not have to sum up to one, or one and only one of the $N$ individuals is endowed with the outcome. $\endgroup$ – Xi'an Aug 13 at 12:45
  • 1
    $\begingroup$ This description doesn't determine a unique model. One common and natural way to model risk is to suppose that it depends on how much time is elapsed: on average, someone with twice the risk will experience an adverse event in half the time. But with this model, after the first $K$ adverse events have been observed in the population, you will tend to have slightly more than the expected numbers of low-risk individuals and slightly fewer than the expected numbers of high-risk individuals. $\endgroup$ – whuber Aug 13 at 15:57
  • 1
    $\begingroup$ @whuber Indeed that is what happened when I tried to sample the indices with varying probabilities. However, I have trouble seeing how it works. $\endgroup$ – Pierre Aug 13 at 17:22
1
$\begingroup$

I think you might want something like this:

set.seed(16)
# we have 4 events, and 10 people of which 5 of them should be twice as likely
# that's the same as having 15 people of which 5 belong to one group and 10 to the other
sample.fun <- function(){sample(1:15, size=4, replace=F)}
samples <- replicate(1e4, sample.fun())

# check whether that worked out ...
sum(colSums(samples < 6)) / (sum(colSums(samples < 6))+sum(colSums(samples >= 6)))
# the proportion of events in individuals 1-5 is 0.3335

We can further check whether the events per individual do make sense

samples <- ifelse(samples==11,6,
                      ifelse(samples==12,7, 
                             ifelse(samples==13,8,
                                    ifelse(samples==14,9, 
                                           ifelse(samples==15,10,samples)))))
table(samples) 
#    1    2    3    4    5    6    7    8    9   10 
# 2674 2653 2704 2590 2719 5267 5370 5388 5330 5305 

I bet the recode part could be done much better than my inelegant ifelse()-structure.

| cite | improve this answer | |
$\endgroup$
1
$\begingroup$

If each individual from the population has a probability $p_i$ to have a certain characteristic, draw the outcome for each of the $N$ individuals and accept the result if exactly $K$ enjoy the outcome. Something like

while(sum(o=(runif(N)<p)!=K){}
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ I think this is called the Fisher noncentral hypergeometric approach, and I worry that it may be inefficient when $N$ is large. More seriously, I think it may be biased, with this example giving probabilities of being sampled of about $\frac{9}{30}$ each for the less likely elements and about $\frac{15}{30}$ for the more likely elements, rather than the $\frac{8}{30}$ and $\frac{16}{30}$ that might be expected $\endgroup$ – Henry Aug 13 at 16:45
0
$\begingroup$

You can do it by computing a randomly variating probability for each individual to have the condition.

n1 = 10 # Number in Group 1
n2 = 10 # Number in Group 2
pop = c(rep(1,n1),rep(2,n2)) # Create a population of 1 and 2
prob = rnorm(length(pop),pop*.33,.2) # For each individual, compute a probability to have the condition
cond = rep(0,length(pop)) # Set conditions to 0
cond[prob>.5]=1 # If prob>0.5: The condition occurs
par(mfrow=c(1,2))
plot(as.factor(pop),prob)
plot(table(pop,cond))

This example gave the output below from which you can see that group 2 has a higher probability of getting the condition and also has visibly more cases. You can change the function for the random variation and the threshold for getting the condition to adjust the simulation to your needs.

Output of the two plotting functions

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.