[R] Re: lme, two random effects, poisson distribution

[Martina Pavlicova, PhD]

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
>