Skip to content

lme, two random effects, poisson distribution

2 messages · Martina Pavlicova, PhD

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
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:
=========================================
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>: