I am wondering how I can simulate Poisson data that is correlated. Let's say, I collect data at two time points. At both time points, the data is Poisson distributed, at time point 1 with $\lambda=4$, at time point 2 with $\lambda=3$ and there should be a strong correlation between measurement at time point 1 and 2, i.e. $\rho = 0.8$. Obviously, for time point 1, I can sample data of size $n$ in R with rpois(n=n,lambda=4). But how can I sample the data for time point 2 that is correlated to time point 1?

Some more details on the specific lab problem: The variable of interest is the number of successful 'catches' of a mouse (it's some kind of coordinative task). We count successes per mouse at time point 1. Then some treatment is applied to the mouse, that is supposed to affect its coordinative ability. We expect a drop in the average number of successes by around 25% in each mouse. Mice that were relatively good at the task before the treatment will probably still be relatively good after the treatment. A statistical comparison between time point 1 and 2 will be performed using the Wilcoxon signed-rank test. We want to know the sample size required to achieve power of 80% (while type I error rate should be 5 %).

Update: It turns out, the number of successes is limited to discrete values between 0 and 7. Further, based on previous experiments we have the following probabilities for each number of successes before treatment: $p=\{0.02,0.03,0.09,0.17,0.19,0.30,0.15,0.05\}$.

Now, I am considering to sample values between 0 and 8 with the given probabilities and use it as the "true ability" of an individual mouse, similarly to a random effects model. Then, I would sample the actual successes for each mouse from a binomial distribution, with 8 trials and probability $p$ corresponding to the "true ability". For the second time point, I would do the same, just use $p_2=p*0.75$ to represent the 25% drop in ability.

Does that make sense to anyone? Any recommendations?

  • 1
    $\begingroup$ We have many threads about this topic: stats.stackexchange.com/…. A glance over them indicates you need to be more specific about the form of dependency between the two variables, for otherwise there are many possible answers. Could you be more specific about what you are trying to simulate and why? $\endgroup$
    – whuber
    Commented Jan 29, 2020 at 16:29
  • 1
    $\begingroup$ @whuber Thanks a lot for your comment! I added some more details on the specific problem to solve. Indeed, there are many threads related to this topic of which many of them were already helpful for me. However, I think this specific problem was not discussed yet. $\endgroup$
    – LuckyPal
    Commented Jan 29, 2020 at 16:40
  • 1
    $\begingroup$ But that doesn't matter, because you are comparing changes in performance. $\endgroup$
    – whuber
    Commented Jan 29, 2020 at 17:17
  • 1
    $\begingroup$ That is a valid point. But what distribution do the changes in performance follow, if the original performances are Poisson distributed? Is it the Skellam distribution, as indicated here math.stackexchange.com/questions/2934636/…? $\endgroup$
    – LuckyPal
    Commented Jan 29, 2020 at 17:21
  • 1
    $\begingroup$ The specific distribution does not matter, actually. We were considering a Poisson distribution simply for the matter of effect size anticipation, since a reduction in means is easier to understand than the rather abstract effect that the Wilcoxon test is actually testing. But any formulation, where we can formulate an effect size in terms of change in means would work. $\endgroup$
    – LuckyPal
    Commented Jan 29, 2020 at 18:45

1 Answer 1


In the end, I came up with the following approach, which is basically a hierarchical model.

The historical data was used as "true" distribution of successful catches in 7 tasks. By sampling from this distribution, I create a set of the true abilities of n mice. Assuming that each mouse has a random component in its coordinative task and won't always perform according to its true ability, I used each mouse's true ability as proportion p in a binomial distribution to sample successes in k=7 trials. Due to the treatment, each mouse's individual p is then lowered by 25% (of course, I could have assumed a random component on the treatment effect, too, but I thought that could be a slight overkill). Finally, a Wilcoxon signed-rank test is used to compare the number of successes before and after treatment.

library(coin) # since this package can handle ties in the data
y <- c(rep(0,2),rep(1,3),rep(2,9),rep(3,17),rep(4,19),rep(5,30),rep(6,15),rep(7,5)) # historical data


  for(i in 1:nsim){
    smp <- sample(y, n, replace=TRUE)
    before.treat <- rbinom(n=n,size=7,prob=smp/7)
    after.treat <- rbinom(n=n,size=7,prob=(smp/7)*0.7)
      p.value[i] <- 1 # if there are no differences at all, we want a non-significant p-value
    } else { 
      test <- wilcoxsign_test(before.treat ~ after.treat)
      p.value[i] <- pvalue(test)
  print(paste("N=",n,", power=",power))

> [1] "N= 15 , power= 0.7182"
> [1] "N= 16 , power= 0.7413"
> [1] "N= 17 , power= 0.7829"
> [1] "N= 18 , power= 0.8077"

I would be happy to take any commments on this approach. Particularly, I noted that the Monte Carlo simulations seem very sensitive to the random seed used and larger simulations do take a lot of time, so I am a bit worried about the validity of the results.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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