Sample size calculation for correlated count data 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?
 A: 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

nsim=1e4
n=14
power=0
p.value<-c()

while(power<0.8){
  n=n+1
  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)
    if(all((before.treat-after.treat)==0)){
      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)
    }
  }
  power<-sum(p.value<0.05)/nsim
  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.
