Dear all,
I would like to know if you think that what I am doing is correct. I 
have always the same self-questionning when working with (generalized) 
linear mixed models with interactions.

First let me introduce you my data.
  A group of monkeys have been experimented within 3 different 
conditions (variable Nb with values 1, 2 and 3).
As a social group of monkeys, the rank among the group of each monkey is 
important in the way of explaining a reaction (variable Rk).
I want to see the effect of the condition, the rank and the interaction 
of these two variables on the outcome of the experiment (variable Hz).

There are 540 observations among 15 individuals (36 per individual)

The selected model :
- the variable Nb have been transformed into a factor
- the square root of the variable Hz was used to be able to use a LME 
instead of a GLMM (valided by a Q-Q plot)
- Hz seemed to be lineary dependent of the rank so I didn't used a 
factor version of Rk (hence Rk is numeric)

I only tried the random intercept based on the identity (ID) of each 
individual.

So the model I am using is :
model.sqrt<-lme(sqrtHz~Rk*Nb, random = ~ 1 | ID,data=data1)

Here is the summary of the fixed part of the model :
Fixed effects: sqrtHz ~ Rk * Nb

                    Value                Std.Error       DF       
t-value         p-value
(Intercept)  0.5031499      0.04741670    521     10.611238     0.0000
Rk             -0.0456437      0.00641704     13       -7.112890    0.0000
Nb2            0.0551886      0.02609983    521       2.114520     0.0349
Nb3            0.0026095      0.02609983    521      0.099982      0.9204
Rk:Nb2      -0.0004107     0.00353217    521     -0.116277     0.9075
Rk:Nb3      0.0144365      0.00353217     521      4.087143     0.0001

To obtain a p-value for the factor Nb and the intercation, I use :
 > anova(model.sqrt)
                       numDF denDF      F-value         p-value
(Intercept)        1          521       113.85619 <.0001
Rk                    1            13         45.33842 <.0001
Nb                     2          521        18.49245 <.0001
Rk:Nb                2         521        11.46233 <.0001


Hence both Nb and the interaction have significant effects.

My first question is do I have the right to do such a thing?

I want to compare the different experiment and the interactions.

First, I tried to use summary(glht(model.sqrt,linfct=mcp(Nb="Tukey"))) 
that give me the result :

                  Estimate    Std. Error z value   Pr(>|z|)
2 - 1 == 0  0.05519    0.02610     2.115      0.0868 .
3 - 1 == 0  0.00261    0.02610     0.100      0.9945
3 - 2 == 0 -0.05258    0.02610     -2.015     0.1088

This output came with the warning message : covariate interactions found 
-- default contrast might be inappropriate

Moreover I can not apply this to have multiple comparisons among the 
interactions.

Browsing the mailing list, I tried to follow the method given by :
http://thebiobucket.blogspot.com/2011/06/glmm-with-custom-multiple-comparisons.html

Hence I did :
 > c <- rbind("Condition 1 - Condition 2" = c(0,0,1,0,0,0),
+             "Condition 1 - Condition 3" = c(0,0,0,1,0,0),
+             "Condition 2 - Condition 3" = c(0,0,1,-1,0,0),
+             "Rank:Condition 1 - Rank:Condition 2" = c(0,0,0,0,1,0),
+              "Rank:Condition 1 - Rank:Condition 2" = c(0,0,0,0,0,1),
+              "Rank:Condition 1 - Rank:Condition 2" = c(0,0,0,0,1,-1)
+ )

and

 > summary(glht(model.sqrt,c))

          Simultaneous Tests for General Linear Hypotheses

Fit: lme.formula(fixed = sqrtHz ~ Rk * Nb, data = data1, random = ~1 |
     ID)

Linear Hypotheses:
                                                                        
   Estimate     Std. Error     z value   Pr(>|z|)
Condition 1 - Condition 2 == 0                     0.0551886  
0.0260998      2.115       0.131
Condition 1 - Condition 3 == 0                     0.0026095  
0.0260998       0.100      1.000
Condition 2 - Condition 3 == 0                     0.0525791  
0.0260998       2.015       0.162
Rank:Condition 1 - Rank:Condition 2 == 0 -0.0004107  0.0035322      
-0.116      1.000
Rank:Condition 1 - Rank:Condition 2 == 0  0.0144365  0.0035322       
4.087 <0.001 ***
Rank:Condition 1 - Rank:Condition 2 == 0 -0.0148472  0.0035322      
-4.203 <0.001 ***

So my real question is : is this correct?
If so, it's quite in contradiction with the results of 
summary(model.sqrt) and anova(model.sqrt), for the multiple comparisons 
among the different conditions. I know that the p-values of 
summary(model.sqrt) can be seen as Dunnett comparison with Condition 1.

Here I don't know what is the correction that is used but I think it 
takes into consideration the fact that I have done 6 different tests. 
These corrections seems to be too conservative as there is no 
significant differences among the different conditions.
Is there a way to select a less conservative correction?

If not I don't know how to explain the differences between the results 
of  both summary(model.sqrt) and anova(model.sqrt) and the last 
displayed one.

Thanks in advance for your answers, and appologies to have written such 
a long message.
Best,
Nicolas
-- 
Nicolas Poulin
Ingénieur de Recherche
Centre de Statistique de Strasbourg (CeStatS)
http://www-math.u-strasbg.fr/CeStatS/ 
<http://www-math.u-strasbg.fr/CentreStat/>
IRMA, UMR 7501
Université de Strasbourg et CNRS
7 rue René-Descartes
67084 Strasbourg Cedex

Tél. 33 (0)3 68 85 01 89

	[[alternative HTML version deleted]]


