R help - Nov 2004 - lme, two random effects, poisson distribution

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

Hello all,

I think that with the help of Mark Irwin, we solved the problem:

We fit the model using glmmPQL and instead of using variable,
'state', we model the independent fixed-effect 'state' as
I(x>10); i.e. it's 0 for resting time slots and 1 for
excited times slots.

Here is the code:
================slugs.lmedata <- groupedData(y ~ I(x>10) | slugs/batch,
data=as.data.frame(data.slugs))

csnew <- corAR1(form = ~ x | slugs/batch)

res.glmm <- glmmPQL(pumps ~ I(x>10), random = ~1|slugs/batch, family
= poisson, data=slugs.lmedata, correlation=csnew)

The summary of the model looks like this:
========================================> summary(res.glmm)
Linear mixed-effects model fit by maximum likelihood
 Data: slugs.lmedata
       AIC      BIC    logLik
  1265.934 1288.157 -626.9672

Random effects:
 Formula: ~1 | slugs
        (Intercept)
StdDev:   0.6959097

 Formula: ~1 | batch %in% slugs
        (Intercept) Residual
StdDev:    1.138034 1.304958

Correlation Structure: AR(1)
 Formula: ~x | slugs/batch
 Parameter estimate(s):
       Phi
-0.2820534
Variance function:
 Structure: fixed weights
 Formula: ~invwt
Fixed effects: y ~ I(x > 10)
                   Value Std.Error  DF   t-value p-value
(Intercept)   -0.5004731 0.2006282 142 -2.494531  0.0138
I(x > 10)TRUE  1.2553644 0.2494294 142  5.032944  0.0000
 Correlation:
              (Intr)
I(x > 10)TRUE -0.32

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-2.0856605 -0.5439464 -0.3606277  0.1855235  5.7606962

Number of Observations: 300
Number of Groups:
           slugs batch %in% slugs
              25              157

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

I think this is the right model. But I am interested very much in
your opinions. I do not use mixed-effect modelling very often.

Thank your for all your help.

Martina

Quoting "Martina Pavlicova, PhD" <mp2370 at columbia.edu>:
>
> 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
>