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
