PVA: Publication Viability Analysis, round 3

[This article was first published on Peter Solymos - R related posts, 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.

A friend and colleague of mine, Péter Batáry
has circulated news from Nature
magazine about the EU freezing innovation funds to Bulgaria.
The article had a figure about publication trends for
Bulgaria, compared with Romania and Hungary.
As I have blogged about such trends in ecology before
(here and
here),
I felt the need to update my PVA models
with two years worth of data from WoS.

After downloading the yearly publications numbers
using filters ADDRESS=HUNGARY; CATEGORIES=ECOLOGY,
I started where I left off few years ago. I fit Ricker growth model
to two time intervals of the data: 1978–1997, and 1998–2017.

The R code below uses the PVAClone package
that I wrote with Khurram Nadeem,
and is based on fitting state-space models using
MCMC and data cloning with JAGS.
The other intrval package is pretty new but handy little helper
(see related posts here)

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

</span><span class="c1">## the data from WoS</span><span class="w">
</span><span class="n">x</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">structure</span><span class="p">(</span><span class="nf">list</span><span class="p">(</span><span class="n">years</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">1973</span><span class="o">:</span><span class="m">2017</span><span class="p">,</span><span class="w"> </span><span class="n">records</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">0</span><span class="p">,</span><span class="w"> </span><span class="m">4</span><span class="p">,</span><span class="w"> </span><span class="m">0</span><span class="p">,</span><span class="w"> </span><span class="m">0</span><span class="p">,</span><span class="w">
    </span><span class="m">6</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">5</span><span class="p">,</span><span class="w"> </span><span class="m">4</span><span class="p">,</span><span class="w"> </span><span class="m">7</span><span class="p">,</span><span class="w"> </span><span class="m">5</span><span class="p">,</span><span class="w"> </span><span class="m">7</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="p">,</span><span class="w"> </span><span class="m">9</span><span class="p">,</span><span class="w"> </span><span class="m">11</span><span class="p">,</span><span class="w"> </span><span class="m">20</span><span class="p">,</span><span class="w"> </span><span class="m">8</span><span class="p">,</span><span class="w"> </span><span class="m">10</span><span class="p">,</span><span class="w"> </span><span class="m">15</span><span class="p">,</span><span class="w"> </span><span class="m">29</span><span class="p">,</span><span class="w"> </span><span class="m">24</span><span class="p">,</span><span class="w"> </span><span class="m">53</span><span class="p">,</span><span class="w">
    </span><span class="m">12</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">30</span><span class="p">,</span><span class="w"> </span><span class="m">32</span><span class="p">,</span><span class="w"> </span><span class="m">36</span><span class="p">,</span><span class="w"> </span><span class="m">45</span><span class="p">,</span><span class="w"> </span><span class="m">39</span><span class="p">,</span><span class="w"> </span><span class="m">42</span><span class="p">,</span><span class="w"> </span><span class="m">43</span><span class="p">,</span><span class="w"> </span><span class="m">50</span><span class="p">,</span><span class="w"> </span><span class="m">62</span><span class="p">,</span><span class="w"> </span><span class="m">95</span><span class="p">,</span><span class="w"> </span><span class="m">106</span><span class="p">,</span><span class="w"> </span><span class="m">113</span><span class="p">,</span><span class="w"> </span><span class="m">83</span><span class="p">,</span><span class="w">
    </span><span class="m">108</span><span class="p">,</span><span class="w"> </span><span class="m">99</span><span class="p">,</span><span class="w"> </span><span class="m">89</span><span class="p">,</span><span class="w"> </span><span class="m">117</span><span class="p">,</span><span class="w"> </span><span class="m">111</span><span class="p">,</span><span class="w"> </span><span class="m">134</span><span class="p">,</span><span class="w"> </span><span class="m">127</span><span class="p">)),</span><span class="w"> </span><span class="n">.Names</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">"years"</span><span class="p">,</span><span class="w"> </span><span class="s2">"records"</span><span class="w">
    </span><span class="p">),</span><span class="w"> </span><span class="n">row.names</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nf">c</span><span class="p">(</span><span class="kc">NA</span><span class="p">,</span><span class="w"> </span><span class="m">45L</span><span class="p">),</span><span class="w"> </span><span class="n">class</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="s2">"data.frame"</span><span class="p">)</span><span class="w">

</span><span class="c1">## fit the 2 models</span><span class="w">
</span><span class="n">ncl</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="m">10</span><span class="w"> </span><span class="c1"># number of clones</span><span class="w">
</span><span class="n">m</span><span class="m">1</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">pva</span><span class="p">(</span><span class="n">x</span><span class="o">$</span><span class="n">records</span><span class="p">[</span><span class="n">x</span><span class="o">$</span><span class="n">years</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">1978</span><span class="p">,</span><span class="w"> </span><span class="m">1997</span><span class="p">)],</span><span class="w"> </span><span class="n">ricker</span><span class="p">(</span><span class="s2">"none"</span><span class="p">),</span><span class="w"> </span><span class="n">ncl</span><span class="p">)</span><span class="w">
</span><span class="n">m</span><span class="m">2</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">pva</span><span class="p">(</span><span class="n">x</span><span class="o">$</span><span class="n">records</span><span class="p">[</span><span class="n">x</span><span class="o">$</span><span class="n">years</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">1998</span><span class="p">,</span><span class="w"> </span><span class="m">2017</span><span class="p">)],</span><span class="w"> </span><span class="n">ricker</span><span class="p">(</span><span class="s2">"none"</span><span class="p">),</span><span class="w"> </span><span class="n">ncl</span><span class="p">)</span><span class="w">

</span><span class="c1">## organize estimates</span><span class="w">
</span><span class="n">cf</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">t</span><span class="p">(</span><span class="n">sapply</span><span class="p">(</span><span class="nf">list</span><span class="p">(</span><span class="n">early</span><span class="o">=</span><span class="n">m</span><span class="m">1</span><span class="p">,</span><span class="w"> </span><span class="n">late</span><span class="o">=</span><span class="n">m</span><span class="m">2</span><span class="p">),</span><span class="w"> </span><span class="n">coef</span><span class="p">)))</span><span class="w">
</span><span class="n">cf</span><span class="o">$</span><span class="n">K</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">with</span><span class="p">(</span><span class="n">cf</span><span class="p">,</span><span class="w"> </span><span class="o">-</span><span class="n">a</span><span class="o">/</span><span class="n">b</span><span class="p">)</span><span class="w">

</span><span class="c1">## growth curve: early period</span><span class="w">
</span><span class="n">yr1</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="m">1978</span><span class="o">:</span><span class="m">1997</span><span class="w">
</span><span class="n">pr1</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">numeric</span><span class="p">(</span><span class="nf">length</span><span class="p">(</span><span class="n">yr1</span><span class="p">))</span><span class="w">
</span><span class="n">pr1</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">log</span><span class="p">(</span><span class="n">x</span><span class="o">$</span><span class="n">records</span><span class="p">[</span><span class="n">x</span><span class="o">$</span><span class="n">years</span><span class="o">==</span><span class="m">1978</span><span class="p">])</span><span class="w">
</span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">i</span><span class="w"> </span><span class="k">in</span><span class="w"> </span><span class="m">2</span><span class="o">:</span><span class="nf">length</span><span class="p">(</span><span class="n">pr1</span><span class="p">))</span><span class="w">
    </span><span class="n">pr1</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">pr1</span><span class="p">[</span><span class="n">i</span><span class="m">-1</span><span class="p">]</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">cf</span><span class="p">[</span><span class="s2">"early"</span><span class="p">,</span><span class="w"> </span><span class="s2">"a"</span><span class="p">]</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">cf</span><span class="p">[</span><span class="s2">"early"</span><span class="p">,</span><span class="w"> </span><span class="s2">"b"</span><span class="p">]</span><span class="o">*</span><span class="nf">exp</span><span class="p">(</span><span class="n">pr1</span><span class="p">[</span><span class="n">i</span><span class="m">-1</span><span class="p">])</span><span class="w">
</span><span class="n">pr1</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="nf">exp</span><span class="p">(</span><span class="n">pr1</span><span class="p">)</span><span class="w">

</span><span class="c1">## growth curve: late period</span><span class="w">
</span><span class="n">yr2</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="m">1998</span><span class="o">:</span><span class="m">2017</span><span class="w">
</span><span class="n">pr2</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">numeric</span><span class="p">(</span><span class="nf">length</span><span class="p">(</span><span class="n">yr2</span><span class="p">))</span><span class="w">
</span><span class="n">pr2</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">log</span><span class="p">(</span><span class="n">x</span><span class="o">$</span><span class="n">records</span><span class="p">[</span><span class="n">x</span><span class="o">$</span><span class="n">years</span><span class="o">==</span><span class="m">1998</span><span class="p">])</span><span class="w">
</span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">i</span><span class="w"> </span><span class="k">in</span><span class="w"> </span><span class="m">2</span><span class="o">:</span><span class="nf">length</span><span class="p">(</span><span class="n">pr2</span><span class="p">))</span><span class="w">
    </span><span class="n">pr2</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">pr2</span><span class="p">[</span><span class="n">i</span><span class="m">-1</span><span class="p">]</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">cf</span><span class="p">[</span><span class="s2">"late"</span><span class="p">,</span><span class="w"> </span><span class="s2">"a"</span><span class="p">]</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">cf</span><span class="p">[</span><span class="s2">"late"</span><span class="p">,</span><span class="w"> </span><span class="s2">"b"</span><span class="p">]</span><span class="o">*</span><span class="nf">exp</span><span class="p">(</span><span class="n">pr2</span><span class="p">[</span><span class="n">i</span><span class="m">-1</span><span class="p">])</span><span class="w">
</span><span class="n">pr2</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="nf">exp</span><span class="p">(</span><span class="n">pr2</span><span class="p">)</span><span class="w">

</span><span class="c1">## and finally the figure using base graphics</span><span class="w">
</span><span class="n">op</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">par</span><span class="p">(</span><span class="n">las</span><span class="o">=</span><span class="m">2</span><span class="p">)</span><span class="w">
</span><span class="n">barplot</span><span class="p">(</span><span class="n">x</span><span class="o">$</span><span class="n">records</span><span class="p">,</span><span class="w"> </span><span class="n">names.arg</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">x</span><span class="o">$</span><span class="n">years</span><span class="p">,</span><span class="w"> </span><span class="n">space</span><span class="o">=</span><span class="m">0</span><span class="p">,</span><span class="w">
    </span><span class="n">ylab</span><span class="o">=</span><span class="s2">"# of publications"</span><span class="p">,</span><span class="w"> </span><span class="n">xlab</span><span class="o">=</span><span class="s2">"years"</span><span class="p">,</span><span class="w">
    </span><span class="n">col</span><span class="o">=</span><span class="n">ifelse</span><span class="p">(</span><span class="n">x</span><span class="o">$</span><span class="n">years</span><span class="w"> </span><span class="o"><</span><span class="w"> </span><span class="m">1998</span><span class="p">,</span><span class="w"> </span><span class="s2">"grey"</span><span class="p">,</span><span class="w"> </span><span class="s2">"gold"</span><span class="p">))</span><span class="w">
</span><span class="n">lines</span><span class="p">(</span><span class="n">yr1</span><span class="o">-</span><span class="nf">min</span><span class="p">(</span><span class="n">x</span><span class="o">$</span><span class="n">years</span><span class="p">)</span><span class="m">+0.5</span><span class="p">,</span><span class="w"> </span><span class="n">pr1</span><span class="p">,</span><span class="w"> </span><span class="n">col</span><span class="o">=</span><span class="m">4</span><span class="p">)</span><span class="w">
</span><span class="n">abline</span><span class="p">(</span><span class="n">h</span><span class="o">=</span><span class="n">cf</span><span class="p">[</span><span class="s2">"early"</span><span class="p">,</span><span class="w"> </span><span class="s2">"K"</span><span class="p">],</span><span class="w"> </span><span class="n">col</span><span class="o">=</span><span class="m">4</span><span class="p">,</span><span class="w"> </span><span class="n">lty</span><span class="o">=</span><span class="m">3</span><span class="p">)</span><span class="w">
</span><span class="n">lines</span><span class="p">(</span><span class="n">yr2</span><span class="o">-</span><span class="nf">min</span><span class="p">(</span><span class="n">x</span><span class="o">$</span><span class="n">years</span><span class="p">)</span><span class="m">+0.5</span><span class="p">,</span><span class="w"> </span><span class="n">pr2</span><span class="p">,</span><span class="w"> </span><span class="n">col</span><span class="o">=</span><span class="m">2</span><span class="p">)</span><span class="w">
</span><span class="n">abline</span><span class="p">(</span><span class="n">h</span><span class="o">=</span><span class="n">cf</span><span class="p">[</span><span class="s2">"late2017"</span><span class="p">,</span><span class="w"> </span><span class="s2">"K"</span><span class="p">],</span><span class="w"> </span><span class="n">col</span><span class="o">=</span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="n">lty</span><span class="o">=</span><span class="m">3</span><span class="p">)</span><span class="w">
</span><span class="n">par</span><span class="p">(</span><span class="n">op</span><span class="p">)</span><span class="w">
</span>

PVA

Here are the model parameters for the two Ricker models:

  a b sigma K
1978–1997 0.38 -0.03 0.60 13.85
1998–2017 0.21 0.00 0.16 119.00

The K carrying capacity used to be 100 based on
1998–2012 data, but now K = 119, which is
a significant improvement — heartfelt kudos to the ecologists in Hungary
(more papers please)!
The growth rate hasn’t changed (a = 0.21).
So we can conclude that if the rate remained constant
but carrying capacity increased, the change must be
related to resource availability
(i.e. increased funding, more jobs, improved infrastructure).

This is good news to me! Let me know what you think by leaving a comment below!

To leave a comment for the author, please follow the link and comment on their blog: Peter Solymos - R related posts.

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)