**Drunks&Lampposts » R**, and kindly contributed to R-bloggers)

A useful property of mixed effects and Bayesian hierarchical models is that lower level estimates are shrunk towards the more stable estimates further up the hierarchy.

To use a time honoured example you might be modelling the effect of a new teaching method on performance at the classroom level. Classes of 30 or so students are probably too small a sample to get useful results. In a hierarchical model the data are pooled so that all classes in a school are modelled together as a hierarchy and even all schools in a district.

At each level in this hierarchy an estimate for the efficiency of the teaching method is obtained. You will get an estimate for the school as a whole and for the district. You will even get estimates for the individual classes. These estimates will be weighted averages of the estimates for the class and the estimate for the school (which in turn is a weighted average of the estimate for the school and the district.) The clever part is that this weighting is itself determined by the data. Where a class is an outlier, and therefore the overall school average is less relevant, the estimate will be weighted towards the class. Where it is typical it will be weighted towards the school. This property is known as shrinkage.

I’m often interested in how much shrinkage is affecting my estimates and I want to *see *it. I’ve created this plot which I find useful. It’s done in R using ggplot and is very simple to code.

The idea is that the non shrunk estimates *b _{i}* (i.e. the estimates that would be obtained by modelling classes individually) are plotted on along the line x=y at the points (

*b*,

_{i}*b*). The estimates they are being shrunk towards

_{i}*a*are plotted at the points (

_{i}*b*,

_{i}*a*). Finally we plot the shrunk estimates

_{i}*s*at (

_{i}*b*,

_{i}*s*) and connect the points with an arrow to illustrate the direction of the shrinkage.

_{i}Here is an example. You can see the extent of the shrinkage by the the distance covered by the arrow towards the higher level estimate.

Note the arrows do sometimes point away from the higher level estimate. This is because this data is for a single coefficient in a hierarchical regression model with multiple coefficients. Where other coefficients have been stabilized by shrinkage this causes this particular coefficient to be revised.

The R code is as follows:

# *-------------------------------------------------------------------- # | FUNCTION: shrinkPlot # | Function for visualising shrinkage in hierarchical models # *-------------------------------------------------------------------- # | Version |Date |Programmer |Details of Change # | 01 |31/08/2013|Simon Raper |first version. # *-------------------------------------------------------------------- # | INPUTS: orig Estimates obtained from individual level # | modelling # | shrink Estimates obtained from hierarchical modelling # | prior Priors in Bayesian model or fixed effects in # | mixed effects model (i.e. what it is shrinking # | towards. # | window Limits for the plot (as a vector) # | # *-------------------------------------------------------------------- # | OUTPUTS: A ggplot object # *-------------------------------------------------------------------- # | DEPENDS: grid, ggplot2 # | # *-------------------------------------------------------------------- library(ggplot) library(grid) shrinkPlot<-function(orig, shrink, prior, window=NULL){ group<-factor(signif(prior,3)) data<-data.frame(orig, shrink, prior, group) g<-ggplot(data=data, aes(x=orig, xend=orig, y=orig, yend=shrink, col=group)) g2<-g+geom_segment(arrow = arrow(length = unit(0.3, "cm"))) +geom_point(data=comp.in, aes(x=coef, y=mean)) g3<-g2+xlab("Estimate")+ylab("Shrinkage")+ ggtitle("Shrinkage Plot") if (is.null(window)==FALSE){ g3<-g3+ylim(window)+xlim(window) } print(g3) }

**leave a comment**for the author, please follow the link and comment on their blog:

**Drunks&Lampposts » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...