How Can Grouping Improve Regression Performance 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. 

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)) 


Would you do the anova(lm(Y ~ X), lmer(Y ~ X + (X|Inudstries))) To see which model is better?
 A: 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:


*

*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().

*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. 
