[R] Problem on model simplification with glmmPQL

[Ronaldo Reis Jr.]

Hi all,

I try to make a split-plot with poisson errors using glmmPQL, but I
have some doubts about the model simplification.

Look my system:

Block = 3 blocks
Xvar1 = 2 levels
Xvar2 = 13 levels
Yvar = Count data Response 

I need know about the behaviour of Var1, Var2 and interaction
Var1:Var2.

Look the levels:
> levels(Xvar1)
[1] "A" "B"
> levels(Xvar2)
 [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m"
> levels(Block)
[1] "B1" "B2" "B3"

My general model in aov is:

aov(Yvar~Xvar1*Xvar2+Error(Block/Xvar1/Xvar2))

In glmmPQL I make this:

create the nested high level variable:

> BlockXvar1 <- as.factor(paste(Block,Xvar1))
> levels(BlockXvar1)
[1] "B1 A" "B1 B" "B2 A" "B2 B" "B3 A" "B3 B"
 
Make the model

> m1 <- 
glmmPQL(Yvar~Xvar1*Xvar2+Block,random=~1|BlockXvar1/Xvar2,family=poisson)
iteration 1 
iteration 2 
iteration 3 
iteration 4 
iteration 5 
iteration 6 
iteration 7 
iteration 8 
iteration 9 
iteration 10 
> anova(m1)
            numDF denDF  F-value p-value
(Intercept)     1    48 697.7770  <.0001
Xvar1           1     2 127.1378  0.0078
Xvar2          12    48 330.5172  <.0001
Block           2     2  22.7960  0.0420
Xvar1:Xvar2    12    48 282.0917  <.0001
 
OK, all variables are significative, now I try to make the model
simplification, I try to amalgamate levels of Xvar2.

Looking the mean:

> tapply(Yvar,list(Xvar2,Xvar1),mean)
           A          B
a  1.6666667  1.6666667
...
h  0.3333333  0.6666667
...
j  0.3333333  0.6666667
...
m 11.0000000  5.6666667

Look that the mean of h and j are exactly the same in A and B,
therefore h and j maybe amalgamated.

Then, I create a new Xvar2 whith these levels in one.

> Xvar2.amalg <- recode(Xvar2,"c('h','j')='hj'")
> levels(Xvar2.amalg)
 [1] "a"  "b"  "c"  "d"  "e"  "f"  "g"  "hj" "i"  "k"  "l"  "m" 

Now I make a new model, in lme I need to fit the model with
method="ML" to make it comparable, but in glmmPQL I dont know.

> m2 <- 
glmmPQL(Yvar~Xvar1*Xvar2.amalg+Block,random=~1|BlockXvar1/Xvar2.amalg,family=poisson)
iteration 1 
iteration 2 
iteration 3 
iteration 4 
iteration 5 
iteration 6 
iteration 7 
iteration 8 
iteration 9 
iteration 10 

and make a model comparison:

> anova(m1,m2)
   Model df      AIC      BIC     logLik   Test  L.Ratio p-value
m1     1 31 164.5596 237.6176  -51.27982                        
m2     2 29 353.6201 421.9647 -147.81006 1 vs 2 193.0605  <.0001

look the high significance, I dont understand this, the behaviour of
this two levels are the same, in this case the models must be no
significance and I get the model more simple, in this case teh model
with levels h and j like one level hj.

I try to make it using lme, but the comparison are significant too.

Where is the error?

Thanks
INte
Ronaldo

--
|   //|\\   [*****************************]
|| ( õ õ )  [Ronaldo Reis Júnior          ]
|     V     [ESALQ/USP-Entomologia, CP-09 ]
||  / l \   [13418-900 Piracicaba - SP    ]
|  /(lin)\  [Fone: 19-429-4199 r.229      ]
||/(linux)\ [chrysopa at insecta.ufv.br      ]
|/ (linux) \[ICQ#: 5692561                ]
||  ( x )   [*****************************]
||| _/ \_ Powered by Gnu/Debian Woody
-----------------------------------
Insecta - Entomologia
Departamento de Biologia Animal
Universidade Federal de Viçosa