I have been trying to figure out how to do a fairly basic repeated measures analysis using linear mixed effects in R, and then analysing it using post-hoc tests. The problem is that I'm not sure whether the output I get is statistically sound?

The response variable: weighted- an index of habitat preference (prop. individuals on habitatA / prop. of total habitat that is A). A value above 1 indicates the habitat is being used more than what you would expect from its availability. this was repeatedly measured on the same colony through time over several weeks

Fixed variables: Type - habitat type (live/dead), weeks - the time variable

Random variables: colony - because each measurement of colony violates independence assumption.

Here's what the data loss like plotted over time (orange=live habitat, blue=dead habitat):

enter image description here

i run the analysis using the lmer() function from the lme4 package:

results_full=lmer(weighted~type*weeks+(weeks|colony), data=Pos, REML=F)

My reasoning is that i have no reason to expect a random intercept, they should all start on 1 at time 0, and then individuals will start avoiding the dead habitat and favouring the live habitat. The (weeks|colony) term allows the slope of each colony to be random across time?

So to my question:

I compare the likelihood of two models with each other, in a likelihood ratio test to get p-values of the fixed effects using a reduced model:

results_null=lmer(weighted~type+weeks+ (weeks|colony), data=Pos, REML=F)
anova(results_null, results_full)

But what I'm really interested in is at what time point (week) do the individuals start avoiding the dead habitat. as you can see from the figure this happens at week 1 so comparing live-dead habitat week by week "should" generate a n/s result at week 0 and sig result from then on (I'm not trying to force a statistically significant result, but the fig is pretty clear...)

I tried converting the weeks into a factor, and then performing

lsmeans(results_full, pairwise~type+weeks)

But it didn't generate anything that seemed meaningful, the output didn't make sense in relation to the data.

Does anyone have any thoughts on A) whether my model and test is appropriate to this data, and B) how I can perform a post hoc test to compare habitat type over time?

Would it be appropriate to use a Dunetts post hoc test to compare preferences to a reference value (=1) rather than to each other?

Grateful for any ideas or pointers!

  • $\begingroup$ (Comment to voters:) Issues like "The problem is that I'm not sure whether the output I get is statistically sound" clearly indicate that this is not simply a coding question. (If you disagree, could you please explain to me how that doesn't require statistical expertise to answer?) $\endgroup$
    – Glen_b
    Commented May 7, 2015 at 10:16
  • $\begingroup$ It appears that, except for week 2, the two weighted values in your plot are dependent. Are they measured independently for the two types, or is weighted some sort of measure of the balance between the two types? If the latter, your model does not make sense. $\endgroup$
    – Russ Lenth
    Commented May 8, 2015 at 18:34
  • $\begingroup$ @Glen_b , not sure what you mean? Isn't CV a place to ask statistical questions? The tone of your comment suggests I've committed some kind of CV faux pas? I'm eternally grateful for help I've received through all the Stack Exchange sites, but we can't all be experts in the etiquette of these sites so please let me know if I've somehow asked my question the wrong way? :) $\endgroup$ Commented May 9, 2015 at 2:37
  • $\begingroup$ @rvl they are a measure of the balance between #fish on dead vs live habitat on a coral colony, but being a mean ratio of two proportions (prop. fish on the habitat/prop. of that habitat available) the index isn't perfectly matched. How do you mean the model doesn't make sense? $\endgroup$ Commented May 9, 2015 at 2:42
  • 1
    $\begingroup$ I think Glen_b is commenting to whoever thinks your question should be closed, and is saying your question DOES have merit. On my question, I was trying to say that if the measures for the two types are completely dependent on each other, then a model for just one of the types could be used to predict results for the other type. But apparently that isn't the situation. I'm thinking though that my earlier answer doesn't get at your question very well. $\endgroup$
    – Russ Lenth
    Commented May 9, 2015 at 2:48

1 Answer 1


It seems to me that your lsmeans statement asks for pairwise comparisons of all combinations of the two factors, which is probably more than you need, and potentially confusing for that reason. A specification like pairwise ~ weeks | type would do pairwise comparisons of the weeks, for each type. But If you're looking specifically at changes, I think a smaller number of custom contrasts would be preferable.

First, I suggest getting the LS means alone in an object, and doing the contrasts separately. This isn't absolutely necessary, it's just cleaner:

lsm = lsmeans(results_full, ~ weeks | type)

I'm assuming there are six different weeks -- if not, you need to adapt the following. One possibility is to test the changes between two consecutive weeks:

jumps = list(j12 = c(-1,1,0,0,0,0), j23 = c(0,-1,1,0,0,0), 
             j34 = c(0,0,-1,1,0,0), j45 = c(0,0,0,-1,1,0),
             j46 = c(0,0,0,0,-1,1))
contrast(lsm, jumps, adjust = "mvt")

The "mvt" adjustment adjusts the P values appropriately for picking out the maximum of these contrasts -- i.e., finding where is the biggest change. (You need to have the mvtnorm package installed; and be warned that the adjusted P values are obtained by simulation, so only the first 3 or so digits are trustworthy.)

An alternative is to compare the average of all data before with the average of all data after:

mean.chg = list(m12 = c(-1, 1/5,1/5,1/5,1/5,1/5),
                m23 = c(-1/2,-1/2, 1/4,1/4,1/4,1/4),
                m34 = c(-1/3,-1/3,-1/3, 1/3,1/3,1/3),
                m45 = c(-1/4,-1/4,-1/4,-1/4, 1/2,1/2),
                m56 = c(-1/5,-1/5,-1/5,-1/5,-1/5, 1))
contrast(lsm, mean.chg, adjust = "mvt")

You probably should pick one of these approaches and go with it, so you're not accused of shopping for the results you want. Or, combine the lists, and the mvt adjustment will penalize you appropriately for testing 10 contrasts instead of 5.

PS - I don't think I made any typos in the untested code above, but it's possible, so be careful.

  • $\begingroup$ Thanks @rvl , this is great, thanks for the time. After having a play with this I'm thinking my last suggestion is more correct, i.e. testing each preference per week against a reference value (=1, no preference). Do you know of a test which would allow me to do that? Dunnets does this for repeated measures ANOVA, but I'm not sure if it's correct to use for these kinds of data? Thanks agin for your help! $\endgroup$ Commented May 9, 2015 at 2:46
  • $\begingroup$ Hmmm, with lsm as defined above, try test(lsm, null =1) $\endgroup$
    – Russ Lenth
    Commented May 9, 2015 at 2:52
  • $\begingroup$ Thanks so much @rvl, this solved my problem. With a bit of tweaking I got your suggested method to do both of the analyses I requested in the original question: lsm = lsmeans(results_full, pairwise~ type|weeks) lsm tests the two types agains each other by week (week=factor). Secondly test(lsm, null =1) tested if each value is different from 1. Thank you so much! $\endgroup$ Commented May 9, 2015 at 4:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.