Recent adventures with lazyeval

[This article was first published on Higher Order Functions, 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.

The lazyeval package is a
tool-set for performing nonstandard evaluation in R.
Nonstandard evaluation refers to any situation where something special happens
with how user input or code is evaluated.

For example, the library function doesn’t evaluate variables. In the
example below, I try to trick library into loading a fake package called
evil_package by assigning to the package name lazyeval. In other words, we
have the expression lazyeval and its value is "evil_package".

<span class="n">print</span><span class="p">(</span><span class="n">.packages</span><span class="p">())</span><span class="w">
</span><span class="c1">#> [1] "stringr"   "dplyr"     "knitr"     "stats"     "graphics"  "grDevices"
#> [7] "utils"     "datasets"  "base"
</span><span class="w">
</span><span class="n">lazyeval</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="s2">"evil_package"</span><span class="w">
</span><span class="n">library</span><span class="p">(</span><span class="n">lazyeval</span><span class="p">)</span><span class="w">

</span><span class="c1"># The lazyeval package is loaded now.
</span><span class="n">print</span><span class="p">(</span><span class="n">.packages</span><span class="p">())</span><span class="w">
</span><span class="c1">#>  [1] "lazyeval"  "stringr"   "dplyr"     "knitr"     "stats"    
#>  [6] "graphics"  "grDevices" "utils"     "datasets"  "base"
</span>

But this gambit doesn’t work because library did something special: It didn’t
evaluate the expression lazyeval. In the source code of library, there is a
line package <- as.character(substitute(package)) which replaces the value of
package with a character version of the expression that the user wrote.

That’s a simple example of nonstandard evaluation, but it’s pervasive. It’s why
you never have to quote column names in dplyr or ggplot2. In this post, I
present some recent examples where I decided to use the lazyeval package to do
something nonstandard. These examples are straight out of the
lazyeval vignette
in terms of complexity, but that’s fine. We all have to start somewhere.

Capturing expressions

Plot titles. While reading the book Machine Learning For Hackers, I wanted
to plot random numbers generated by probability distributions discussed by the
book. I used the lazyeval::expr_text function to capture the command used to
generate the numbers and write it as the title of the plot.

<span class="n">library</span><span class="p">(</span><span class="n">dplyr</span><span class="p">,</span><span class="w"> </span><span class="n">warn.conflicts</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">FALSE</span><span class="p">)</span><span class="w">
</span><span class="n">library</span><span class="p">(</span><span class="n">ggplot2</span><span class="p">)</span><span class="w">

</span><span class="n">plot_dist</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="k">function</span><span class="p">(</span><span class="n">xs</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
  </span><span class="n">data</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">data_frame</span><span class="p">(</span><span class="n">x</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">xs</span><span class="p">)</span><span class="w">
  </span><span class="n">ggplot</span><span class="p">(</span><span class="n">data</span><span class="p">)</span><span class="w"> </span><span class="o">+</span><span class="w">
    </span><span class="n">aes</span><span class="p">(</span><span class="n">x</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">x</span><span class="p">)</span><span class="w"> </span><span class="o">+</span><span class="w">
    </span><span class="n">geom_density</span><span class="p">()</span><span class="w"> </span><span class="o">+</span><span class="w">
    </span><span class="n">ggtitle</span><span class="p">(</span><span class="n">lazyeval</span><span class="o">::</span><span class="n">expr_text</span><span class="p">(</span><span class="n">xs</span><span class="p">))</span><span class="w"> 
</span><span class="p">}</span><span class="w">

</span><span class="n">plot_dist</span><span class="p">(</span><span class="n">rcauchy</span><span class="p">(</span><span class="n">n</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">250</span><span class="p">,</span><span class="w"> </span><span class="n">location</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0</span><span class="p">,</span><span class="w"> </span><span class="n">scale</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">1</span><span class="p">))</span><span class="w">
</span>

center

<span class="n">plot_dist</span><span class="p">(</span><span class="n">rgamma</span><span class="p">(</span><span class="n">n</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">25000</span><span class="p">,</span><span class="w"> </span><span class="n">shape</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">5</span><span class="p">,</span><span class="w"> </span><span class="n">rate</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">.5</span><span class="p">))</span><span class="w">
</span>

center

<span class="n">plot_dist</span><span class="p">(</span><span class="n">rexp</span><span class="p">(</span><span class="n">n</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">25000</span><span class="p">,</span><span class="w"> </span><span class="n">rate</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">.5</span><span class="p">))</span><span class="w">
</span>

center

Less fussy warning messages. I recently inherited some code where there were
custom warning messages based on the input. The code threw a warning whenever a
duplicate participant ID was found in a survey. It went something like this:

<span class="c1"># some dummy data
</span><span class="n">study1</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">data_frame</span><span class="p">(</span><span class="w">
  </span><span class="n">id</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nf">c</span><span class="p">(</span><span class="m">1</span><span class="p">,</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="m">3</span><span class="p">,</span><span class="w"> </span><span class="m">4</span><span class="p">),</span><span class="w"> 
  </span><span class="n">response</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nf">c</span><span class="p">(</span><span class="s2">"b"</span><span class="p">,</span><span class="w"> </span><span class="s2">"c"</span><span class="p">,</span><span class="w"> </span><span class="s2">"a"</span><span class="p">,</span><span class="w"> </span><span class="s2">"b"</span><span class="p">))</span><span class="w">

</span><span class="n">study2</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">data_frame</span><span class="p">(</span><span class="w">
  </span><span class="n">id</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nf">c</span><span class="p">(</span><span class="m">1</span><span class="p">,</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="m">3</span><span class="p">,</span><span class="w"> </span><span class="m">1</span><span class="p">),</span><span class="w"> 
  </span><span class="n">response</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nf">c</span><span class="p">(</span><span class="s2">"a"</span><span class="p">,</span><span class="w"> </span><span class="s2">"a"</span><span class="p">,</span><span class="w"> </span><span class="s2">"a"</span><span class="p">,</span><span class="w"> </span><span class="s2">"b"</span><span class="p">,</span><span class="w"> </span><span class="s2">"c"</span><span class="p">))</span><span class="w">

</span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">anyDuplicated</span><span class="p">(</span><span class="n">study1</span><span class="o">$</span><span class="n">id</span><span class="p">))</span><span class="w"> </span><span class="p">{</span><span class="w">
  </span><span class="n">warning</span><span class="p">(</span><span class="s2">"Duplicate IDs found in Study1"</span><span class="p">,</span><span class="w"> </span><span class="n">call.</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">FALSE</span><span class="p">)</span><span class="w">
</span><span class="p">}</span><span class="w">

</span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">anyDuplicated</span><span class="p">(</span><span class="n">study2</span><span class="o">$</span><span class="n">id</span><span class="p">))</span><span class="w"> </span><span class="p">{</span><span class="w">
  </span><span class="n">warning</span><span class="p">(</span><span class="s2">"Duplicate IDs found in Study2"</span><span class="p">,</span><span class="w"> </span><span class="n">call.</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">FALSE</span><span class="p">)</span><span class="w">
</span><span class="p">}</span><span class="w">
</span><span class="c1">#> Warning: Duplicate IDs found in Study2
</span>

To extend this code to a new study, one would just copy-and-paste and update the
if statement’s condition and warning messages. Like so:

<span class="n">study3</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">data_frame</span><span class="p">(</span><span class="w">
  </span><span class="n">id</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nf">c</span><span class="p">(</span><span class="m">1</span><span class="p">,</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="m">3</span><span class="p">,</span><span class="w"> </span><span class="m">2</span><span class="p">),</span><span class="w"> 
  </span><span class="n">response</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nf">c</span><span class="p">(</span><span class="s2">"b"</span><span class="p">,</span><span class="w"> </span><span class="s2">"c"</span><span class="p">,</span><span class="w"> </span><span class="s2">"a"</span><span class="p">,</span><span class="w"> </span><span class="s2">"b"</span><span class="p">))</span><span class="w">

</span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">anyDuplicated</span><span class="p">(</span><span class="n">study3</span><span class="o">$</span><span class="n">id</span><span class="p">))</span><span class="w"> </span><span class="p">{</span><span class="w">
  </span><span class="n">warning</span><span class="p">(</span><span class="s2">"Duplicate IDs found in Study2"</span><span class="p">,</span><span class="w"> </span><span class="n">call.</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">FALSE</span><span class="p">)</span><span class="w">
</span><span class="p">}</span><span class="w">
</span><span class="c1">#> Warning: Duplicate IDs found in Study2
</span>

Wait, that’s not right! I forgot to update the warning message…

This setup is too brittle for me, so I abstracted the procedure into a
function. First, I wrote a helper function to print out duplicates elements in a
vector. This helper will let us make the warning message a little more
informative.

<span class="c1"># Print out duplicated elements in a vector
</span><span class="n">print_duplicates</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="k">function</span><span class="p">(</span><span class="n">xs</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
  </span><span class="n">duplicated</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">xs</span><span class="p">[</span><span class="n">duplicated</span><span class="p">(</span><span class="n">xs</span><span class="p">)]</span><span class="w">
  </span><span class="n">duplicated</span><span class="w"> </span><span class="o">%>%</span><span class="w"> </span><span class="n">sort</span><span class="w"> </span><span class="o">%>%</span><span class="w"> </span><span class="n">unique</span><span class="w"> </span><span class="o">%>%</span><span class="w"> </span><span class="n">paste0</span><span class="p">(</span><span class="n">collapse</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">", "</span><span class="p">)</span><span class="w">
</span><span class="p">}</span><span class="w">

</span><span class="n">print_duplicates</span><span class="p">(</span><span class="n">study2</span><span class="o">$</span><span class="n">id</span><span class="p">)</span><span class="w">
</span><span class="c1">#> [1] "1, 2"
</span>

Next, I wrote a function to issue the warnings. I used lazyeval::expr_label
convert the user-inputted expression into a string wrapped in backticks.

<span class="c1"># Print a warning if duplicates are found in a vector
</span><span class="n">warn_duplicates</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="k">function</span><span class="p">(</span><span class="n">xs</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
  </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">anyDuplicated</span><span class="p">(</span><span class="n">xs</span><span class="p">))</span><span class="w"> </span><span class="p">{</span><span class="w">
    </span><span class="c1"># Get what the user wrote for the xs argument
</span><span class="w">    </span><span class="n">actual_xs</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">lazyeval</span><span class="o">::</span><span class="n">expr_label</span><span class="p">(</span><span class="n">xs</span><span class="p">)</span><span class="w">
    </span><span class="n">msg</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">paste0</span><span class="p">(</span><span class="s2">"Duplicate entries in "</span><span class="p">,</span><span class="w"> </span><span class="n">actual_xs</span><span class="p">,</span><span class="w"> </span><span class="s2">": "</span><span class="p">,</span><span class="w">
                  </span><span class="n">print_duplicates</span><span class="p">(</span><span class="n">xs</span><span class="p">))</span><span class="w">
    </span><span class="n">warning</span><span class="p">(</span><span class="n">msg</span><span class="p">,</span><span class="w"> </span><span class="n">call.</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">FALSE</span><span class="p">)</span><span class="w">
  </span><span class="p">}</span><span class="w">
  </span><span class="nf">invisible</span><span class="p">(</span><span class="kc">NULL</span><span class="p">)</span><span class="w">
</span><span class="p">}</span><span class="w">

</span><span class="n">warn_duplicates</span><span class="p">(</span><span class="n">study1</span><span class="o">$</span><span class="n">id</span><span class="p">)</span><span class="w">
</span><span class="n">warn_duplicates</span><span class="p">(</span><span class="n">study2</span><span class="o">$</span><span class="n">id</span><span class="p">)</span><span class="w">
</span><span class="c1">#> Warning: Duplicate entries in `study2$id`: 1, 2
</span><span class="n">warn_duplicates</span><span class="p">(</span><span class="n">study3</span><span class="o">$</span><span class="n">id</span><span class="p">)</span><span class="w">
</span><span class="c1">#> Warning: Duplicate entries in `study3$id`: 2
</span>

The advantage of this approach is that the warning is a generic message that can
work on any input. But in a funny way, the warning is also customized for the
input because it includes the input printed verbatim.

An aside: In plotting, I used lazyeval::expr_text, but here I used
lazyeval::expr_label. The two differ slighty. Namely, expr_label surrounds
the captured expression with backticks to indicate that expression is code:

<span class="n">lazyeval</span><span class="o">::</span><span class="n">expr_text</span><span class="p">(</span><span class="n">stop</span><span class="p">())</span><span class="w">
</span><span class="c1">#> [1] "stop()"
</span><span class="n">lazyeval</span><span class="o">::</span><span class="n">expr_label</span><span class="p">(</span><span class="n">stop</span><span class="p">())</span><span class="w">
</span><span class="c1">#> [1] "`stop()`"
</span>

Asking questions about a posterior distribution

I fit regression models with RStanARM. It lets me fit Bayesian models in
Stan by writing conventional R modeling code. (Btw, I’m
giving a tutorial on RStanARM in a month.)

Here’s a model about some famous flowers.

<span class="n">library</span><span class="p">(</span><span class="n">rstanarm</span><span class="p">)</span><span class="w">

</span><span class="n">model</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">stan_glm</span><span class="p">(</span><span class="w">
  </span><span class="n">Petal.Width</span><span class="w"> </span><span class="o">~</span><span class="w"> </span><span class="n">Petal.Length</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">Species</span><span class="p">,</span><span class="w">
  </span><span class="n">data</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">iris</span><span class="p">,</span><span class="w">
  </span><span class="n">family</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">gaussian</span><span class="p">(),</span><span class="w"> 
  </span><span class="n">prior</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">normal</span><span class="p">(</span><span class="m">0</span><span class="p">,</span><span class="w"> </span><span class="m">1</span><span class="p">))</span><span class="w">
</span>

Here’s the essential relationship being explored.

<span class="n">ggplot</span><span class="p">(</span><span class="n">iris</span><span class="p">)</span><span class="w"> </span><span class="o">+</span><span class="w"> 
  </span><span class="n">aes</span><span class="p">(</span><span class="n">x</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">Petal.Length</span><span class="p">,</span><span class="w"> </span><span class="n">y</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">Petal.Width</span><span class="p">,</span><span class="w"> </span><span class="n">color</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">Species</span><span class="p">)</span><span class="w"> </span><span class="o">+</span><span class="w"> 
  </span><span class="n">geom_point</span><span class="p">()</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">stat_smooth</span><span class="p">(</span><span class="n">method</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"lm"</span><span class="p">)</span><span class="w">
</span>

center

The model gives me 4000 samples from the posterior distribution of the model.

<span class="n">summary</span><span class="p">(</span><span class="n">model</span><span class="p">)</span><span class="w">
</span><span class="c1">#> stan_glm(formula = Petal.Width ~ Petal.Length * Species, family = gaussian(), 
#>     data = iris, prior = normal(0, 1))
#> 
#> Family: gaussian (identity)
#> Algorithm: sampling
#> Posterior sample size: 4000
#> Observations: 150
#> 
#> Estimates:
#>                                  mean   sd   2.5%   25%   50%   75%
#> (Intercept)                     0.0    0.2 -0.4   -0.1   0.0   0.1 
#> Petal.Length                    0.2    0.1  0.0    0.1   0.2   0.3 
#> Speciesversicolor              -0.1    0.3 -0.6   -0.3  -0.1   0.1 
#> Speciesvirginica                1.1    0.3  0.5    0.9   1.1   1.3 
#> Petal.Length:Speciesversicolor  0.1    0.1 -0.1    0.1   0.1   0.2 
#> Petal.Length:Speciesvirginica   0.0    0.1 -0.2   -0.1   0.0   0.1 
#> sigma                           0.2    0.0  0.2    0.2   0.2   0.2 
#> mean_PPD                        1.2    0.0  1.2    1.2   1.2   1.2 
#> log-posterior                  32.9    1.9 28.2   31.8  33.2  34.3 
#>                                  97.5%
#> (Intercept)                     0.3   
#> Petal.Length                    0.4   
#> Speciesversicolor               0.5   
#> Speciesvirginica                1.7   
#> Petal.Length:Speciesversicolor  0.4   
#> Petal.Length:Speciesvirginica   0.2   
#> sigma                           0.2   
#> mean_PPD                        1.2   
#> log-posterior                  35.7   
#> 
#> Diagnostics:
#>                                mcse Rhat n_eff
#> (Intercept)                    0.0  1.0   704 
#> Petal.Length                   0.0  1.0   711 
#> Speciesversicolor              0.0  1.0   983 
#> Speciesvirginica               0.0  1.0   960 
#> Petal.Length:Speciesversicolor 0.0  1.0   699 
#> Petal.Length:Speciesvirginica  0.0  1.0   638 
#> sigma                          0.0  1.0  2593 
#> mean_PPD                       0.0  1.0  3365 
#> log-posterior                  0.1  1.0  1204 
#> 
#> For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
</span>

At the 2.5% quantile, the Petal.Length effect looks like zero or less than
zero. What proportion of the Petal.Length effects (for setosa flowers) is
positive?

To answer questions like this one in a convenient way, I wrote a function that
takes a boolean expression about a model’s parameters and evaluates it inside of
the data-frame summary of the model posterior distribution. lazyeval::f_eval
does the nonstandard evaluation: It evaluates an expression captured by a
formula like ~ 0 < Petal.Length inside of a list or data-frame. (Note that the
mean of a logical vector is the proportion of the elements that are true.)

<span class="c1"># Get proportion of posterior samples satisfying some inequality
</span><span class="n">posterior_proportion_</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="k">function</span><span class="p">(</span><span class="n">model</span><span class="p">,</span><span class="w"> </span><span class="n">inequality</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
  </span><span class="n">draws</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">as.data.frame</span><span class="p">(</span><span class="n">model</span><span class="p">)</span><span class="w">
  </span><span class="n">mean</span><span class="p">(</span><span class="n">lazyeval</span><span class="o">::</span><span class="n">f_eval</span><span class="p">(</span><span class="n">inequality</span><span class="p">,</span><span class="w"> </span><span class="n">data</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">draws</span><span class="p">))</span><span class="w">
</span><span class="p">}</span><span class="w">

</span><span class="n">posterior_proportion_</span><span class="p">(</span><span class="n">model</span><span class="p">,</span><span class="w"> </span><span class="o">~</span><span class="w"> </span><span class="m">0</span><span class="w"> </span><span class="o"><</span><span class="w"> </span><span class="n">Petal.Length</span><span class="p">)</span><span class="w">
</span><span class="c1">#> [1] 0.9325
</span>

But all those tildes… The final underscore in posterior_proportion_
follows a convention for distinguishing between nonstandard evaluation functions
that require formulas and those that do not. In the dplyr package, for example,
there is select_/select/, mutate_/mutate, and so on. We can do the
formula-free form of this function by using lazyeval::f_capture to capture the
input expression as a formula. We then hand that off to the version of the
function that takes a formula.

<span class="n">posterior_proportion</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="k">function</span><span class="p">(</span><span class="n">model</span><span class="p">,</span><span class="w"> </span><span class="n">expr</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
  </span><span class="n">posterior_proportion_</span><span class="p">(</span><span class="n">model</span><span class="p">,</span><span class="w"> </span><span class="n">lazyeval</span><span class="o">::</span><span class="n">f_capture</span><span class="p">(</span><span class="n">expr</span><span class="p">))</span><span class="w">
</span><span class="p">}</span><span class="w">

</span><span class="n">posterior_proportion</span><span class="p">(</span><span class="n">model</span><span class="p">,</span><span class="w"> </span><span class="m">0</span><span class="w"> </span><span class="o"><</span><span class="w"> </span><span class="n">Petal.Length</span><span class="p">)</span><span class="w">
</span><span class="c1">#> [1] 0.9325
</span>

Here’s another question: What proportion of the posterior of the Petal.Length
effect for virginica flowers is positive? In classical models, we would change
the reference level for the categorical variable, refit the model, and check the
significance. But I don’t want to refit this model because that would repeat the
MCMC sampling. (It takes about 30 seconds to fit this model after all!) I’ll
just ask the model for the sum of Petal.Length and
Petal.Length:Speciesversicolor effects. That will give me the estimated
Petal.Length effect but adjusted for the versicolor species.

<span class="n">posterior_proportion</span><span class="p">(</span><span class="n">model</span><span class="p">,</span><span class="w"> </span><span class="m">0</span><span class="w"> </span><span class="o"><</span><span class="w"> </span><span class="n">Petal.Length</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">`Petal.Length:Speciesversicolor`</span><span class="p">)</span><span class="w">
</span><span class="c1">#> [1] 1
</span><span class="w">
</span><span class="n">posterior_proportion</span><span class="p">(</span><span class="n">model</span><span class="p">,</span><span class="w"> </span><span class="m">0</span><span class="w"> </span><span class="o"><</span><span class="w"> </span><span class="n">Petal.Length</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">`Petal.Length:Speciesvirginica`</span><span class="p">)</span><span class="w">
</span><span class="c1">#> [1] 0.99975
</span>

(The backticks around Petal.Length:Speciesversicolor here prevent the :
symbol from being evaluated as an operator.)

To leave a comment for the author, please follow the link and comment on their blog: Higher Order Functions.

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)