you can use the glht function in package multcomp, as here:

summary(glht(YOUR MODEL, linfct=mcp(YOUR FIXED FACTOR="Tukey")))

what is written in capitals should be replaced by your object names!!

See also here:

https://stat.ethz.ch/pipermail/r-help/2012-January/299623.html

Cheers, Ivan

Dear Simon,

I am not used to work with the R software, but I think you can use HSDTukey.

Furthermore, this file may help you.

Regards.

Thanks Jurandy, I tried with Tukey, but this test cannot be applied to an object of class "lme".

(Thanks for pdf)

package agricolae in R has many post hoc test , maybe it can be useful for you

you can use the glht function in package multcomp, as here:

summary(glht(YOUR MODEL, linfct=mcp(YOUR FIXED FACTOR="Tukey")))

what is written in capitals should be replaced by your object names!!

See also here:

https://stat.ethz.ch/pipermail/r-help/2012-January/299623.html

Cheers, Ivan

Ivan, apparently glht function is deprecated in R version 2.15.1. I'll try with agricolae pack. thanks!

Post Hoc tests are just different ways to adjust p-value regarding the number of comparisons performed. So, if you have two factors and only one is significant (I assume that there is no significant interaction either), you actually have four groups, one to each level of your significant factor. You can calculate, now, the means of your four groups (each group is formed by all individuals with that level, regardless the other factor).

Now, calculate all difference between groups (there are six differences). Which one are significant? All difference that is greater than the Honestly Significant Difference (HSD), which can be calculated as:

t(a/2) * sqrt(MSR / n),

where t(a/2) is the quantile of the Tukey distribution, with "m" means and "df" degrees of freedom, MSR is the residual mean square from ANOVA table and n is the size of the groups. How to do that at R????

qtukey(0.95, nmeans = m, df = df) * sqrt(MSR / n)

Replace m, df, MSR and n properly...

If you have different group sizes, replace "n" by 2*n1*n2 / (n1+n2), where n1 and n2 are the sizes of the two groups under evaluation. In this case, you have the Tukey-Kramer procedure, suitable to comparisons of groups with different sizes. In R,

(2*length(G1)*length(G2)/(length(G1)+length(G2)))

Easier than look for a package... =D

I agree the easiest way is using TukeyHSD or even the pairwise.test functions. However, for me, I get stuck using R when I want to extract the LSMEANS for the fixed effects. Does anyone has an idea how to get this? Thanks

You can use lsmeans function (library lsmeans).

Try

lsmeans(YOUR MODEL, pairwise~FIXED FACTOR, adjust="tukey")

If you want to test the interaction try

lsmeans(YOUR MODEL, pairwise~FACTOR1*FACTOR2, adjust="tukey")

Normally I solve it with glht and some little trick for interactions. For example with both solutions you can compute the differences by the following way.

But, on the other hand not is easy to find a good tutorial about this kind of statistics. Under my own point of view, you must interpret this results carefully

Kind regards and I hope that you enjoy the example

Javier

############################################

require(LMERConvenienceFunctions)

require(multcomp)

require(lsmeans)

data

Thank you, Carlos Lara-Romero. I have found: lsmeans(YOUR MODEL, pairwise~FACTOR1*FACTOR2, adjust="tukey") useful for my work.

You're welcome!

Carlos Lara-Romero, thank you for that code, it has helped me in my analysis

Great!

Hi Carlos,

Is there any trick to install the library lsmeans? I have tried hardly but I can´t

Thank you!

Hi Natalia, 

I've had no problems with the library.

This problems could happen due to incompatibilities between versions. I am working with R version 3.0.3 and lsmeans version 2.05. 

I installed the library using RStudio. Maybe, you could try to use the same software for managing the library.

Good luck!

I have installed R 3.0.3 and it works!!

Thank you so much

Dear all,

Are you sure that using LSmean to compute post-hoc for a lmer model takes into account the random effect? if, not, I presume it is unecessary to use lmer but instead lm...?

thanks a lot

carole

Hi Carole,

I guess the answer is yes as the post-hoc comparison using lsmeans is calculated based on LMM model as you can see in the above replies.

Cheers,

Patchanok

Hi Carole,

lsmean does account for the random effect in the multicomp (however, you have to be sure your model is fitted using ML and not REML which would provide false results). PLus, if I may say, that's the nature of your sampling design (or your data) that drives the fact that you'll do lm or lmer (eg if your sampling desing is nested) not the possibility to do post-hoc test or not.

Another solution is to do bootstrap or MCMCsamp to estimate the SE, but this will not result in p-values and the decision for marginal results will be all yours...

Hope this helps.

Hi Yoan,

Thank you for your message. Could I ask if you have a reference for the using of ML with lsmean? It will be useful for my writing. I cannot find it anywhere.

Take care,

Patchanok

 thank you so much. indeed i have use glmer for my data. I'm currently searching how to do contrast comparisons. I have used mcp package but undortunately I have an interaction between a continuous and a categorical variables. I'm not sure that contrast is relevant when interaction are presents (I have a nice message "

"covariate interactions found -- default contrast might be inappropriate

2: In RET$pfunction("adjusted", ...) : Completion with error > abseps

3: In RET$pfunction("adjusted", ...) : Completion with error > abseps'.

I have read a lot of information about contrast possibilities and actually I have any idea how to perform them correctly regarding my design (repeated design with crossed effects).

thanks a lot

carole

I hope that this can help: summary(glht(model_lmer4, linfct=mcp(a="Tukey", interaction_average=TRUE)))

Hi Federico, thanks a lot for your quick response.. I have tried this but I obtain only simple effect and some kindly messages :

In mcp2matrix2(model, linfct = linfct, interaction_average = ia, :

covariate interactions found -- please choose appropriate contrast

2: In RET$pfunction("adjusted", ...) : Completion with error > abseps

3: In RET$pfunction("adjusted", ...) : Completion with error > abseps

the interaction variable is a continuous one, not a categorical. i'am wonering f this is not a problem here..

do you have any idea about that??

thanks

carole

Hi Carole And Federico,

Carole, if you have such interaction (factor vs. continuous variable), it does not make sense to do a multicomp as, what you have, is actually n equations (one for each level of the factor) and the significance of the interaction does already tell you whether the slope is significantly influenced by your factor or not.

Federico, I do not have a reference, I might have read it in the help of lsmeans but could not find it again (it's monday morning though...). There are some elements (but dealing with AIC model comparisons) in the MuMIn help (http://cran.r-project.org/web/packages/MuMIn/MuMIn.pdf, p.5 for example) and in the AICcmodavg help also (http://cran.r-project.org/web/packages/AICcmodavg/AICcmodavg.pdf  p. 10). 

Hope this helps.

HI yoan

thank you, I agree with that, it means that I couldn't go further? I also thought to LSMEANS because pairwaise comparisons are important for my hypothesis.

best

carole

For me, it means that you have all the answers you need in the summary of your model (except maybe if your factor comprises more thant 2 levels). In this latter case, you should try to represent your equations to have the trends, and comment this, or fix your continuous variable and compare means and se.

Best, Yoan.

Indeed, I have one categorical with two levels but another one with 4 levels. What do you mean by "fix" the continuous variable? thanks a lot

Hi Yoan, thank you for the reply, re the reference for lsmeans. :-)

Hi Patchanok,

I think that you don't need an specific reference. You can just cite lsmeans as Yoan suggested. I agree with Yoan it is a good approach representing the levels of the categorical variable (i.e. a slope per level) versus the continous variable.

Hi to all,

lsmeans is a "standard" You don't need a specific reference. But, does anyone know any way to compute the effect size of this test? A referee ask me for it and I don't know how to compute it.

Any suggestion?

Cheers

In the lsmeans output, it provides you mean and SE. In my study, I use Cohen's d effect size which requires mean, SD and number of subjects in each group. So I change SE to SD in Excel, then get the effect size based on pairwise comparison computed based on LMM.

I am also analysing data from mix modelling using lme function. I have encountered a problem like when I used 'method=ML' and 'method=REML'. the results produced different significance for the treatments during pairwise comparison by lsmeans function. Why this difference? and which method is best? anyone can suggest me!

Hi,

I read somewhere (should be in the lsmeans help?) that for multicomparisons you should use the "ML". I have in mind that with "REML" there are uncertainities about the number of degrees of freedom but cannot assert that for sure.

Otherwise, to know how to cite a package, try citation("package") ! ^^

Cheers.

The book by Pinheiro and Bates (Mixed-effects models in S and S-PLUS) discusses multiple comparisons in the LME model.  I don't have my copy handy, but Yoan Paillet's comment that the model should be fit using ML if you want to perform multiple comparisons is correct if I am remembering Pinheiro and Bates.

Hi Yoan,

For mix effect model in R, you should use method 'ML' for model comparison, BUT you should use method 'REML' during multiple comparison using lsmeans.

Hi Khagendra,

could you develop your answer a bit as it is not fully clear to me : you mean generalized (not mixed) models should use ML (which make sense if you do not have a random component) but that lsmeans uses REML, but for which kind of model ? Fixed only ?

That does not make sense to me as in this case lsmeans only performs a pure tukey (or other type of multiple comparison test) and your glm is equivalent to an anova or a lm if your error disctribution is Gaussian.

Any reference and/or reason to do so ?

Thanks in advance !

My clarification was only for linear mixed effect model , where you used 'lme' function under 'nlme' package, and not for GLM and ML. However, I can say  'ML' is biased for the estimation of variance components, but if you have larger sample sizes then the bias gets smaller. 

Hi,

You can follow the steps below, using multcomp:

> library(multcomp)

> com=glht(your model,linfct=mcp(you fixed factor="Tukey"))

> cld(com)

Cheers.

Hi Yoan and Khagendra,

I am quite confused with your answers about using ML or REML. I used  lmer fit by REML like: mymodel=lmer(A~B*C+(1|D), where D is a random factor. Then, I used lsmeans  to compare for each factor and their interaction which levels are significant, for ex.: lsmeans(mymodel, pairwise~A,ajust="tukey"). So, does the results gave by lsmeans (using the lmer model fitting with REML) is right ?

Correction:

Hi Simón P. Castillo ,

lsmeans seems to be a correct way to perform pairwase interractions.

But don't forget to correct your P-values accoring to the interractions you get.

Also, if the FACTOR1*FACTOR2 is not significant, you shoud put a bar | instead of a star * like this:

lsmeans(model, pairwise ~ FACTOR1 | FACTOR2, adjust = "tukey").

Cheers

But what about: lsmeans(model, pairwise ~ FACTOR1 * FACTOR2 | DAY, adjust =  "tukey") where DAY is continuous variable? Should DAY be treated as factor to get multiple comparison for different days (lets say from DAY 1 to DAY 10?).

You might try the Dunnett's test or running the longitudinal data analysis, depending on your data and purposes

Plus, if your treat day like a factor, you ignore the fact that your observations are ordinated (day one comes before day 2) and non-independant... This should first be integrated in the model (e.g. in the variance-covariance structure) before trying to perform post-hoc tests. And even with that, I don't even know post-hoc tests would be meaningful, because your model already tests whether the effect of "day" is significant (when it's coded as a continuous variable).

I need both: I want to see the effect of day on the response variable, plus I also want to day as a repeated variance on which variance-covariance structure will be built. I make new variable called DAYFAC to know the effects of day (how it differs from DAY1 to DAY10) and I use day as a repeated measures for variance-covariance structure! Is not it right strategy?

summary(glht(YOUR MODEL, linfct=mcp(YOUR FIXED FACTOR="Tukey")))

If I have two fixed factors，and I want to get their post hoc of  interaction, how to modify the code?

Hi,

You can try this : lsmeans(model, pairwise ~ FACTOR1 * FACTOR2 | DAY, at = list(DAY = c(1:10)), adjust = "tukey")

Best regards

Hi,

you can try with this: summary(glht(model, lsm(pairwise ~ Factor1|Factor2)))

Otherwise with package phia:

testInteractions(model, fixed="Factor1", across="Factpor2")

testInteractions(model, pairwise=c("Factor1","Factor2"), idata=data_frame)

Ciao!

I have this model:

fit6

Hi Megan Earney, you can use cld(lsmeans(fit6, pairwise ~ beach | method)) for compact letter display of pairwise comparisons in lsmeans package

hello,

does someone know the difference between Anova type II, Anova type III following a lmer model? Indeed, the results are very different regarding the significance of a variable or interactions. I use Reml (small sample size : 148 subjects with 5 repetitions for each) with kenward rogers df approximation. my fixed variables are all categorical and there is a three way interaction. I only use random intercept. thanks a lot :-)

Hi Fantini Carole,

You may find some responses here:

https://stats.stackexchange.com/questions/60362/choice-between-type-i-type-ii-or-type-iii-anova

https://stats.stackexchange.com/questions/223626/r-anova-vs-anova-for-test-of-categorical-predictor-from-glmer-or-glm-nb-ob

Best regards

TukeyHSD(aov(yourmodel))

provides: contrast(diff) /lwr/upr/p adj.

I'd use lsmeans way though

The package emmeans: install.packages("emmeans") supports lme/lmer models.

name

Hi Dear all,

In addation to the previous answers, I would like to add some codes that allow us to get post hoc test for an interaction term as well as the classification of its different levels .

library(lmer)

library(lemerTest)

library(lsmeans) # for lsmeans () and lsm()

library(multcomp) # for glht () and cld()

> model summary(glht(model, linfct = mcp(FIXED FACTOR ="Tukey")))

# -- 1) Least square difference (LSD) for the interaction "FACTOR 1:FACTOR 2":

> difflsmeans(model, test.effs = "FACTOR 1:FACTOR 2")

# -- 2) Using Turkey adjustement:

> lsmo = means::lsmeans(model, pairwise ~FACTOR 1|FACTOR 2, adjust = "Tukey")

> summary(as.glht(pairs(lmer.lsmo), by=NULL))

# or alternatively:

> summary(glht(model, lsm(pairwise ~ FACTOR 1|FACTOR 2), by=NULL))

# -- 3) Clasification of means by factor's levels:

# A single factor

> tuk.glht cld(tuk.glht, decreasing=T)

# An interaction term

> lsmo cld(lsmo, decreasing =T, by =NULL)

Best wishes !

As someone has reopened this topic I think for most applications I'd recommend the emmeans package.

https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html

List of vignette topics:

https://cran.r-project.org/web/packages/emmeans/vignettes/vignette-topics.html