Visualizing bivariate shrinkage

[This article was first published on Ecology in silico, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Inspired by this post about visualizing shrinkage on Coppelia, and this thread about visualizing mixed models on Stack Exchange, I started thinking about how to visualize shrinkage in more than one dimension. One might find themselves in this situation with a varying slope, varying intercept hierarichical (mixed effects) model, a model with two varying intercepts, etc. Then I found a great bivariate shrinkage plot in Doug Bates’ lme4 book, Figure 3.8.

Here is a snippet to reproduce a similar bivariate shrinkage plot in ggplot2, adding a color coded probability density surface and contours for the estimated multivariate normal distribution of random effects, using the same sleep study data that Bates used.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
<span class="line"><span class="c1">## 2d shrikage for random slope, random intercept model</span>
</span><span class="line">library<span class="p">(</span>lme4<span class="p">)</span>
</span><span class="line">library<span class="p">(</span>ggplot2<span class="p">)</span>
</span><span class="line">library<span class="p">(</span>grid<span class="p">)</span>
</span><span class="line">library<span class="p">(</span>mvtnorm<span class="p">)</span>
</span><span class="line">library<span class="p">(</span>reshape2<span class="p">)</span>
</span><span class="line">library<span class="p">(</span>ellipse<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># fit models</span>
</span><span class="line">m0 <span class="o"><-</span> lm<span class="p">(</span>Reaction <span class="o">~</span> Days <span class="o">*</span> Subject<span class="p">,</span> data<span class="o">=</span>sleepstudy<span class="p">)</span>
</span><span class="line">m1 <span class="o"><-</span> lmer<span class="p">(</span>Reaction <span class="o">~</span> Days <span class="o">+</span> <span class="p">(</span>Days <span class="o">|</span> Subject<span class="p">),</span> sleepstudy<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># extract fixed effect estimates</span>
</span><span class="line">intercepts <span class="o"><-</span> grepl<span class="p">(</span><span class="s">"Days"</span><span class="p">,</span> names<span class="p">(</span>coef<span class="p">(</span>m0<span class="p">)))</span> <span class="o">==</span> <span class="kc">FALSE</span>
</span><span class="line">m0_intercepts <span class="o"><-</span> coef<span class="p">(</span>m0<span class="p">)[</span>intercepts<span class="p">]</span>
</span><span class="line">m0_intercepts<span class="p">[</span><span class="m">-1</span><span class="p">]</span> <span class="o"><-</span> m0_intercepts<span class="p">[</span><span class="m">-1</span><span class="p">]</span> <span class="o">+</span> m0_intercepts<span class="p">[</span><span class="m">1</span><span class="p">]</span>
</span><span class="line">slopes <span class="o"><-</span> <span class="o">!</span>intercepts
</span><span class="line">m0_slopes <span class="o"><-</span> coef<span class="p">(</span>m0<span class="p">)[</span>slopes<span class="p">]</span>
</span><span class="line">m0_slopes<span class="p">[</span><span class="m">-1</span><span class="p">]</span> <span class="o"><-</span> m0_slopes<span class="p">[</span><span class="m">-1</span><span class="p">]</span> <span class="o">+</span> m0_slopes<span class="p">[</span><span class="m">1</span><span class="p">]</span>
</span><span class="line">
</span><span class="line"><span class="c1"># extract random effect estimates</span>
</span><span class="line">m1_intercepts <span class="o"><-</span> coef<span class="p">(</span>m1<span class="p">)</span><span class="o">$</span>Subject<span class="p">[,</span> <span class="m">1</span><span class="p">]</span>
</span><span class="line">m1_slopes <span class="o"><-</span> coef<span class="p">(</span>m1<span class="p">)</span><span class="o">$</span>Subject<span class="p">[,</span> <span class="m">2</span><span class="p">]</span>
</span><span class="line">
</span><span class="line">d <span class="o"><-</span> data.frame<span class="p">(</span>interceptF <span class="o">=</span> m0_intercepts<span class="p">,</span>
</span><span class="line">                slopeF <span class="o">=</span> m0_slopes<span class="p">,</span>
</span><span class="line">                interceptR <span class="o">=</span> m1_intercepts<span class="p">,</span>
</span><span class="line">                slopeR <span class="o">=</span> m1_slopes
</span><span class="line">                <span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># 95% bivariate normal ellipse for random effects</span>
</span><span class="line">df_ell <span class="o"><-</span> data.frame<span class="p">(</span>ellipse<span class="p">(</span>VarCorr<span class="p">(</span>m1<span class="p">)</span><span class="o">$</span>Subject<span class="p">,</span> centre<span class="o">=</span>fixef<span class="p">(</span>m1<span class="p">)))</span>
</span><span class="line">names<span class="p">(</span>df_ell<span class="p">)</span> <span class="o"><-</span> c<span class="p">(</span><span class="s">"intercept"</span><span class="p">,</span> <span class="s">"slope"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># bivariate normal density surface</span>
</span><span class="line">lo <span class="o"><-</span> <span class="m">200</span> <span class="c1"># length.out of x and y grid</span>
</span><span class="line">xvals <span class="o"><-</span> seq<span class="p">(</span>from <span class="o">=</span> min<span class="p">(</span>df_ell<span class="o">$</span>intercept<span class="p">)</span> <span class="o">-</span> <span class="m">.1</span> <span class="o">*</span> abs<span class="p">(</span>min<span class="p">(</span>df_ell<span class="o">$</span>intercept<span class="p">)),</span>
</span><span class="line">             to <span class="o">=</span> max<span class="p">(</span>df_ell<span class="o">$</span>intercept<span class="p">)</span> <span class="o">+</span> <span class="m">.05</span> <span class="o">*</span> abs<span class="p">(</span>max<span class="p">(</span>df_ell<span class="o">$</span>intercept<span class="p">)),</span>
</span><span class="line">             length.out <span class="o">=</span> lo<span class="p">)</span>
</span><span class="line">yvals <span class="o"><-</span> seq<span class="p">(</span>from <span class="o">=</span> min<span class="p">(</span>df_ell<span class="o">$</span>slope<span class="p">)</span> <span class="o">-</span> <span class="m">.4</span> <span class="o">*</span> abs<span class="p">(</span>min<span class="p">(</span>df_ell<span class="o">$</span>slope<span class="p">)),</span>
</span><span class="line">             to <span class="o">=</span> max<span class="p">(</span>df_ell<span class="o">$</span>slope<span class="p">)</span> <span class="o">+</span> <span class="m">.1</span> <span class="o">*</span> abs<span class="p">(</span>max<span class="p">(</span>df_ell<span class="o">$</span>slope<span class="p">)),</span>
</span><span class="line">             length.out <span class="o">=</span> lo<span class="p">)</span>
</span><span class="line">
</span><span class="line">z <span class="o"><-</span> matrix<span class="p">(</span><span class="m">0</span><span class="p">,</span> lo<span class="p">,</span> lo<span class="p">)</span>
</span><span class="line"><span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>lo<span class="p">)</span> <span class="p">{</span>
</span><span class="line">  x <span class="o">=</span> xvals
</span><span class="line">  y <span class="o">=</span> yvals<span class="p">[</span>i<span class="p">]</span>
</span><span class="line">  z<span class="p">[,</span>i<span class="p">]</span> <span class="o">=</span> dmvnorm<span class="p">(</span>cbind<span class="p">(</span>x<span class="p">,</span>y<span class="p">),</span>
</span><span class="line">                  mean <span class="o">=</span> fixef<span class="p">(</span>m1<span class="p">),</span>
</span><span class="line">                  sigma <span class="o">=</span> VarCorr<span class="p">(</span>m1<span class="p">)</span><span class="o">$</span>Subject<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line">mv_ranef <span class="o"><-</span> melt<span class="p">(</span>z<span class="p">)</span>
</span><span class="line">names<span class="p">(</span>mv_ranef<span class="p">)</span> <span class="o"><-</span> c<span class="p">(</span><span class="s">"x"</span><span class="p">,</span> <span class="s">"y"</span><span class="p">,</span> <span class="s">"z"</span><span class="p">)</span>
</span><span class="line">mv_ranef<span class="o">$</span>x <span class="o"><-</span> xvals<span class="p">[</span>mv_ranef<span class="o">$</span>x<span class="p">]</span>
</span><span class="line">mv_ranef<span class="o">$</span>y <span class="o"><-</span> yvals<span class="p">[</span>mv_ranef<span class="o">$</span>y<span class="p">]</span>
</span><span class="line">
</span><span class="line">
</span><span class="line">p <span class="o"><-</span> ggplot<span class="p">(</span>d<span class="p">)</span> <span class="o">+</span>
</span><span class="line">  geom_raster<span class="p">(</span>aes<span class="p">(</span>x<span class="o">=</span>x<span class="p">,</span> y<span class="o">=</span>y<span class="p">,</span> fill<span class="o">=</span>z<span class="p">),</span> data<span class="o">=</span>mv_ranef<span class="p">)</span> <span class="o">+</span>
</span><span class="line">  scale_fill_gradient<span class="p">(</span>low<span class="o">=</span><span class="s">"white"</span><span class="p">,</span> high<span class="o">=</span><span class="s">"red"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  guides<span class="p">(</span>fill<span class="o">=</span><span class="kc">FALSE</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  geom_path<span class="p">(</span>data<span class="o">=</span>df_ell<span class="p">,</span> aes<span class="p">(</span>x<span class="o">=</span>intercept<span class="p">,</span> y<span class="o">=</span>slope<span class="p">),</span> size<span class="o">=</span><span class="m">.5</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  geom_contour<span class="p">(</span>aes<span class="p">(</span>x<span class="o">=</span>x<span class="p">,</span> y<span class="o">=</span>y<span class="p">,</span> z<span class="o">=</span>z<span class="p">),</span> data<span class="o">=</span>mv_ranef<span class="p">,</span> size<span class="o">=</span><span class="m">.1</span><span class="p">,</span> color<span class="o">=</span><span class="s">"black"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  geom_segment<span class="p">(</span>aes<span class="p">(</span>x<span class="o">=</span>interceptF<span class="p">,</span> y<span class="o">=</span>slopeF<span class="p">,</span>
</span><span class="line">                   xend<span class="o">=</span>interceptR<span class="p">,</span> yend<span class="o">=</span>slopeR<span class="p">),</span>
</span><span class="line">               arrow <span class="o">=</span> arrow<span class="p">(</span>length <span class="o">=</span> unit<span class="p">(</span><span class="m">0.3</span><span class="p">,</span> <span class="s">"cm"</span><span class="p">)),</span>
</span><span class="line">               alpha<span class="o">=</span><span class="m">.7</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  xlab<span class="p">(</span><span class="s">"Estimated intercepts"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  ylab<span class="p">(</span><span class="s">"Estimated slopes"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  theme<span class="p">(</span>legend.position<span class="o">=</span><span class="s">"none"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  ggtitle<span class="p">(</span><span class="s">"Bivariate shrinkage plot"</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">  theme_bw<span class="p">()</span>  <span class="o">+</span>
</span><span class="line">  theme<span class="p">(</span>panel.grid.major <span class="o">=</span> element_blank<span class="p">(),</span> panel.grid.minor <span class="o">=</span> element_blank<span class="p">())</span>
</span><span class="line">p
</span>

To leave a comment for the author, please follow the link and comment on their blog: Ecology in silico.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)