Skip to content

GLM: Difference between treatment groups for each colony (Tukey posthoc test)

5 messages · Sophie Waegebaert, paul debes, Peter Claussen

Original
Sophie Waegebaert
# 
Hi,

I'm still learning how to use R and I have some trouble making using Tukey
posthoc tests. I have a dataset with 3 colonies (A, B and C). Each colony
is divided into 2 treatments: control and DWV. I want to run a GLM to test
wether there is a difference in life expectancy (last.scan) between the
treatment groups for each of the colonies, but I do not know if I am using
the right strategy.

I have taken 'treatment' as a fixed factor and 'colony' as a random factor:

fit_life = lmer(last.scan~treatment + (1|colony), data =
data)Anova(fit_life, type = 3) # Type of treatment has a significant
effect on on the life expectancy.
    Response: last.scan
          Chisq Df Pr(>Chisq)
    (Intercept) 106.976  1  < 2.2e-16 ***
    treatment    25.373  1  4.724e-07 ***

And this is the code I use to do a Tukey posthoc test:

mcp = glht(fit_life, linfct = mcp(treatment = "Tukey"))
summary(mcp)# DWV treatment significantly changes life expectancy (z =
-9.734, p = < 2e-16)

Is it possible to find the difference for each colony?

Thanks a lot for an explanation or hint!

Cheers,

Sophie
paul debes
# 
Hi Sophie and List,

If you are interested in the treatment contrasts for each level of the  
Colony factor (with three levels), it may be better to view your  
experiment as a factorial and specify a simple linear model rather than  
the presently specified linear mixed model. A random Colony term with only  
three levels may also not provide the best estimate for the correlation of  
colony data between treatments (i.e., as intraclass correlation) that is  
taken into account when testing the general treatment term.

Would this work for you?:

fit_life = lm(last.scan ~ 1 + treatment*colony, data = data)

In case the treatment:colony term is significant you could conduct the  
additional pairwise tests of interest.

Best,
Paul



On Mon, 07 Dec 2015 13:02:36 +0200, Sophie Waegebaert
<sophie.waegebaert at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Sophie Waegebaert
    

    

      
      #
    

  


  

      
Hi Paul,

Thank you very much for the hint!

The treatment:colony term is not significant, so I will not include it in
the model.

However, I get the same results for the Tukey posthoc test as I had before.
I am able to say that the DWV treatment significantly reduces the life
expectancy, but I do not know the difference between the colonies. With the
results of the Tukey test I want to make a figure (see attachment), but I
do not know how to do it. The asterisks are based on Tukey posthoc test.

Do you know a way to do this in R?

Thanks a lot!

Cheers,
Sophie



2015-12-07 12:30 GMT+01:00 paul debes <paul.debes at utu.fi>:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      paul debes
    

    

      
      #
    

  


  

      
Hi Sophie,

If the treatment:colony term is NS, what about the colony term?

If the colony term is NS, you could also remove the colony term from the  
model and only leave the treatment term in the model (no need for a  
post-hoc test for the treatment term: it has only two levels: F-test/LRT  
results are sufficient).

If the colony term is also important (but the interaction is NS), the  
treatment effects are the same for each level of colony, and the colony  
effects differ from each other similarly within each treatment level. To  
find what colony levels differ from each other you will need to specify  
the colony term for your post-hoc test.
You can plot the means in two ways, depending on what is your study focus:  
1) two means for the two treatments averaged across colonies AND/OR 3  
means for the three colonies averaged across treatments (in 2 figures). 2)  
six means of the three colonies within each treatment, or vice versa (1  
figure).

To predict model means for many model types, I personally like using the  
library "lsmeans".
This might work (use "summary" to facilitate getting the results as  
data.frame for easier plotting):

library(lsmeans)
model.means.treatment = data.frame(summary(lsmeans(model, ~ treatment)))
model.means.colony = data.frame(summary(lsmeans(model, ~ colony)))
model.means.inter = data.frame(summary(lsmeans(model, ~ treatment:colony)))

Hope this helps,
Paul


On Mon, 07 Dec 2015 14:12:37 +0200, Sophie Waegebaert

        
<sophie.waegebaert at gmail.com> wrote:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Peter Claussen
    

    

      
      #
    

  


  

      
Sophie,

As I read this, you are mixing analysis.

You?ve specified ?treatment? as fixed and ?colony? as random. In this analysis, it is reasonable to perform posthoc comparisons on treatment, but not on colony. Colony is a random variable and each colony is simply drawn from a larger pool, one just as likely as the other. 

The Anova table based on a mixed model only reports fixed effects, and that is proper. If you want to get an estimate of how much variation is associated with colonies, use 

summary(fit_life)
or
VarCorr(fit_life)

This will report the variance associated with colonies, and you might judge from the magnitude of this term how different are the colonies - how variable is your population of colonies. The caveat is that with only 3 colonies, you might not get an accurate estimate.

Peter