I have a dataset where my Predictor variable (IV) does a reasonable job at predicting the response variable (DV). I have a second IV that is a grouping variable. In this plot there seems to be some grouping going on but I am not to sure if I am correct. I have also plotted a big circle which is the centroid for each group.

My question is 2 fold. How can I be certain that adding this grouping feature improves my model? Secondly, if the grouping is significant, how do I go about including it in my model (I program in R)?

Thank you kindly for your time.

enter image description here

Below is what I think @Erik is trying to tell me to do. Erick I ran a partial pool model: I scaled and centered data here which is why the X axis is different for the 2 graphs. I also added a no-pooling and complete pooling model to the graph

mod = lmer(Y ~ X + (X|Inudstries)) 

enter image description here

Would you do the anova(lm(Y ~ X), lmer(Y ~ X + (X|Inudstries))) To see which model is better?

  • $\begingroup$ (i) What kind of model do you use? Is this linear regression? (ii) What does the grouping represent? $\endgroup$
    – Roland
    Dec 17, 2019 at 10:04
  • $\begingroup$ @Roland this a simple regression Y~ X. I then coloured by a grouping factor (Industry). This is financial data but unfortunately can not go deeper then that. Does it look like grouping by industry adds value? $\endgroup$ Dec 17, 2019 at 14:54
  • 1
    $\begingroup$ You can test that easily: anova(lm(y ~ x, data = DF), lm(y ~ x * g, data = DF)) $\endgroup$
    – Roland
    Dec 17, 2019 at 15:00
  • $\begingroup$ Thanks Roland. The 2 models are statistically different with the RSS model 1 = .102 and RSS model 2 = .900. However, do you think a linear model is the best choice here? For future prediction, would it be better to shrink to the centroid? If so, how would you go about that? I was thinking using Partial Pooling/Multilevel Modeling $\endgroup$ Dec 17, 2019 at 15:32
  • $\begingroup$ How many industries are represented in your data? $\endgroup$
    – Erik Ruzek
    Dec 17, 2019 at 21:09

1 Answer 1


First, given the data you plotted, your groups look to have quite different means on your outcome. If the rest of your data is similar, then I can guarantee that a model that somehow accounts for groups is going to be a better fit to your data than one that does not do so. You have a couple of options with 15 groups:

  1. Fit a multilevel model that allows each group to have its own intercept, but that imposes a distribution on the values of the intercepts (normal, mean of 0 with an estimated variance). In so doing, this is also referred to as a random intercept model (or a partial pooling model). This is handy if you have different numbers of observations for each group as the intercept for groups with fewer data points will be pulled back toward the overall average outcome for the entire dataset. Since you have only 15 groups, you want to make sure and estimate your model using REML (restricted maximum likelihood), which is the default in lme4().
  2. You could enter your grouping variable as a factor in a OLS model. This is sometimes called a fixed effect model (or no-pooling model). This model also estimates an intercept for each group, but does not adjust the intercept based on the amount of data for a given group. Hypothetically, if a group had only one data point, the value of that group's intercept in this model would be that data point. In contrast, the multilevel model would almost certainly pull that group's intercept value back to the mean value on the outcome in the data.

In terms of a formal test of whether the partial pooling (multilevel) model is a better fit to your data than an OLS model without accounting for groups, I will edit this post with some R code for determining that.

Edit based on updated original post: Testing whether the partial pooling model is a better fit to the data than the complete pooling or no pooling models requires you to use the exactLRT() function in the RLRsim() package. Assuming you called the complete pooling model cpm the no pooling model npm, and the partial pooling model ppm, the code is as such:

# need to use full maximum likelihood on your partial pooling model for exactLRT() to work - I had suggested using restricted ML b/c of sample size
ppmA <- refitML(ppm) 

#likelihood ratio test comparing models
RLRsim::exactLRT(m=ppmA, m0=cpm) #test of partial vs complete pooling models
RLRsim::exactLRT(m=ppmA, m0=npm) #test of partial vs no pooling models

A significant test would indicate that the less parsimonious model is a better fit for the data. In your case, the ppm has more parameters than the cpm while the npm has more parameters than the ppm. I would expect the ppm to win out here, but will be interested to know what you find.

The anova() command works for comparing nested lmer() models. Note that when using REML, you can only use likelihood ratio testing to compare models with different random effects specifications.

  • $\begingroup$ Erik, I edited my original post to do what you mentioned (I think). I also added another question at the bottom. Thank you for taking the time to help me $\endgroup$ Dec 17, 2019 at 22:01
  • $\begingroup$ One thing to note is that the lmer() code you posted allows both each group's intercept and the slope of X on Y to vary across industries. I would suggest starting with just allowing intercepts to vary, and testing that model relative to the cpm and npm models. Then I would use anova to test whether the multilevel model with varying slopes and intercepts fits your data better than the model with just varying intercepts. $\endgroup$
    – Erik Ruzek
    Dec 17, 2019 at 22:46
  • $\begingroup$ Thanks Eric. Just what I was looking for. I'll post results tn $\endgroup$ Dec 17, 2019 at 22:51

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.