2
$\begingroup$

I am analysing a mixed effect model:

library(lme4)
lmm_treat <- lmer(total ~ TREAT * POS + (1|NO_UNIT), data=data)

DESCRIPTION OF THE MODEL: "total" is the value I am comparing between the treatments (TREAT) and temperatures (POS). TREAT is used as a factor and POS as numeric (and I cannot change it to factor due to problems with the degrees of freedom).

AIM: I want to test whether the value "total" is significantly different between treatments (TREAT) and temperatures (POS) and the interaction between them.

H0: The total value is not different between treatments (TREAT), temperatures (POS) and the interaction between them

ANALYSIS OF THE MODEL: I analyse the model using Anova type II. When the differences between TREAT are significant I run a Tukey-Kramer post hoc test.

library(car)

Anova(lmm_treat)

#to analyse the differences between treatments I do this:

library(ghlt)
summary(glht(lmm_treat,linfct=mcp(TREAT="Tukey")))

When the interaction between POS:TREAT is significant I do a bootstrap analysis to see which TREATs have significantly different interactions. As TREAT is a factor, I have to change the baseline and repeat the fit of the model everytime to have the results complete.

library(boot)

#check the levels of TREAT

levels(data$TREAT)

[1] "BA5" "BNA" "NBA1" "NBA5" "NBNA"

mySumm <- function(.) { s <- sigma(.)
c(beta =getME(., "beta"), sigma = s, sig01 = unname(s * getME(., "theta"))) }

# run bootstrap analysis for calculation of confidence intervals of parameter estimates
mod_lmm1_boot <- bootMer(lmm_treat,mySumm, nsim=300)
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=2)
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=3)
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=4)
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=5)
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=6)
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=7)
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=8)
boot.ci(mod_lmm1_boot,type="perc",conf=.95,index=9)

#then I put other TREAT as the baseline and I repeat the method
data$TREAT <- relevel(data$TREAT, "NBNA")

TWO QUESTIONS:

  1. Is it correct to use Tukey not including the TREAT:POS interaction as it was defined in the model? If not, which other analysis is better?
  2. Is there any other way of analysing the interaction that is not bootstrapping in the way I am doing it? I have the feeling that changing the baseline every time to have a full analysis is not the best....

EDIT: I added some data in case somebody wants to try out something

data<-structure(list(POS = c(1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 
    4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 
    5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 
    1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5), 
        TREAT = structure(c(6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
        6L, 6L, 6L, 6L, 6L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
        4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
        5L, 5L, 5L, 5L, 5L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
        3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
        2L, 2L, 2L, 2L, 2L), .Label = c("BA1", "BA5", "BNA", "NBA1", 
        "NBA5", "NBNA"), class = "factor"), total = c(45157490.15, 
        30714210, 37549097.17, 26735871.89, 39617487.82, 47906908.53, 
        48642615.88, 33871012.13, 26352389.87, 25953030.37, 37550336.85, 
        46137989.05, 34008559.85, 23769345.12, 30966059.08, 59472042.02, 
        61960980.86, 40160217.98, 36441957.94, 28162593.61, 54707657.67, 
        50555556.05, 42669534.21, 35217029.22, 27028853.79, 43400224.72, 
        50721934.09, 34940532.34, 32281963.49, 33546024.53, 30752701.93, 
        40799531.72, 38421451.33, 28135608.63, 40474465.41, 30164687.1, 
        32075304.33, 24093615.88, 25185393.21, 32423802.51, 31883418.58, 
        39548752.26, 45952327.77, 37053647.85, 41470953.52, 27692166.81, 
        71865066.04, 55702824.79, 47587488.02, 40221235.39, 38034718.14, 
        37895296.57, 27487907.01, 34127936.28, 45487761.41, 33154895.21, 
        43223297.18, 45125276.85, 47481387.34, 59212779.4, 8576024.475, 
        47314612.41, 27390001.83, 31459563.39, 37896690.43, 31734599.76, 
        55331352.05, 54359675.18, 45311774.87, 49534712.91, 11681184.43, 
        59984689.01, 51226923.15, 32707623.96, 28300471.74), NO_UNIT = structure(c(1L, 
        1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 
        4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 
        7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 10L, 
        10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 
        12L, 12L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 14L, 
        15L, 15L, 15L, 15L, 15L), .Label = c("1", "2", "3", "4", 
        "5", "6", "7", "8", "9", "10", "11", "12", "16", "17", "18"
        ), class = "factor")), .Names = c("POS", "TREAT", "total", 
    "NO_UNIT"), row.names = c(NA, -75L), class = "data.frame")
$\endgroup$
  • $\begingroup$ I do not use the Tukey test anymore, instead I just use the bootstrapping to see all the differences; however changing the baseline everytime does not look like the best option... $\endgroup$ – kumbu Mar 1 '17 at 8:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.