[R] lme, two random effects, poisson distribution

[Martina Pavlicova, PhD]

Hello,

I have a dataset concerning slugs. For each slug, the number of
pumps per one time slot was counted. The number of pumps follows
Bi(30, p) where p is very small, thus could be approximated by
Poisson dist. (# of pumps is very often = 0)

The slugs were observed during 12 time slots which are correlated in
time as AR(1). The time slots are divided into two categories:
      Resting time slots (the first 10)
      Excited time slots (the last 2)

I used model:

************pumps_ti = state_t + slugs_i + error_ti

slugs and error are normaly distributed

      pumps_ti - # of pumps for i-th animal and t-th time slot
      x_t - order of the time slot (x_1 = 1, ..., x_12 = 12)
      state_t - state_t = 0 for resting time slots (t=1,...,10)
              state_t = 1 for excited time slots (t=11,12)
      slugs_i - ith animal, where i = 1,...,25

I would like to find out if the # of pumps depends on the variable
state, assuming the correlation AR(1) between x_t and slugs being a
random-effect on intercept.

slugs.lmedata <- groupedData(pumps ~ state | slugs,
data=as.data.frame(data.slugs))
cs <- corAR1(form= ~ x|slugs)
res1 <- lme(pumps ~ state, random = ~1 | slugs, data=slugs.lmedata,
cor=cs)

---------------------------------

Now, I would like to add a complication to the model: The slugs were
observed in batches:
Batch_1 = {slugs_1, slugs_2, slugs_3}
Batch_2 = {slugs_5, slugs_6}
Batch_3 = {slugs_7, slugs_8, slugs_9, slugs_10}
Batch_4 = {slugs_11}
.
.
.
Batch_12 = {slugs_24, slugs_25}
Notice that there are 12 batches, and the number of slugs in each
batch differ, from 1 slug to 4 slugs.

I consider batch to be another random-effect on intercept. Thus I
fit model:

************pumps_tij = state_t + slugs_i + batch_ij + error_tij

Slugs, batch and error are normally distributed, but slugs and batch
are not nested factors.

I had fit following (however I'm not sure if that is right):

slugs.lmedataB <- groupedData(pumps ~ state | slugs/batch,
data=as.data.frame(data.slugs))
csB <- corAR1(form= ~ x|slugs/batch)
res1B <- lme(pumps ~ state, random = ~1 | slugs/batch,
data=slugs.lmedataB, cor=csB)

----------------------------------------

QUESTIONS:
1) Are my models right? Do I model the res1B model properly?

2) Until now, I have assumed that the number of pumps follow the
normal distribution. However I know that the variable pumps is
distributed along Poisson distribution. How can I model that? I
would like to use LOG or SQRT transformation, but I don't know how.

*****************************************

Thank you very much for all your help.

Martina Pavlicova

-- 
Department of Biostatistics
Columbia University
722 W. 168th Street, 6th floor
New York, NY 10032

Phone: (212) 305-9405
Fax: (212) 305-9408
Email: mp2370 at columbia.edu