Conditional ggplot2 geoms in functions (QTL plots)

[This article was first published on Shirin's playgRound, 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.

When running an analysis, I am usually combining functions from multiple packages. Most of these packages come with their own plotting functions. And while they are certainly convenient in that they allow me to get a quick glance at the data or the output, they all have their own style. If I want to prepare a report, proposal or a paper though, I want all my plots to come from a single cast so that they give a consistent feel to the story I want to tell with my data.

I have grown very fond of ggplot2 for producing plots, because of it’s tidy grammar, flexibility and easy customization. Here, I want to show an example of my own ggplot2 function to produce QTL plots for Karl Broman’s qtl package.

Specifically, I want to show how to incorporate conditional geoms when using ggplot2 in a function call.

I will show how my function can be used with the output from qtl’s scanone() function by following the examples from the package tutorial. The tutorial does an excellent job describing the workflow, so I am not going to explain it here.

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

For defining the function and preparing the data for plotting, I am working with ggplot2, ggrepel (for adding text labels that don’t overlap), tidyr and dplyr. This is my function code:

<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">library</span><span class="p">(</span><span class="n">ggrepel</span><span class="p">)</span><span class="w">
</span><span class="n">library</span><span class="p">(</span><span class="n">tidyr</span><span class="p">)</span><span class="w">
</span><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">qtl_plot</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">input</span><span class="p">,</span><span class="w">              </span><span class="c1"># data frame input from scanone
</span><span class="w">                     </span><span class="n">mult.pheno</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="c1"># multiple phenotypes?
</span><span class="w">                     </span><span class="n">model</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"normal"</span><span class="p">,</span><span class="w">   </span><span class="c1"># model used in scanone
</span><span class="w">                     </span><span class="n">chrs</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">NA</span><span class="p">,</span><span class="w">          </span><span class="c1"># chromosomes to display
</span><span class="w">                     </span><span class="n">lod</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">NA</span><span class="p">,</span><span class="w">           </span><span class="c1"># LOD threshold
</span><span class="w">                     </span><span class="n">rug</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="c1"># plot marker positions as rug?
</span><span class="w">                     </span><span class="n">ncol</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">NA</span><span class="p">,</span><span class="w">          </span><span class="c1"># number of columns for facetting
</span><span class="w">                     </span><span class="n">labels</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">NA</span><span class="w">         </span><span class="c1"># optional dataframe to plot QTL labels
</span><span class="w">                     </span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
                     
                  </span><span class="c1"># if we have multiple phenotypes and/or a 2part model, gather input
</span><span class="w">                  </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">mult.pheno</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="n">model</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="s2">"2part"</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
                    </span><span class="n">input</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">gather</span><span class="p">(</span><span class="n">input</span><span class="p">,</span><span class="w"> </span><span class="n">group</span><span class="p">,</span><span class="w"> </span><span class="n">lod</span><span class="p">,</span><span class="w"> </span><span class="n">grep</span><span class="p">(</span><span class="s2">"pheno"</span><span class="p">,</span><span class="w"> </span><span class="n">colnames</span><span class="p">(</span><span class="n">input</span><span class="p">)))</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="k">else</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">mult.pheno</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
                    </span><span class="n">input</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">gather</span><span class="p">(</span><span class="n">input</span><span class="p">,</span><span class="w"> </span><span class="n">group</span><span class="p">,</span><span class="w"> </span><span class="n">lod</span><span class="p">,</span><span class="w"> </span><span class="n">grep</span><span class="p">(</span><span class="s2">"pheno"</span><span class="p">,</span><span class="w"> </span><span class="n">colnames</span><span class="p">(</span><span class="n">input</span><span class="p">)))</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="k">else</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">model</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="s2">"2part"</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
                    </span><span class="n">input</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">gather</span><span class="p">(</span><span class="n">input</span><span class="p">,</span><span class="w"> </span><span class="n">method</span><span class="p">,</span><span class="w"> </span><span class="n">lod</span><span class="p">,</span><span class="w"> </span><span class="n">lod.p.mu</span><span class="o">:</span><span class="n">lod.mu</span><span class="p">)</span><span class="w">
                  </span><span class="p">}</span><span class="w">
    
                  </span><span class="c1"># if not all chromosomes should be displayed, subset input
</span><span class="w">                  </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">is.na</span><span class="p">(</span><span class="n">chrs</span><span class="p">)[</span><span class="m">1</span><span class="p">])</span><span class="w"> </span><span class="p">{</span><span class="w">
                    </span><span class="n">input</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">input</span><span class="p">[</span><span class="nf">as.character</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">chr</span><span class="p">)</span><span class="w"> </span><span class="o">%in%</span><span class="w"> </span><span class="n">chrs</span><span class="p">,</span><span class="w"> </span><span class="p">]</span><span class="w">
                  </span><span class="p">}</span><span class="w">

                  </span><span class="c1"># if there is more than one LOD column, gather input
</span><span class="w">                  </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">any</span><span class="p">(</span><span class="n">colnames</span><span class="p">(</span><span class="n">input</span><span class="p">)</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="s2">"lod"</span><span class="p">))</span><span class="w"> </span><span class="p">{</span><span class="w">
                    </span><span class="n">input</span><span class="o">$</span><span class="n">lod</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">input</span><span class="p">[,</span><span class="w"> </span><span class="n">grep</span><span class="p">(</span><span class="s2">"lod"</span><span class="p">,</span><span class="w"> </span><span class="n">colnames</span><span class="p">(</span><span class="n">input</span><span class="p">))]</span><span class="w">
                  </span><span class="p">}</span><span class="w">

                  </span><span class="c1"># if no number of columns for facetting is defined, plot all in one row
</span><span class="w">                  </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="nf">is.na</span><span class="p">(</span><span class="n">ncol</span><span class="p">))</span><span class="w"> </span><span class="p">{</span><span class="w">
                    </span><span class="n">ncol</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="nf">length</span><span class="p">(</span><span class="n">unique</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">chr</span><span class="p">))</span><span class="w">
                  </span><span class="p">}</span><span class="w">

                  </span><span class="c1"># if labels are set and there is no name column, set from rownames
</span><span class="w">                  </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">is.na</span><span class="p">(</span><span class="n">labels</span><span class="p">)[</span><span class="m">1</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="nf">is.null</span><span class="p">(</span><span class="n">labels</span><span class="o">$</span><span class="n">name</span><span class="p">))</span><span class="w"> </span><span class="p">{</span><span class="w">
                      </span><span class="n">labels</span><span class="o">$</span><span class="n">name</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">rownames</span><span class="p">(</span><span class="n">labels</span><span class="p">)</span><span class="w">
                    </span><span class="p">}</span><span class="w">
                  </span><span class="p">}</span><span class="w">

                  </span><span class="c1"># plot input data frame position and LOD score
</span><span class="w">                  </span><span class="n">plot</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">ggplot</span><span class="p">(</span><span class="n">input</span><span class="p">,</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">pos</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">lod</span><span class="p">))</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="p">{</span><span class="w">

                    </span><span class="c1"># if LOD threshold is given, plot as horizontal line
</span><span class="w">                    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">is.na</span><span class="p">(</span><span class="n">lod</span><span class="p">)[</span><span class="m">1</span><span class="p">]</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="nf">length</span><span class="p">(</span><span class="n">lod</span><span class="p">)</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><span class="n">geom_hline</span><span class="p">(</span><span class="n">yintercept</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">lod</span><span class="p">,</span><span class="w"> </span><span class="n">linetype</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"dashed"</span><span class="p">)</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</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="o">!</span><span class="nf">is.na</span><span class="p">(</span><span class="n">lod</span><span class="p">)[</span><span class="m">1</span><span class="p">]</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="nf">length</span><span class="p">(</span><span class="n">lod</span><span class="p">)</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><span class="n">geom_hline</span><span class="p">(</span><span class="n">data</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">lod</span><span class="p">,</span><span class="w"> </span><span class="n">aes</span><span class="p">(</span><span class="n">yintercept</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">lod</span><span class="p">,</span><span class="w"> </span><span class="n">linetype</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">group</span><span class="p">))</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="p">{</span><span class="w">

                    </span><span class="c1"># plot rug on bottom, if TRUE
</span><span class="w">                    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">rug</span><span class="p">)</span><span class="w"> </span><span class="n">geom_rug</span><span class="p">(</span><span class="n">size</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.1</span><span class="p">,</span><span class="w"> </span><span class="n">sides</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"b"</span><span class="p">)</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="p">{</span><span class="w">

                    </span><span class="c1"># if input has column method but not group, plot line and color by method
</span><span class="w">                    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">method</span><span class="p">)</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">group</span><span class="p">))</span><span class="w"> </span><span class="n">geom_line</span><span class="p">(</span><span class="n">aes</span><span class="p">(</span><span class="n">color</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">method</span><span class="p">),</span><span class="w"> </span><span class="n">size</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="n">alpha</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.6</span><span class="p">)</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="p">{</span><span class="w">

                    </span><span class="c1"># if input has column group but not method, plot line and color by group
</span><span class="w">                    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">group</span><span class="p">)</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">method</span><span class="p">))</span><span class="w"> </span><span class="n">geom_line</span><span class="p">(</span><span class="n">aes</span><span class="p">(</span><span class="n">color</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">group</span><span class="p">),</span><span class="w"> </span><span class="n">size</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="n">alpha</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.6</span><span class="p">)</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="p">{</span><span class="w">

                    </span><span class="c1"># if input has columns method and group, plot line and color by method & linetype by group
</span><span class="w">                    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">group</span><span class="p">)</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="o">!</span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">method</span><span class="p">))</span><span class="w"> </span><span class="n">geom_line</span><span class="p">(</span><span class="n">aes</span><span class="p">(</span><span class="n">color</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">method</span><span class="p">,</span><span class="w"> </span><span class="n">linetype</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">group</span><span class="p">),</span><span class="w"> </span><span class="n">size</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="n">alpha</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.6</span><span class="p">)</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="p">{</span><span class="w">

                    </span><span class="c1"># set linetype, if input has columns method and group
</span><span class="w">                    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">group</span><span class="p">)</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="o">!</span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">method</span><span class="p">))</span><span class="w"> </span><span class="n">scale_linetype_manual</span><span class="p">(</span><span class="n">values</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">"solid"</span><span class="p">,</span><span class="w"> </span><span class="s2">"twodash"</span><span class="p">,</span><span class="w"> </span><span class="s2">"dotted"</span><span class="p">))</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="p">{</span><span class="w">

                    </span><span class="c1"># if input has neither columns method nor group, plot black line
</span><span class="w">                    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">group</span><span class="p">)</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="nf">is.null</span><span class="p">(</span><span class="n">input</span><span class="o">$</span><span class="n">method</span><span class="p">))</span><span class="w"> </span><span class="n">geom_line</span><span class="p">(</span><span class="n">size</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="n">alpha</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.6</span><span class="p">)</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="p">{</span><span class="w">

                    </span><span class="c1"># if QTL positions are given in labels df, plot as point...
</span><span class="w">                    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">is.na</span><span class="p">(</span><span class="n">labels</span><span class="p">)[</span><span class="m">1</span><span class="p">])</span><span class="w"> </span><span class="n">geom_point</span><span class="p">(</span><span class="n">data</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">labels</span><span class="p">,</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">pos</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">lod</span><span class="p">))</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="p">{</span><span class="w">

                    </span><span class="c1"># ... and plot name as text with ggrepel to avoid overlapping
</span><span class="w">                    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="nf">is.na</span><span class="p">(</span><span class="n">labels</span><span class="p">)[</span><span class="m">1</span><span class="p">])</span><span class="w"> </span><span class="n">geom_text_repel</span><span class="p">(</span><span class="n">data</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">labels</span><span class="p">,</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">pos</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">lod</span><span class="p">,</span><span class="w"> </span><span class="n">label</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">name</span><span class="p">),</span><span class="w"> </span><span class="n">nudge_y</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.5</span><span class="p">)</span><span class="w">
                  </span><span class="p">}</span><span class="w"> </span><span class="o">+</span><span class="w"> 
                    </span><span class="c1"># facet by chromosome
</span><span class="w">                    </span><span class="n">facet_wrap</span><span class="p">(</span><span class="o">~</span><span class="w"> </span><span class="n">chr</span><span class="p">,</span><span class="w"> </span><span class="n">ncol</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">ncol</span><span class="p">,</span><span class="w"> </span><span class="n">scales</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"free_x"</span><span class="p">)</span><span class="w"> </span><span class="o">+</span><span class="w">
                    </span><span class="c1"># minimal plotting theme
</span><span class="w">                    </span><span class="n">theme_minimal</span><span class="p">()</span><span class="w"> </span><span class="o">+</span><span class="w">
                    </span><span class="c1"># increase strip title size
</span><span class="w">                    </span><span class="n">theme</span><span class="p">(</span><span class="n">strip.text</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">element_text</span><span class="p">(</span><span class="n">face</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"bold"</span><span class="p">,</span><span class="w"> </span><span class="n">size</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">12</span><span class="p">))</span><span class="w"> </span><span class="o">+</span><span class="w">
                    </span><span class="c1"># use RcolorBrewer palette
</span><span class="w">                    </span><span class="n">scale_color_brewer</span><span class="p">(</span><span class="n">palette</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"Set1"</span><span class="p">)</span><span class="w"> </span><span class="o">+</span><span class="w">
                    </span><span class="c1"># Change plot labels
</span><span class="w">                    </span><span class="n">labs</span><span class="p">(</span><span class="n">x</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"Chromosome"</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="s2">"LOD"</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="s2">""</span><span class="p">,</span><span class="w">
                         </span><span class="n">linetype</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="n">print</span><span class="p">(</span><span class="n">plot</span><span class="p">)</span><span class="w">
</span><span class="p">}</span><span class="w">
</span>

The first example uses the hyper data set and builds a simple QTL model with three modeling functions: the EM algorithm, Haley-Knott regression and multiple imputation. The genome wide LOD threshold is calculated with permutation. Feeding this LOD threshold into the summary output gives us the markers with a significant phenotype association (i.e. the QTL).

<span class="n">data</span><span class="p">(</span><span class="n">hyper</span><span class="p">)</span><span class="w">

</span><span class="n">hyper</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">est.rf</span><span class="p">(</span><span class="n">hyper</span><span class="p">)</span><span class="w">

</span><span class="n">hyper</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">calc.genoprob</span><span class="p">(</span><span class="n">hyper</span><span class="p">,</span><span class="w"> </span><span class="n">step</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><span class="n">error.prob</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.01</span><span class="p">)</span><span class="w">
</span><span class="n">out.em</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">scanone</span><span class="p">(</span><span class="n">hyper</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="s2">"normal"</span><span class="p">)</span><span class="w">
</span><span class="n">out.hk</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">scanone</span><span class="p">(</span><span class="n">hyper</span><span class="p">,</span><span class="w"> </span><span class="n">method</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"hk"</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="s2">"normal"</span><span class="p">)</span><span class="w">

</span><span class="n">hyper</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">sim.geno</span><span class="p">(</span><span class="n">hyper</span><span class="p">,</span><span class="w"> </span><span class="n">step</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="n">n.draws</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">16</span><span class="p">,</span><span class="w"> </span><span class="n">error.prob</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.01</span><span class="p">)</span><span class="w">
</span><span class="n">out.imp</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">scanone</span><span class="p">(</span><span class="n">hyper</span><span class="p">,</span><span class="w"> </span><span class="n">method</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"imp"</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="s2">"normal"</span><span class="p">)</span><span class="w">

</span><span class="n">operm.hk</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">scanone</span><span class="p">(</span><span class="n">hyper</span><span class="p">,</span><span class="w"> </span><span class="n">pheno.col</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><span class="n">model</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"normal"</span><span class="p">,</span><span class="w"> </span><span class="n">method</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"hk"</span><span class="p">,</span><span class="w"> </span><span class="n">n.perm</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">1000</span><span class="p">,</span><span class="w"> </span><span class="n">verbose</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">lod_threshold</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">summary</span><span class="p">(</span><span class="n">operm.hk</span><span class="p">,</span><span class="w"> </span><span class="n">alpha</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.05</span><span class="p">)</span><span class="w">
</span><span class="n">labels_df</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">summary</span><span class="p">(</span><span class="n">out.hk</span><span class="p">,</span><span class="w"> </span><span class="n">perms</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">operm.hk</span><span class="p">,</span><span class="w"> </span><span class="n">alpha</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.05</span><span class="p">,</span><span class="w"> </span><span class="n">pvalues</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">TRUE</span><span class="p">))</span><span class="w">
</span>

The most basic version of my plotting function takes as input a data frame from the scanone() function: It needs to have a column with chromosome information, position and LOD score. Without defining anything else, the function will return a line plot of the LOD scores for each marker position on each chromosome in the input data.

<span class="n">qtl_plot</span><span class="p">(</span><span class="n">input</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">out.hk</span><span class="p">)</span><span class="w">
</span>

If we want to compare the output from e.g. different QTL models, we can follow the tidyverse principle by combining the output data frames in long format and adding a column with the respective method names. If the input has a method column, it will color the lines accordingly. We can also specify, which chromosomes to display (here, 1 and 4) and give a LOD threshold to plot as a horizontal dashed line. If we want to see each marker position, we can specify rug = TRUE to plot them underneath the line plots.

And we can also plot the markers with a significant phenotype association (the QTL) by adding a data frame with the labeling information to the function call. This data frame needs to contain the chromosome, position and LOD score information. If there is no name column, it will take the row names (which usually give the marker name) as labels.

<span class="n">qtl_plot</span><span class="p">(</span><span class="n">input</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">rbind</span><span class="p">(</span><span class="n">data.frame</span><span class="p">(</span><span class="n">out.em</span><span class="p">,</span><span class="w"> </span><span class="n">method</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"EM algorithm"</span><span class="p">),</span><span class="w"> 
                       </span><span class="n">data.frame</span><span class="p">(</span><span class="n">out.hk</span><span class="p">,</span><span class="w"> </span><span class="n">method</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"Haley-Knott regression"</span><span class="p">),</span><span class="w"> 
                       </span><span class="n">data.frame</span><span class="p">(</span><span class="n">out.imp</span><span class="p">,</span><span class="w"> </span><span class="n">method</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"Multiple imputation"</span><span class="p">)),</span><span class="w"> 
         </span><span class="n">chrs</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">4</span><span class="p">),</span><span class="w"> 
         </span><span class="n">lod</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">lod_threshold</span><span class="p">[</span><span class="m">1</span><span class="p">],</span><span class="w"> 
         </span><span class="n">rug</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">TRUE</span><span class="p">,</span><span class="w"> 
         </span><span class="n">labels</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">labels_df</span><span class="p">)</span><span class="w">
</span>

The second example I want to show is from the listeria data. This time, a slightly more complex 2-part model is built.

<span class="n">data</span><span class="p">(</span><span class="n">listeria</span><span class="p">)</span><span class="w">

</span><span class="n">listeria</span><span class="o">$</span><span class="n">pheno</span><span class="o">$</span><span class="n">logSurv</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="nf">log</span><span class="p">(</span><span class="n">listeria</span><span class="o">$</span><span class="n">pheno</span><span class="p">[,</span><span class="m">1</span><span class="p">])</span><span class="w">
</span><span class="n">listeria</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">est.rf</span><span class="p">(</span><span class="n">listeria</span><span class="p">)</span><span class="w">

</span><span class="n">newmap</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">est.map</span><span class="p">(</span><span class="n">listeria</span><span class="p">,</span><span class="w"> </span><span class="n">error.prob</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.01</span><span class="p">)</span><span class="w">
</span><span class="n">listeria</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">replace.map</span><span class="p">(</span><span class="n">listeria</span><span class="p">,</span><span class="w"> </span><span class="n">newmap</span><span class="p">)</span><span class="w">

</span><span class="n">listeria</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">calc.errorlod</span><span class="p">(</span><span class="n">listeria</span><span class="p">,</span><span class="w"> </span><span class="n">error.prob</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.01</span><span class="p">)</span><span class="w">

</span><span class="n">listeria</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">calc.genoprob</span><span class="p">(</span><span class="n">listeria</span><span class="p">,</span><span class="w"> </span><span class="n">step</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">2</span><span class="p">)</span><span class="w">
</span><span class="n">out.2p</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">scanone</span><span class="p">(</span><span class="n">listeria</span><span class="p">,</span><span class="w"> </span><span class="n">pheno.col</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">3</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="s2">"2part"</span><span class="p">,</span><span class="w"> </span><span class="n">upper</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">TRUE</span><span class="p">)</span><span class="w">

</span><span class="n">operm.2p</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">scanone</span><span class="p">(</span><span class="n">listeria</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="s2">"2part"</span><span class="p">,</span><span class="w"> </span><span class="n">pheno.col</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">3</span><span class="p">,</span><span class="w"> </span><span class="n">upper</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">TRUE</span><span class="p">,</span><span class="w"> </span><span class="n">n.perm</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">25</span><span class="p">,</span><span class="w"> </span><span class="n">verbose</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">lod_threshold</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">summary</span><span class="p">(</span><span class="n">operm.2p</span><span class="p">,</span><span class="w"> </span><span class="n">alpha</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.05</span><span class="p">)</span><span class="w">

</span><span class="n">labels_df</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">summary</span><span class="p">(</span><span class="n">out.2p</span><span class="p">,</span><span class="w"> </span><span class="n">format</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"allpeaks"</span><span class="p">,</span><span class="w"> </span><span class="n">perms</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">operm.2p</span><span class="p">,</span><span class="w"> </span><span class="n">alpha</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.05</span><span class="p">,</span><span class="w"> </span><span class="n">pvalues</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">TRUE</span><span class="p">))</span><span class="w">
</span>

The output from a 2-part model looks different than from the other model functions: Instead of only one LOD score column, we now get three columns named lod.pu.mu, lod.p and lod.mu (check back with the tutorial if you want to know what they mean). The plotting function is built so that if you define the model parameter as “2part”, each of the three LOD columns will be plotted in a different color.

If we also want to plot the QTL labels, we need to prepare the label data frame before feeding it into the plotting function:

<span class="n">p.mu</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">labels_df</span><span class="p">[,</span><span class="w"> </span><span class="m">1</span><span class="o">:</span><span class="m">4</span><span class="p">]</span><span class="w">
</span><span class="n">colnames</span><span class="p">(</span><span class="n">p.mu</span><span class="p">)[</span><span class="m">3</span><span class="p">]</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="s2">"lod"</span><span class="w">
</span><span class="n">p.mu</span><span class="o">$</span><span class="n">name</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="s2">"p.mu"</span><span class="w">

</span><span class="n">p</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">labels_df</span><span class="p">[,</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">5</span><span class="o">:</span><span class="m">7</span><span class="p">)]</span><span class="w">
</span><span class="n">colnames</span><span class="p">(</span><span class="n">p</span><span class="p">)[</span><span class="m">3</span><span class="p">]</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="s2">"lod"</span><span class="w">
</span><span class="n">p</span><span class="o">$</span><span class="n">name</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="s2">"p"</span><span class="w">

</span><span class="n">mu</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">labels_df</span><span class="p">[,</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">8</span><span class="o">:</span><span class="m">10</span><span class="p">)]</span><span class="w">
</span><span class="n">colnames</span><span class="p">(</span><span class="n">mu</span><span class="p">)[</span><span class="m">3</span><span class="p">]</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="s2">"lod"</span><span class="w">
</span><span class="n">mu</span><span class="o">$</span><span class="n">name</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="s2">"mu"</span><span class="w">

</span><span class="n">labels_df</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">rbind</span><span class="p">(</span><span class="n">p.mu</span><span class="p">,</span><span class="w"> </span><span class="n">p</span><span class="p">,</span><span class="w"> </span><span class="n">mu</span><span class="p">)</span><span class="w">
</span>

If the chromosomes are not in the right order, we can reorder the factor labels:

<span class="n">f</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">"1"</span><span class="p">,</span><span class="w"> </span><span class="s2">"2"</span><span class="p">,</span><span class="w"> </span><span class="s2">"3"</span><span class="p">,</span><span class="w"> </span><span class="s2">"4"</span><span class="p">,</span><span class="w"> </span><span class="s2">"5"</span><span class="p">,</span><span class="w"> </span><span class="s2">"6"</span><span class="p">,</span><span class="w"> </span><span class="s2">"7"</span><span class="p">,</span><span class="w"> </span><span class="s2">"8"</span><span class="p">,</span><span class="w"> </span><span class="s2">"9"</span><span class="p">,</span><span class="w"> </span><span class="s2">"10"</span><span class="p">,</span><span class="w"> </span><span class="s2">"11"</span><span class="p">,</span><span class="w"> </span><span class="s2">"12"</span><span class="p">,</span><span class="w"> </span><span class="s2">"13"</span><span class="p">,</span><span class="w"> </span><span class="s2">"14"</span><span class="p">,</span><span class="w"> </span><span class="s2">"15"</span><span class="p">,</span><span class="w"> </span><span class="s2">"16"</span><span class="p">,</span><span class="w"> </span><span class="s2">"17"</span><span class="p">,</span><span class="w"> </span><span class="s2">"18"</span><span class="p">,</span><span class="w"> </span><span class="s2">"19"</span><span class="p">,</span><span class="w"> </span><span class="s2">"X"</span><span class="p">)</span><span class="w">
</span><span class="n">out.2p</span><span class="o">$</span><span class="n">chr</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">factor</span><span class="p">(</span><span class="n">out.2p</span><span class="o">$</span><span class="n">chr</span><span class="p">,</span><span class="w"> </span><span class="n">levels</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">f</span><span class="p">)</span><span class="w">
</span><span class="n">labels_df</span><span class="o">$</span><span class="n">chr</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">factor</span><span class="p">(</span><span class="n">labels_df</span><span class="o">$</span><span class="n">chr</span><span class="p">,</span><span class="w"> </span><span class="n">levels</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">f</span><span class="p">)</span><span class="w">
</span>

The parameter ncol let’s you define how many facet columns you want to have. Here, I want two rows and three columns to display the six chromosomes I chose for plotting:

<span class="n">qtl_plot</span><span class="p">(</span><span class="n">input</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">out.2p</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="s2">"2part"</span><span class="p">,</span><span class="w">
         </span><span class="n">chrs</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">3</span><span class="p">,</span><span class="w"> </span><span class="m">5</span><span class="o">:</span><span class="m">6</span><span class="p">,</span><span class="w"> </span><span class="m">13</span><span class="p">,</span><span class="w"> </span><span class="m">15</span><span class="p">),</span><span class="w">
         </span><span class="n">lod</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">lod_threshold</span><span class="p">[</span><span class="m">1</span><span class="p">],</span><span class="w"> 
         </span><span class="n">rug</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">TRUE</span><span class="p">,</span><span class="w"> 
         </span><span class="n">ncol</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">3</span><span class="p">,</span><span class="w">
         </span><span class="n">labels</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">labels_df</span><span class="p">)</span><span class="w">
</span>

Binary models are again treated the same way as e.g. normal models:

<span class="n">y</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">listeria</span><span class="o">$</span><span class="n">pheno</span><span class="o">$</span><span class="n">logSurv</span><span class="w">
</span><span class="n">my</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="nf">max</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="w"> </span><span class="n">na.rm</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">TRUE</span><span class="p">)</span><span class="w">
</span><span class="n">z</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="nf">as.numeric</span><span class="p">(</span><span class="n">y</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="n">my</span><span class="p">)</span><span class="w">
</span><span class="n">y</span><span class="p">[</span><span class="n">y</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="n">my</span><span class="p">]</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="kc">NA</span><span class="w">
</span><span class="n">listeria</span><span class="o">$</span><span class="n">pheno</span><span class="o">$</span><span class="n">logSurv2</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">y</span><span class="w">
</span><span class="n">listeria</span><span class="o">$</span><span class="n">pheno</span><span class="o">$</span><span class="n">binary</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">z</span><span class="w">

</span><span class="n">out.p</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">scanone</span><span class="p">(</span><span class="n">listeria</span><span class="p">,</span><span class="w"> </span><span class="n">pheno.col</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">model</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"binary"</span><span class="p">)</span><span class="w">
</span><span class="n">out.np1</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">scanone</span><span class="p">(</span><span class="n">listeria</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="s2">"np"</span><span class="p">,</span><span class="w"> </span><span class="n">ties.random</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="kc">TRUE</span><span class="p">)</span><span class="w">
</span><span class="n">out.np2</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">scanone</span><span class="p">(</span><span class="n">listeria</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="s2&qu...

To leave a comment for the author, please follow the link and comment on their blog: Shirin's playgRound.

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)