I'm dealing with a dataset that comprises ~30 groups of observations. In each group, we measured independent and dependent variables a certain number of times. If I plot the relationship between independent and dependent variables for each group separately, I get something like this (showing just 4 of the groups here, with the independent variable on the X and dependent on the Y):

enter image description here

My ultimate goal is to determine whether the effect of the independent variable (i.e., the slope of a linear regression dep ~ ind) is significantly different between groups, and to quantify which groups have the largest/smallest effects.

To do this, I've put every data point into a single dataframe, with a column indicating the group it came from. Then, I ran the following model, calculating an interaction term for group.

lm1 = lm(dep ~ ind * group, data=df)

To figure out whether a model that accounts for differences between groups is a better fit than a model that doesn't, I also fit the following model:

lm2 = lm(dep ~ ind, data=df)

And then I compared the two models using ANOVA:

anova(lm1, lm2)

The p-value of the ANOVA is about 1e-5, suggesting that the model that takes group into account is the better fit.

However, I've also done the following to figure out if effects in particular groups are larger/smaller than in others:

library(lsmeans) lst <- lstrends(lm1, "group", var="ind") pairs(lst)

Here, I'm performing pair-wise comparisons of effects between all groups, and after correcting for multiple testing, none of the pairwise comparisons are significantly different.

To me, this suggests that overall, a model that takes group into account is a better fit. But because the number of data points in each group is small (~8-10), the per-group regressions produce large confidence intervals, making it difficult to find significant differences between groups.

Am I thinking about this correctly?

(Apologies if a similar question has been asked before, I couldn't find any questions that directly addressed my concerns.)

| cite | improve this question | | | | |
  • $\begingroup$ I don't know what lsmeans does but why are you not looking at the lm1 output? The linear model will give you the p-values for all the tested hypotheses. $\endgroup$ – user2974951 Sep 17 '18 at 6:28
  • $\begingroup$ It is possible for an overall effect to be significant and yet have no significant pairwise comparisons. However, another aspect of this is an adjustment is made to the P values for multiplicity, and it is pretty drastic when protecting against making any type I errors among those 435 comparisons, as is done with the default adjustment method (Tukey). $\endgroup$ – Russ Lenth Sep 17 '18 at 12:35
  • $\begingroup$ PS As of Oct 10, lsmeans will be archived on CRAN. Use emmeans instead. $\endgroup$ – Russ Lenth Sep 17 '18 at 12:40

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.