Visualizing a multilevel model (HLM) in ggplot2 I have longitudinal data on several countries, looking at GDP and CO2 Emissions. In ggplot2, it is easy to make the software do something HLM-ish by plotting relationships separately for every country: 
 ggplot(dat, aes(x=CO2.Emissions, y=GDP, color=as.factor(Country))) + 
   geom_point(shape=20) + 
   geom_smooth(method=lm) + 
   theme(legend.position="none") + 
   scale_y_log10(name="Log10(GDP)") +
   scale_x_log10(name="Log10(CO2 Emissions)")

I get the following plot:

However, this is not a true plot of a multilevel model. I would love to do something LIKE this but visualizing results from a multilevel model. Specifically, the model is: 
 lmer(GDP ~ 1 + CO2.Emissions + (1 + CO2.Emissions | Country), data=dat )

This generates a random slope and intercept for each country. QUESTION: can I plot these and get something similar to (and as aesthetically pleasing as) the ggplot above? I want to visualize the relationships as depicted in the model, which ggplot2 is not doing.
Any help is appreciated!
 A: When you fit an lmer model, you can use the coef() function to extract the coefficients from the model. Your code would look something like:
mod1 <- lmer(GDP ~ 1 + CO2.Emissions + (1 + CO2.Emissions | Country), data=dat)

Then you can call coef() and extract the coefficients for each group by specifying:
coef(mod1)$Country

This will give you a vector of intercepts (the '1' you specified in the random term) and slopes (for 'C02.Emissions').
You can then save each of these into their own vector:
intercepts <- coef(mod1)$Country[,1] # Specifying the first column only
slopes <- coef(mod1)$Country[,2] # Specifying the second column only

Instead of calling geom_smooth(), you could then specify specific slopes and intercepts by adding this to your plot:
geom_abline(slope=slopes, intercept=intercepts)

The positive of this is that it is using the model-implied slopes and intercepts. The downside is that it will extrapolate the lines beyond the values for each cluster (in this case, 'Country').
I would then add another geom_abline that is the average slope and intercept, which you can get from:
summary(mod1)$coef

Another way of doing this—that is not using the model-implied slopes and intercepts—is by specifying group=cluster. Using this, it fits a different OLS line for each cluster (which is not what the multilevel model fit with lmer() is doing, obviously).
I've done this before, and it adapted to your variables would look something like:
ggplot(datalong, aes(x=CO2.Emissions, y=GDP, group=Country))+
  stat_smooth(method="lm", se=FALSE, size=.5, color="springgreen") + # slopes for different countries
  stat_smooth(aes(group=1), method="lm", color="blue", size=1.5) + # average slope with SE

This might be slightly easier to do, but it will not match up with the model-implied slopes and intercepts you get from the coef(lmer(...))$cluster approach.
A: When you want to plot lmer()objects, I find it easiest to use predict(). First you fit your model:  
random.coef.model <- lmer(GDP ~ 1 + CO2.Emissions + (1 + CO2.Emissions | Country), 
                          data=dat)

Then you predict GDP values corresponding to your predictor variable (CO2.Emissions):  
dat$random.coefficients.predictions <- predict(random.coef.model)

And then you are free to plot them either using geom_smooth(se=FALSE) or geom_line(). If you want to have a scatterplot at the same time, you need to feed geom_point() new aes(), as the existing y-values would be the predicted values.  So:  
random.coef.graph <- ggplot(aes(x = CO2.Emissions, y = random.coefficients.predictions, 
                                color = as.factor(Country)), data = dat) +
geom_line(size=.3) +
geom_point(aes(y = GDP)) +
ggthemes::theme_tufte() #just to make it nice!

