Functions through Poisson distribution in R?

Question

I am trying to run a trial that runs a couple of changing variables through a poisson distribution and gives an output. It goes through three years and there are 12 clusters. All of the groups 1-6 are control and 7 to 12 are intervention.

Year one would have clusters 1-12. with 6 control and 6 intervention. same for year two and three. So in year one you get a v and that v is the same for each cluster but each cluster has a different lambda.

This is the basic code and equation.

mu=3.4
lambda=rnorm(1,0,.27)
v=rnorm(1,0,.53)
e=rnorm(1,0,.74)
x=rpois(1, exp(mu+lambda+v+e))

Lambda and v are normally distributed with those st. deviations for each group. I would like in say year 1 for it to take v=rnorm(1,0,.53) and run that through the equation. At the same time I want it to take lambda=rnorm(1,0,.27) and run that for each of the 12 clusters with the same v from year one. So basically one v and 12 lambdas per year. I believe this can be done with functions but I am having a difficult time with it. I tried using for loops but I don't think it is doing exactly what I want. Here is the code with the for loops:

mu=3.4
lambda=numeric(length=12)
v=numeric(length=3)
for (i in 1:120{
lambda[i]=rnorm(1,0,.27)
}
for (j in 1:3){
v[j]=rnorm(1,0,.53)
}
e=rnorm(1,0,.74)
x=rpois(1, exp(mu+lambda+v+e))

I am not really sure if this is doing that correct thing or not so if someone with more experience could explain what my code is doing and how to make it do what I want with functions for lambda and v I would greatly appreciate it. Thanks for the help.

Answer 1

It sounds like you have 12 groups and 3 years for a total of 36 data points. For each year i you have a single value v[i] and 12 group-specific values lambda[1], lambda[2], ..., lambda[12]. You didn't really specify how each e is drawn, so let's assume there's a global e value.

You can build up a data frame with the relevant parameters for every group and year:

dat <- expand.grid(group=1:12, year=1:3)

Then you can specify each parameter in your data frame:

set.seed(144)
dat$lambda <- rnorm(36, 0, .27)
dat$v <- ave(dat$year, dat$year, FUN=function(x) rnorm(1, 0, .53))
dat$e <- rnorm(1, 0, .74)
dat$mu <- 3.4

Finally you can compute the x values from the poisson distribution based on all the other specified parameters.

dat$x <- apply(dat, 1, function(x) rpois(1, exp(sum(x[3:5]))))
dat
#    group year       lambda         e           v x
# 1      1    1 -0.445650165 0.4905398 -0.68348553 0
# 2      2    1  0.162758849 0.4905398 -0.68348553 0
# 3      3    1 -0.127948650 0.4905398 -0.68348553 0
# 4      4    1 -0.485355518 0.4905398 -0.68348553 0
# 5      5    1 -0.383702670 0.4905398 -0.68348553 1
# 6      6    1  0.042899625 0.4905398 -0.68348553 1
# 7      7    1  0.035036875 0.4905398 -0.68348553 3
# 8      8    1 -0.339165570 0.4905398 -0.68348553 0
# 9      9    1  0.039870111 0.4905398 -0.68348553 0
# 10    10    1  0.264376944 0.4905398 -0.68348553 2
# 11    11    1 -0.158722934 0.4905398 -0.68348553 2
# 12    12    1  0.065626713 0.4905398 -0.68348553 0
# 13     1    2 -0.119740653 0.4905398 -0.94943812 1
# 14     2    2 -0.273056956 0.4905398 -0.94943812 0
# 15     3    2  0.082479733 0.4905398 -0.94943812 1
# 16     4    2  0.004055817 0.4905398 -0.94943812 3
# 17     5    2 -0.136651972 0.4905398 -0.94943812 0
# 18     6    2 -0.379565000 0.4905398 -0.94943812 1
# 19     7    2 -0.347022844 0.4905398 -0.94943812 0
# 20     8    2 -0.025990597 0.4905398 -0.94943812 0
# 21     9    2  0.136014382 0.4905398 -0.94943812 3
# 22    10    2  0.063815822 0.4905398 -0.94943812 0
# 23    11    2  0.067932957 0.4905398 -0.94943812 1
# 24    12    2  0.099834989 0.4905398 -0.94943812 0
# 25     1    3 -0.135519422 0.4905398 -0.04616012 3
# 26     2    3 -0.226890654 0.4905398 -0.04616012 2
# 27     3    3 -0.024936233 0.4905398 -0.04616012 1
# 28     4    3 -0.182543225 0.4905398 -0.04616012 2
# 29     5    3  0.537843637 0.4905398 -0.04616012 2
# 30     6    3  0.291349032 0.4905398 -0.04616012 3
# 31     7    3  0.255455171 0.4905398 -0.04616012 0
# 32     8    3 -0.053874069 0.4905398 -0.04616012 5
# 33     9    3  0.052474650 0.4905398 -0.04616012 0
# 34    10    3 -0.317645288 0.4905398 -0.04616012 0
# 35    11    3  0.162754030 0.4905398 -0.04616012 1
# 36    12    3  0.337018922 0.4905398 -0.04616012 8