# How much of the world is woody?

<h2 class="pm-node nj-subtitle">Simulation &amp; Analysis</h2>

<p class="pm-node nj-authors">Richard G. FitzJohn, Matthew W. Pennell, Amy E. Zanne, Peter F. Stevens, David C. Tank, William K. Cornwell</p>

*This notebook was generated by the [Data/Processing notebook](https://nextjournal.com/mpd/wood-analysis), and runs on an environment generated therein.  Most plots converted to Plotly for interactivity.*

This is the full analysis that underlies our paper.

```r id=de97d7f8-ad8c-4e0e-af05-766d73ac6905
knitr::opts_chunk$set(tidy=FALSE)
```

```r id=d2e9b42b-efd4-43cc-b168-1b20ea6532e8
invisible(suppressWarnings(sample(1:250, 1, pr=rep(1, 250), replace=TRUE)))
```

# Preliminaries

Ensure that you have all the required packages installed by running

`make deps
`

Before running this script, be sure to have the required data in place.  If you have `make` installed, this can be done by running

`make data-processed
`

at the command line.  See `README.md` for more details.

See the file `data/README.md` for more information about the data that we use, and alternative ways of running things if you don't have `make` installed.

```r id=09c35830-3eae-4c8b-a561-3a7a07325e1c
library(diversitree)
source("R/load.R")
source("R/util.R")
```

Colours used throughout:

```r id=da9ce426-335a-4ba6-9e38-12d0c7d3c9a0
cols.methods <- c(binomial="#a63813",
                  hypergeometric="#4d697f")
cols.tree <- c(Monilophytes="#a63813",
               Gymnosperms="#21313b",
               BasalAngiosperms="#eeb911",
               Monocots="#204d14",
               Eudicots="#4d697f",
               Rest="gray15")
cols.woody <- c(Woody="#533d0c",
                Herbaceous="#799321")
cols.woody <- c(Woody="black",
                Herbaceous="white")
cols.shading <- "#eeb911"
cols.tropical <- c(tropical="#ea7518",
                   temperate="#62a184")
```

The 'load.woodiness.data.genus' function hides a reasonable amount of data cleaning required to load the data.  This mostly involves matching the woodiness data set from Zanne et al. to our species list (derived from The Plant List), cleaning up synonomies, and then collapsing down to genus.

The final object has columns

* genus, family, order -- taxonomic information
* W, V, H              -- number of species scored as woody, variable, herbaceous (respectively)
* N                    -- number of species in the genus
* K                    -- number of species with known state, after dropping all "variable" species
* p                    -- fraction of known species that are woody, after dropping all "variable" species

```r id=3ebb6dbe-6da3-46d2-806a-baa7e51bb2bd
dat.g <- load.woodiness.genus()
```

# Imputing the state of missing species

For each genus:

A. If the genus has a valid fraction of species (i.e. K > 0 so p is defined), then sample the number of species that are woody from either

* the binomial distribution (strong prior; assuming known species were sampled with replacement from the pool of species); or
* the hypergeometric distribution (weak prior; assuming that known species were sampled without replacement from the pool of species).

B. If the genus has no valid fraction of species (i.e., K == 0 so p is undefined), then sample from the emperical distribution of per-genus fractions.  We're going to feed data into this by taxonomic order, so this will come from the per-order distribution.

```r id=301f2420-acf7-405a-b971-dc5102143ebf
sim <- function(x, nrep, with.replacement=TRUE, p=1/20) {
  ## First, focus on cases where we have a valid estimate of the
  ## fraction of species that are woody (i.e., at least one known
  ## species).
  ok <- !is.na(x$p)

  w <- matrix(NA, nrow(x), nrep)

  ## A: genera with any known species
  if (with.replacement)
    w[ok,] <- x$W[ok] + rbinom(sum(ok), x$N[ok]-x$K[ok], x$W[ok]/x$K[ok])
  else
    w[ok,] <- t(sapply(which(ok), function(i)
                       rhyper2(nrep, x$H[i], x$W[i], x$N[i])))

  ## B: genera with no known species
  n.unk <- sum(!ok)
  w[!ok,] <- apply(w[ok,,drop=FALSE] / x$N[ok], 2, function(y)
                   rbinom(n.unk, x$N[!ok], quantile(y, runif(n.unk))))
  
  rownames(w) <- x$genus

  summarise.sim(w, x[c("order", "family", "genus",
                       "W", "V", "H", "N", "K")])
}
```

This collects up the results at different taxonomic levels.

```r id=f28dc23d-1cda-4e4e-966c-9d2cabef57cf
summarise <- function(x, p=1/20)
  structure(c(mean(x), quantile(x, c(p/2, 1-p/2))),
            names=c("mean", "lower", "upper"))
summarise.sim <- function(w, info) {
  order <- info$order[[1]]

  info.cols <- c("W", "V", "H", "N", "K")

  ## Genus is easy;
  w.g <- cbind(info, t(apply(w, 1, summarise)))

  ## Family is a pain:
  w.f <- do.call(rbind,
                 lapply(split(as.data.frame(w), info$family), colSums))
  w.f <- t(apply(w.f, 1, summarise))
  w.f <- data.frame(order=order,
                    aggregate(info[info.cols], info["family"], sum),
                    w.f, stringsAsFactors=TRUE)
  rownames(w.f) <- NULL
  
  ## Order is easy; we are guaranteed to have just one order here, so:
  w.o <- data.frame(order=order,
                    as.data.frame(t(colSums(info[info.cols]))),
                    t(summarise(colSums(w))), stringsAsFactors=FALSE)

  ret <- list(genus=w.g, family=w.f, order=w.o)
  attr(ret, "total") <- colSums(w)
  ret
}

rhyper2 <- function(nn, s0, s1, xn, fraction=FALSE) {
  x1 <- seq(s1, xn - s0)
  x0 <- xn - x1
  p1 <- dhyper(s1, x1, x0, s0+s1)
  p1 <- p1 / sum(p1)
  x1[sample(length(p1), nn, TRUE, p1)]
}

do.simulation <- function(dat.g, nrep, with.replacement) {
  f <- function(level) {
    ret <- do.call(rbind, lapply(res, "[[", level))
    rownames(ret) <- NULL
    ret[c("p.mean", "p.lower", "p.upper")] <-
      ret[c("mean", "lower", "upper")] / ret[["N"]]
    ret
  }

  res <- lapply(split(dat.g, dat.g$order),
                sim, nrep, with.replacement)
  total <- rowSums(sapply(res, attr, "total"))
  overall <- summarise(total)
  overall.p <- overall / sum(dat.g$N)

  list(genus=f("genus"), family=f("family"), order=f("order"),
       overall=overall, overall.p=overall.p, total=total)
}
```

```r id=88695a89-75ca-4475-8bf2-82f29530265b
set.seed(1)
if (!(exists("res.b") && exists("res.h"))) {
  res.b <- do.simulation(dat.g, 1000, TRUE)  # binomial - with replacement
  res.h <- do.simulation(dat.g, 1000, FALSE) # hypergeometric - without
}
```

Repeat this sampling proceedure under a more extreme classification; all variable species scored as woody:

```r id=493966e3-c829-46ba-9615-58571ef44711
if (!(exists("res.b.w") && exists("res.h.w"))) {
  dat.g.w <- load.woodiness.genus(extreme="woody")
  res.b.w <- do.simulation(dat.g.w, 1000, TRUE)  # binomial       & woody
  res.h.w <- do.simulation(dat.g.w, 1000, FALSE) # hypergeometric & woody
}
```

...and with all variable species scored as herbaceous:

```r id=ca8543ca-4ec0-44e8-a12f-7678c05de9bf
if (!(exists("res.b.h") && exists("res.h.h"))) {
  dat.g.h <- load.woodiness.genus(extreme="herbaceous")
  res.b.h <- do.simulation(dat.g.h, 1000, TRUE)  # binomial       & herby
  res.h.h <- do.simulation(dat.g.h, 1000, FALSE) # hypergeometric & herby
}
```

# Distribution of estimates of woodiness

This is the raw distribution; i.e., our estimate of the fraction of species that are woody and its estimate.

```r id=aa125966-03eb-4e80-ae86-e6818ebcdc17
suppressMessages(library(dplyr, quietly=TRUE))
suppressMessages(library(plotly, quietly=TRUE))

fig.distribution.raw <- function(res.b, res.h) {
  n.spp <- sum(res.b$order$N)
  p.b <- res.b$total / n.spp * 100
  p.h <- res.h$total / n.spp * 100
  
  r <- range(p.b, p.h)
  br <- seq(r[1], r[2], length.out=30)

  h.b <- hist(p.b, br, plot=FALSE)
  h.h <- hist(p.h, br, plot=FALSE)
  
  hg.b <- data.frame(c(h.b[c("mids","density")]))
  hg.h <- data.frame(c(h.h[c("mids","density")]))

  xlim <- c(42, 50)
  ylim <- range(h.b$density, h.h$density)

  cols <- cols.methods
  
  plot_ly(opacity = 0.6) %>%
  add_bars(data = hg.b, x = ~mids, y = ~density, 
    name = "binomial w/replacement",
    marker = list(color = cols[1],
                  line=list(color='black',width=1))) %>%
  add_bars(data = hg.h, x = ~mids, y = ~density,
    name = "hypergeometric w/o replacement",
    marker = list(color = cols[2],
                  line=list(color='black',width=1))) %>%
  layout(barmode = "overlay", 
         legend = list(x = 0.55, y = 0.95),
         yaxis = list(range=ylim, title = "Probability density"),
         xaxis = list(range=xlim, 
           title = "Percentage of woody species among all vascular plants"))
}
```

```r id=bdbc00be-b887-4a47-ba1b-e6b69fce4280
fig.distribution.raw(res.b, res.h)
```

[__out-0.plotly.json][nextjournal#output#bdbc00be-b887-4a47-ba1b-e6b69fce4280#__out-0.plotly.json]

Next, expand that plot to include the more extreme estimates (variable species as woody and variable species as herbaceous). This shows how much the classification errors affect the analysis.

```r id=8f34d550-b5f1-4714-897b-bea40be282d1
fig.distribution.raw.errors <- function(res.b,   res.h,
                                        res.b.w, res.h.w,
                                        res.b.h, res.h.h) {
  n.spp <- sum(res.b$order$N)
  res <- list(b  =res.b,   h  =res.h,
              b.w=res.b.w, h.w=res.h.w,
              b.h=res.b.h, h.h=res.h.h)
  p <- lapply(res, function(x) x$total / n.spp * 100)

  r <- range(unlist(p))
  br <- seq(r[1], r[2], length.out=30)

  h <- lapply(p, hist, br, plot=FALSE)

  xlim <- c(42, 50)
  ylim <- range(unlist(lapply(h, function(x) x$density)))

  hdf <- lapply(h, function(x) {
    foo = data.frame(c(x[c("mids","density")]))
    foo$widths = diff(x[["breaks"]])
    foo
  })

  cols <- cols.methods

  p1 <- plot_ly() %>% 
    add_bars(data = hdf$b.w, x = ~mids, y = ~density, width = ~widths,
      marker=list(color=cols[[1]]),opacity = .5) %>%
    add_bars(data = hdf$b.h, x = ~mids, y = ~density, width = ~widths,
      marker=list(color=cols[[1]]),opacity = .5) %>%
    add_bars(data = hdf$b, x = ~mids, y = ~density, width = ~widths,
      marker=list(color=cols[[1]],opacity=.7,
                  line=list(color='black',width=1))) %>%
    layout(showlegend=FALSE,
      xaxis=list(range=xlim,title=NA,showticklabels=FALSE),
      yaxis=list(range=ylim))

  p2 <- plot_ly() %>% 
    add_bars(data = hdf$h.w, x = ~mids, y = ~density, width = ~widths,
      marker=list(color=cols[[2]]),opacity = .5) %>%
    add_bars(data = hdf$h.h, x = ~mids, y = ~density, width = ~widths,
      marker=list(color=cols[[2]]),opacity = .5) %>%
    add_bars(data = hdf$h, x = ~mids, y = ~density, width = ~widths,
      marker=list(color=cols[[2]],opacity=.7,
                  line=list(color='black',width=1))) %>%
    layout(showlegend=FALSE,
      xaxis=list(range=xlim),
      yaxis=list(range=ylim))
    
  subplot(p1,p2,nrows=2) %>% # subplots need axis titles set here
    layout(xaxis2=list(title="Percentage of woody species among all vascular plants"),
           yaxis=list(title="Probability density"),
           yaxis2=list(title="Probability density"))
}
```

```r id=d586d1e8-fcb9-49d7-a444-c8bc369c9c8b
fig.distribution.raw.errors(res.b, res.h, res.b.w, res.h.w, res.b.h, res.h.h)
```

[__out-0.plotly.json][nextjournal#output#d586d1e8-fcb9-49d7-a444-c8bc369c9c8b#__out-0.plotly.json]

This draws the upper boundary of a bar chart in Plotly: x = breaks, y = heights.

```r id=f36f1c8d-4e0b-4abd-a248-41d0ec7b5fef
add_bline <- function(p, x = NULL, y = NULL, ...,
                      data = NULL, inherit = TRUE) {
    
  xvert <- c(rbind(x,x))
  yvert <- c(0,rbind(y,y),0)
  
  add_trace(p, x = ~xvert, y = ~yvert, type = "scatter", mode = "lines",
            ..., data = data, inherit = inherit)
}
```

Woodiness is structured among genera and other taxonomic groups. Make a plot of per genus/family/order estimates of woodiness:

```r id=41d716e7-6216-497b-94ec-6d838fe84dce
fig.fraction.by.group <- function(res.b, res.h, dat.g, level="genus") {
  op <- par(mfrow=c(2, 1), mar=c(2, 2, .5, .5), oma=c(2, 0, 0, 0))
  on.exit(par(op))
  lwd <- 1.5

  n.br <- c(genus=50, family=40, order=30)[[level]]
  tmp <- aggregate(dat.g[c("W", "K")], dat.g[level], sum)
  tmp <- tmp[tmp$K >= 10,] # at least 10 records per group
  h <- hist(100 * tmp$W / tmp$K, n.br, plot=FALSE)

  hdf <- data.frame(c(h[c("mids","density")]))
             
	p0 <- plot_ly() %>%
    add_bline(x = h[["breaks"]], y = h[["density"]],color=I('black'),
              showlegend=FALSE) %>%
    layout(
      xaxis=list(range=c(-2,102),zeroline=FALSE),
      yaxis=list(showticklabels=FALSE))
  
  cols <- cols.methods

  x.b <- res.b[[level]]
  x.h <- res.h[[level]]
  h.b <- hist(100*x.b$p.mean[x.b$N >= 10], n=n.br, plot=FALSE)
  h.h <- hist(100*x.h$p.mean[x.h$N >= 10], n=n.br, plot=FALSE)
  ylim <- range(h.b$density, h.h$density)
  
  p1 <- plot_ly() %>%
    add_bline(x = h.b[["breaks"]], y = h.b[["density"]],
      color=I(cols[[1]]),name="Strong prior (binomial)") %>%
    add_bline(x = h.h[["breaks"]], y = h.h[["density"]],
      color=I(cols[[2]]),name="Weak prior (hypergeometric)") %>%
    layout(xaxis=list(range=c(-2,102),zeroline=FALSE),
           yaxis=list(showticklabels=FALSE,title=NA))
  
  subplot(p0, p1, nrows = 2) %>%
    layout(legend = list(x = 0.1, y = 0.40),
           xaxis2=list(title="Percentage of woody species in genus"),
            yaxis=list(title="Probability density"),
           yaxis2=list(title="Probability density"))
}
```

```r id=69f94f73-7688-4a80-b964-ad9d0da42ada
fig.fraction.by.group(res.b, res.h, dat.g, "genus")
```

[__out-0.plotly.json][nextjournal#output#69f94f73-7688-4a80-b964-ad9d0da42ada#__out-0.plotly.json]

```r id=fc85cc1b-d812-4cec-b116-bfc2f2eb414f
fig.fraction.by.group(res.b, res.h, dat.g, "family")
```

[__out-0.plotly.json][nextjournal#output#fc85cc1b-d812-4cec-b116-bfc2f2eb414f#__out-0.plotly.json]

```r id=485e1f69-20de-4421-8d91-9e13bc555c8a
fig.fraction.by.group(res.b, res.h, dat.g, "order")
```

[__out-0.plotly.json][nextjournal#output#485e1f69-20de-4421-8d91-9e13bc555c8a#__out-0.plotly.json]

```r id=c496c6b1-0cb9-4ab2-a74a-fb523471aa7c
dat.g
```

[__out-0.csv][nextjournal#output#c496c6b1-0cb9-4ab2-a74a-fb523471aa7c#__out-0.csv]

Woodiness also varies over the tree; plot the per-order estimate of woodiness around the tree (the code to do this is not very pretty, so is kept separately in `R/plot-tree.R`)

```r id=63c79e7b-d4a8-4f34-a61a-5f9e1bd92bcc
source("R/plot-tree.R")
```

Phylogeny at the level of order:

```r id=13ef8a7c-b67a-4d86-a9ee-3e8e0660f516
phy.o <- load.phylogeny.order()
```

An option tweak is needed to get a little more height for a nice-looking plot.

```r id=055e2f80-0282-46f6-9d12-f79d4503038c
options(nextjournal.ggplot.device.default.height=6)

fig.fraction.on.phylogeny(phy.o, res.b)
```

![__out-0][nextjournal#output#055e2f80-0282-46f6-9d12-f79d4503038c#__out-0]

```r id=5c90104a-7127-41e9-848c-780d87851deb
fig.fraction.on.phylogeny(phy.o, res.h)
```

![__out-0][nextjournal#output#5c90104a-7127-41e9-848c-780d87851deb#__out-0]

Reset plot height.

```r id=8dc1918b-0604-4df5-8c0a-e82f8098398b
options(nextjournal.ggplot.device.default.height=4.471428)
```

# Survey:

Raw results of the survey:

```r id=31a53de1-8b52-4e61-a19f-35043b50c5be
d.survey <- load.survey()
```

Convert estimates to normal using logit transformation:

```r id=b8acdd3d-3851-485c-8e9b-495366ccd7f2
d.survey$Estimate.logit <- boot::logit(d.survey$Estimate / 100)
```

Model with training and familiarity as factors:

```r id=aa5fc0f3-738e-4de0-9aaf-ae9fc6baa00e
res <- lm(Estimate.logit ~ Training + Familiarity, d.survey)
summary(res)
anova(res)
```

[__out-0.csv][nextjournal#output#aa5fc0f3-738e-4de0-9aaf-ae9fc6baa00e#__out-0.csv]

Regression against |latitude|:

```r id=8bfe096b-aadb-45c3-ab99-3aa6e5546b1a
res.lat <- lm(Estimate.logit ~ abs(Lat), d.survey)
anova(res.lat)
summary(res.lat)
```

[__out-0.csv][nextjournal#output#8bfe096b-aadb-45c3-ab99-3aa6e5546b1a#__out-0.csv]

Here is the fitted result:

```r id=8e81ef83-7e75-49d5-a7b4-8a0618fea60f
s.res <- data.frame(
  fitted(res.lat),
  residuals(res.lat),
  res.lat$model["abs(Lat)"]
)
colnames(s.res) <- c("fit","res","lat")

plot_ly(s.res, type='scatter',
        x = ~lat, y = ~res, name='Data',mode='markers') %>%
  add_trace(x =~ lat, y = ~fit, name='Fit',mode='lines') %>%
  layout(xaxis=list(title="abs(Lat)"),
         yaxis=list(title="Estimate.logit"),
         legend = list(orientation = 'h'))
```

[__out-0.plotly.json][nextjournal#output#8e81ef83-7e75-49d5-a7b4-8a0618fea60f#__out-0.plotly.json]

As a categorical tropical/non-tropical variable:

```r id=1309da62-146a-4279-9f7a-d7fe8144e298
res.tro <- lm(Estimate.logit ~ Tropical, d.survey)
anova(res.tro)
summary(res.tro)
```

[__out-0.csv][nextjournal#output#1309da62-146a-4279-9f7a-d7fe8144e298#__out-0.csv]

Distribution of estimates vs different levels of botanical familiarity and education, with the estimates from the database overlaid.

```r id=0152f4e7-2e1e-4167-b562-8267e39e2426
fig.survey.results <- function(fxn, cats, d.survey, res.b, res.h, ...) {
  ci <- 100*cbind(res.b$overall.p, res.h$overall.p)
  col <- cols.methods
  
  frame = model.frame(fxn, d.survey)
  
  # For some reason add_ribbons() needs this number of data points
  # despite being on a separate x axis.
  ndat = nrow(frame) 

  plot_ly(frame) %>%
    add_boxplot(y = ~Estimate, x = cats, 
      name="Survey Responses",legendgroup="group1") %>%
    add_lines(x = ~c(0,1), y = ~c(ci["mean",1],ci["mean",1]),
      xaxis="x2",line=list(color='black',width=0.5),
      name="Binomial Mean",legendgroup="group2") %>%
    add_polygons(xaxis="x2",
      x=~c(0,0,1,1), name="Binomial Range",legendgroup="group2",
      y=~c(ci["lower",1],ci["upper",1],ci["upper",1],ci["lower",1]),
      fillcolor=col[1],line=list(color=col[1]),opacity=0.5) %>%
    add_lines(x = ~c(0,1), y = ~c(ci["mean",2],ci["mean",2]),
      xaxis="x2", line=list(color='black',width=0.5),
      name="Hypergeometric Mean",legendgroup="group3") %>%
    add_polygons(xaxis="x2",
      x=~c(0,0,1,1), name="Hypergeometric Range",legendgroup="group3",
      y=~c(ci["lower",2],ci["upper",2],ci["upper",2],ci["lower",2]),
      fillcolor=col[2],line=list(color=col[2]),opacity=0.5) %>%
    layout(
      yaxis=list(title="Estimate of percentage woodiness",
        zeroline=FALSE,showline=TRUE,range=c(0,100)),categoryorder="array",
                              categoryarray=c(1,2,3,4,5),
                              ticktext=c("1","2","3","4","5"),
      xaxis=list(title=NA,showline=TRUE),
      xaxis2=list(overlaying = "x",showticklabels=FALSE),
      ...
    )
}
```

```r id=e0997a7e-dfec-4346-ba44-4a54b6dfe419
fig.survey.results(Estimate ~ Familiarity, ~Familiarity, d.survey, res.b, res.h)

# rename training levels for display, filter NA
library(forcats)
lels <- levels(d.survey$Training)
sh_training <- fct_recode(d.survey$Training[!is.na(d.survey$Training)],
  "Postgrad"=lels[1], "Part postgrad"=lels[2], "Undergrad"=lels[3],
  "Part undergrad"=lels[4], "None"=lels[5])

fig.survey.results(Estimate ~ Training, sh_training, d.survey, res.b, res.h)
```

[__out-0.plotly.json][nextjournal#output#e0997a7e-dfec-4346-ba44-4a54b6dfe419#__out-0.plotly.json]

[__out-1.plotly.json][nextjournal#output#e0997a7e-dfec-4346-ba44-4a54b6dfe419#__out-1.plotly.json]

And the raw distribution of survey results, showing the general tendency for relatively low estimates, again overlaid with the estimates from the database:

```r id=2c21f7c8-07e3-45eb-be38-ccf7d69cf04f
fig.survey.distribution <- function(d.survey, res.b, res.h) {
  lwd <- 1.5
  col <- cols.methods

  ci <- 100*cbind(res.b$overall.p, res.h$overall.p)
  
  hist.full <- hist(d.survey$Estimate, plot=FALSE)
  hist.trop <- hist(d.survey$Estimate[d.survey$Tropical],
                    breaks=seq(0,90,10),plot=FALSE)
  hist.temp <- hist(d.survey$Estimate[!d.survey$Tropical],
                    breaks=seq(0,90,10),plot=FALSE)
  hist.data <- data.frame(c(
    hist.full[c("mids","counts")],
    hist.trop[c("mids","counts")],
    hist.temp[c("mids","counts")]
  ))
  colnames(hist.data) <- c(
    "mids","counts","mids.trop","counts.trop","mids.temp","counts.temp"
  )
  ndat = nrow(hist.data) 
  
  plot.top = 10*ceiling(max(hist.full$counts)/10)
  

  plot_ly(hist.data) %>%
    add_bars(x = ~mids, y = ~counts, type = 'bar',name="Total Counts",
      width = 10,marker = list(line = list(color = 'black',width=1)),
            legendgroup="group1") %>%
    add_polygons(name="Binomial Range",
      x = ~c(ci["lower",1],ci["lower",1],ci["upper",1],ci["upper",1]),
      y = ~c(0,plot.top,plot.top,0),legendgroup="group2",
      color=I(col[1]),opacity=0.5,line=list(width=0)) %>%
    add_lines(y = ~c(0,plot.top), x = ~c(ci["mean",1],ci["mean",1]),
      mode="lines",marker = list(opacity=0), name = "Binomial Mean",
      line=list(color='black',width=0.5),legendgroup="group2") %>%
    add_polygons(name="Hypergeometric Range",
      x = ~c(ci["lower",2],ci["lower",2],ci["upper",2],ci["upper",2]),
      y = ~c(0,plot.top,plot.top,0),legendgroup="group3",
      color=I(col[2]),opacity=0.5,line=list(width=0)) %>%
    add_lines(y = ~c(0,plot.top), x = ~c(ci["mean",2],ci["mean",2]),
      mode="lines",marker = list(opacity=0), name = "Hypergeometric Mean",
      line=list(color='black',width=0.5),legendgroup="group3") %>%
    add_bline(x = hist.trop[["breaks"]], y = hist.trop[["counts"]],
      name="Tropical Counts", color=I(cols.tropical[1]),
      legendgroup="group1") %>%
    add_bline(x = hist.temp[["breaks"]], y = hist.temp[["counts"]],
      name="Temperate Counts", color=I(cols.tropical[2]),
      legendgroup="group1") %>%
    layout(barmode='stack',
      xaxis=list(title="Estimate of percentage woodiness",
        zeroline=FALSE,showline=TRUE,range=c(0,100)),
      yaxis=list(title="Number of responses",
        showline=TRUE, range=c(0,plot.top)),
      legend = list(x = 0.65, y = 1)
    )
}
```

```r id=27d01133-edf9-46d7-b383-8a60e2e25c26
fig.survey.distribution(d.survey, res.b, res.h)
```

[__out-0.plotly.json][nextjournal#output#27d01133-edf9-46d7-b383-8a60e2e25c26#__out-0.plotly.json]

# Variability of a genus vs size

These plots look at the how variable a genus is vs its size (left column) or the number of species in a known state (right column). The top row looks at how variable the genus is (a single type or a mixed type), the middle row at the probability that a genus is varaible.  The bottom row looks at the proportion of species that are woody in a genus vs its size, testing if woody genera are relatively larg or relatively small (woody genera are relatively smaller).

```r id=841ce35e-3bb7-4afa-8456-ee1651aa1d92
fig.variability <- function(dat.g) {
  dat.g$p.rare <- (0.5 - abs(dat.g$p - 1/2)) * 2
  dat.g$variable <- dat.g$p.rare > 0

  sub <- dat.g[!is.nan(dat.g$p),]

  ## Breaks for the moving average:
  br.N <- log.seq.range(sub$N, 20)
  br.K <- log.seq.range(sub$K, 15)

  ## Classify points:
  i.N <- findInterval(sub$N, br.N, all.inside=TRUE)  
  i.K <- findInterval(sub$K, br.K, all.inside=TRUE)  

  ## Midpoints for plotting.
  mid.N <- (br.N[-1] + br.N[-length(br.N)])/2
  mid.K <- (br.K[-1] + br.K[-length(br.K)])/2

  m.N <- tapply(sub$p.rare, i.N, mean)
  m.K <- tapply(sub$p.rare, i.K, mean)
  p.N <- tapply(sub$variable, i.N, mean)
  p.K <- tapply(sub$variable, i.K, mean)
  f.N <- tapply(sub$p, i.N, mean)
  f.K <- tapply(sub$p, i.K, mean)

  pch <- 19
  cex <- 0.5
  col <- "#00000066"

  op <- par(oma=c(4.1, 4.1, 0, 0),
            mar=c(1.1, 1.1, .5, .5),
            mfrow=c(3, 2))
  
  plot(p.rare ~ N, sub, pch=pch, cex=cex, col=col, log="x",
       axes=FALSE)
  lines(m.N ~ mid.N, col="red")
  axis(1, labels=FALSE) 
  axis(2, c(0, 1), c("Single type", "50:50"), las=1)
  mtext("Variability", 2, 3)
  box(bty="l")
  label(.02, .96, 1)

  plot(p.rare ~ K, sub, pch=pch, cex=cex, col=col, log="x",
       axes=FALSE)
  lines(m.K ~ mid.K, col="red")
  axis(1, labels=FALSE) 
  axis(2, c(0, 1), labels=FALSE)
  box(bty="l")
  label(.02, .96, 2)

  plot(variable ~ N, sub, pch=pch, cex=cex, col=col, log="x",
       bty="l", las=1, axes=FALSE)
  lines(p.N ~ mid.N, col="red")
  axis(1, labels=FALSE)
  axis(2, las=1)
  mtext("Probability genus is variable", 2, 3) 
  box(bty="l") 
  label(.02, .96, 3)

  plot(variable ~ K, sub, pch=pch, cex=cex, col=col, log="x",
       bty="l", yaxt="n", las=1, axes=FALSE)
  lines(p.K ~ mid.K, col="red")
  axis(1, labels=FALSE)
  axis(2, labels=FALSE) 
  box(bty="l") 
  label(.02, .96, 4)

  plot(p ~ N, sub, pch=pch, cex=cex, col=col, log="x",
       bty="l", las=1)
  lines(f.N ~ mid.N, col="red")
  mtext("Number of species in genus", 1, 3)
  mtext("Proportion of species woody", 2, 3)
  label(.02, .96, 5)

  plot(p ~ K, sub, pch=pch, cex=cex, col=col, log="x",
       bty="l", yaxt="n", las=1)
  lines(f.K ~ mid.K, col="red")
  axis(2, labels=FALSE)
  mtext("Number of species with known state", 1, 3)
  label(.02, .96, 6)
}
```

This grid needs twice as much height to come out readable.

```r id=769d6d65-6da6-43a5-ac6f-50b13f0f55b5
options(nextjournal.ggplot.device.default.height=4.471428*2)

fig.variability(dat.g)
```

![__out-0][nextjournal#output#769d6d65-6da6-43a5-ac6f-50b13f0f55b5#__out-0]

Reset default height.

```r id=f0dae82c-ac58-4ff9-8c26-21716fa78d93
options(nextjournal.ggplot.device.default.height=4.471428)
```

# Write out partly processed data sets

To save having to rerun everything, these are the estimates by genus, family and order for the two different sampling approaches and the three different treatment of variable species (so 3 x 2 x 3 = 18 files)

*Note: locked on Nextjournal, since we don't need these with checkpoints enabled.*

```r id=e188eb12-f8c8-4d82-80f2-3cd75d288caf
write.output <- function(d, type) {
  for (tax in c("genus", "family", "order"))
    write.csv(d[[tax]],
              file.path("output/results", sprintf("%s-%s.csv", tax, type)),
              row.names=FALSE)
}
```

Core data sets:

```r id=8b745dbf-dd4a-4e55-88ff-e0092344c922
dir.create("output/results", FALSE)
write.output(res.b, "binomial-strong-prior")
write.output(res.h, "hypergeometric-weak-prior")
```

Extra data sets:

```r id=245ec53b-f1be-44b0-8f4e-9361a21963fd
write.output(res.b.w, "binomial-strong-prior-wood-biased")
write.output(res.h.w, "hypergeometric-weak-prior-wood-biased")
write.output(res.b.h, "binomial-strong-prior-herb-biased")
write.output(res.h.h, "hypergeometric-weak-prior-herb-biased")
```

Meta data for interpreting these files.

```r id=4bdb5870-5187-47a8-ab84-c8426a731e17
metadata <- 
  list(order="Taxonomic order",
       family="Taxonomic family",
       genus="Taxonomic genus",
       W="Number of species known to be woody",
       V="Number of species known to be variable",
       H="Number of species known to be herbaceous",
       N="Estimate of the number of species in the genus/family/order",
       K="Number of species with known state (W + H)",
       mean="Mean estimated number of woody species",
       lower="0.025 quantile of estimated number of woody species",
       upper="0.975 quantile of estimated number of woody species",
       p.mean="Mean estimated fration of woody species",
       p.lower="0.025 quantile of estimated fraction of woody species",
       p.upper="0.975 quantile of estimated fraction of woody species")
metadata <- data.frame(column=names(metadata),
                       description=unlist(metadata),
                       stringsAsFactors=FALSE)
write.csv(metadata, "output/results/metadata.csv", row.names=FALSE)
```

# Graphical abstract:

```r id=a222099e-ea7c-4c3f-9653-725fbd9b21f9
fig.graphical.abstract <- function(res.b, res.h, dat.g, d.survey) {
  p.raw <- sum(dat.g$W) / sum(dat.g$K)
  p.survey <- mean(d.survey$Estimate) / 100
  p.data <- mean(c(res.b$overall.p[["mean"]], res.h$overall.p[["mean"]]))
  
  f <- function(p, title) {
    pie(c(p, 1-p), c("Woody", "Herbaceous"), col=cols.woody)
    text(0, par("usr")[2], title, adj=c(0.5, 0), cex=1.2)
  }
  par(mfrow=c(1, 3), mar=rep(1, 4))
  f(p.raw,    "Raw data")
  f(p.survey, "Survey estimate")
  text(0, 1.5, "How much of the world is woody?", cex=1.5, xpd=NA)
  f(p.data,   "Bias corrected")
}
```

```r id=14063e8c-6d26-48b5-bc23-d037a03f1290
fig.graphical.abstract(res.b, res.h, dat.g, d.survey)
```

![__out-0][nextjournal#output#14063e8c-6d26-48b5-bc23-d037a03f1290#__out-0]

# Produce PDF versions of figures for publication:

```r id=d4f39b27-3337-497b-8ff8-644e73a54d12
if (interactive()) {
  to.pdf("/results/fraction-by-genus.pdf", 6, 6,
         fig.fraction.by.group(res.b, res.h, dat.g, "genus"))
  to.pdf("/results/fraction-by-family.pdf", 6, 6,
         fig.fraction.by.group(res.b, res.h, dat.g, "family"))
  to.pdf("/results/fraction-by-order.pdf", 6, 6,
         fig.fraction.by.group(res.b, res.h, dat.g, "order"))

  to.pdf("/results/fraction-on-phylogeny.pdf", 6, 6,
         fig.fraction.on.phylogeny(phy.o, res.b))

  to.pdf("/results/fraction-on-phylogeny-supp.pdf", 6, 6,
         fig.fraction.on.phylogeny(phy.o, res.h))

  to.pdf("/results/distribution-raw.pdf", 6, 4,
         fig.distribution.raw(res.b, res.h))
  to.pdf("/results/distribution-raw-errors.pdf", 6, 4,
         fig.distribution.raw.errors(res.b, res.h, res.b.w, res.h.w,
                                     res.b.h, res.h.h))

  to.pdf("/results/survey-results.pdf", 6, 4,
         fig.survey.results(d.survey, res.b, res.h))

  to.pdf("/results/survey-distribution.pdf", 6, 5,
         fig.survey.distribution(d.survey, res.b, res.h))

  to.pdf("/results/variability.pdf", 7, 8,
         fig.variability(dat.g))

  to.pdf("/results/graphical-abstract.pdf", 7, 3.5,
         fig.graphical.abstract(res.b, res.h, dat.g, d.survey))
}
```

# Version information:

```r id=df23f5a5-ba32-497d-a13f-f1f04230cf75
sessionInfo()
```

[nextjournal#output#bdbc00be-b887-4a47-ba1b-e6b69fce4280#__out-0.plotly.json]:
<https://nextjournal.com/data/QmaKq1KVZWNkgsHufp9cXc6e3PZAfavdi8zANk9Ed2TQjQ?content-type=application/json&node-id=bdbc00be-b887-4a47-ba1b-e6b69fce4280&filename=__out-0.plotly.json&node-kind=output>

[nextjournal#output#d586d1e8-fcb9-49d7-a444-c8bc369c9c8b#__out-0.plotly.json]:
<https://nextjournal.com/data/QmWi39YJzNcnHt33MoxRwNhFbp5Cm2mCy5rJLrpjge1gea?content-type=application/json&node-id=d586d1e8-fcb9-49d7-a444-c8bc369c9c8b&filename=__out-0.plotly.json&node-kind=output>

[nextjournal#output#69f94f73-7688-4a80-b964-ad9d0da42ada#__out-0.plotly.json]:
<https://nextjournal.com/data/QmZaT4jj1eYCEEiXqrtz2MZDEsqry3xcASARe6AYyPkvNJ?content-type=application/json&node-id=69f94f73-7688-4a80-b964-ad9d0da42ada&filename=__out-0.plotly.json&node-kind=output>

[nextjournal#output#fc85cc1b-d812-4cec-b116-bfc2f2eb414f#__out-0.plotly.json]:
<https://nextjournal.com/data/QmPaa2zVuVfMo4f4KR36EN4UcGjjEo4vN8TeyZtW1UDQKP?content-type=application/json&node-id=fc85cc1b-d812-4cec-b116-bfc2f2eb414f&filename=__out-0.plotly.json&node-kind=output>

[nextjournal#output#485e1f69-20de-4421-8d91-9e13bc555c8a#__out-0.plotly.json]:
<https://nextjournal.com/data/QmTDVG9N44HTLR3jMKhUf3cBWcNQtzSuQzEcCcdLgipJ2M?content-type=application/json&node-id=485e1f69-20de-4421-8d91-9e13bc555c8a&filename=__out-0.plotly.json&node-kind=output>

[nextjournal#output#c496c6b1-0cb9-4ab2-a74a-fb523471aa7c#__out-0.csv]:
<https://nextjournal.com/data/QmQgXDxD5Myg65wEiFXDh7yDNTFYbStpVuRaBi7zhqJf4j?content-type=text/csv&node-id=c496c6b1-0cb9-4ab2-a74a-fb523471aa7c&filename=__out-0.csv&node-kind=output>

[nextjournal#output#055e2f80-0282-46f6-9d12-f79d4503038c#__out-0]:
<https://nextjournal.com/data/QmddGHHWkf9Y1e2KwLQc8xkyLYQbnkKZUpMdeZ5yGbheF6?content-type=image/svg%2Bxml&node-id=055e2f80-0282-46f6-9d12-f79d4503038c&filename=__out-0&node-kind=output>

[nextjournal#output#5c90104a-7127-41e9-848c-780d87851deb#__out-0]:
<https://nextjournal.com/data/Qmb49JxJtjnnJvVux8jcLsgz8acQqmhtVKCaBhBLDn4qiT?content-type=image/svg%2Bxml&node-id=5c90104a-7127-41e9-848c-780d87851deb&filename=__out-0&node-kind=output>

[nextjournal#output#aa5fc0f3-738e-4de0-9aaf-ae9fc6baa00e#__out-0.csv]:
<https://nextjournal.com/data/QmTTNWcvcFALoGLMjkZjU338S3iszuEgedjQR6t6gbYLuR?content-type=text/csv&node-id=aa5fc0f3-738e-4de0-9aaf-ae9fc6baa00e&filename=__out-0.csv&node-kind=output>

[nextjournal#output#8bfe096b-aadb-45c3-ab99-3aa6e5546b1a#__out-0.csv]:
<https://nextjournal.com/data/QmWPBsX2S1WjWHvVYLjDEN25tQcoL1qAUNL7CvNa3d1FfG?content-type=text/csv&node-id=8bfe096b-aadb-45c3-ab99-3aa6e5546b1a&filename=__out-0.csv&node-kind=output>

[nextjournal#output#8e81ef83-7e75-49d5-a7b4-8a0618fea60f#__out-0.plotly.json]:
<https://nextjournal.com/data/QmP8xeDPS4h6XMCxUHkCx2it3R9qH7BtnQjMC5qUGTEXzc?content-type=application/json&node-id=8e81ef83-7e75-49d5-a7b4-8a0618fea60f&filename=__out-0.plotly.json&node-kind=output>

[nextjournal#output#1309da62-146a-4279-9f7a-d7fe8144e298#__out-0.csv]:
<https://nextjournal.com/data/QmaeUyWnC6eUYuqskkz21UU8hzQTY8zj5WTHpj2axF5bth?content-type=text/csv&node-id=1309da62-146a-4279-9f7a-d7fe8144e298&filename=__out-0.csv&node-kind=output>

[nextjournal#output#e0997a7e-dfec-4346-ba44-4a54b6dfe419#__out-0.plotly.json]:
<https://nextjournal.com/data/QmcmnxmLScmQ5mcaocBGNzJqHPJ6e9mXfDPGffFXJy3ibj?content-type=application/json&node-id=e0997a7e-dfec-4346-ba44-4a54b6dfe419&filename=__out-0.plotly.json&node-kind=output>

[nextjournal#output#e0997a7e-dfec-4346-ba44-4a54b6dfe419#__out-1.plotly.json]:
<https://nextjournal.com/data/QmSZWR2qkiBrCSHAtNDjVzwFBzGPCxXikZGRXBsmxGDUXa?content-type=application/json&node-id=e0997a7e-dfec-4346-ba44-4a54b6dfe419&filename=__out-1.plotly.json&node-kind=output>

[nextjournal#output#27d01133-edf9-46d7-b383-8a60e2e25c26#__out-0.plotly.json]:
<https://nextjournal.com/data/QmZUnujdXYqdpzZtSZ1v3gizt91p562W2T6ZiNifTUTLbu?content-type=application/json&node-id=27d01133-edf9-46d7-b383-8a60e2e25c26&filename=__out-0.plotly.json&node-kind=output>

[nextjournal#output#769d6d65-6da6-43a5-ac6f-50b13f0f55b5#__out-0]:
<https://nextjournal.com/data/QmdQUAX1bN64tH9L4FKtFeosdW4f2EvCfeSNNo74FxeXTT?content-type=image/svg%2Bxml&node-id=769d6d65-6da6-43a5-ac6f-50b13f0f55b5&filename=__out-0&node-kind=output>

[nextjournal#output#14063e8c-6d26-48b5-bc23-d037a03f1290#__out-0]:
<https://nextjournal.com/data/Qmc4LdAMXp47pzMFTA6xE5CMaba5KD7RLM6UNJJgZLUPqT?content-type=image/svg%2Bxml&node-id=14063e8c-6d26-48b5-bc23-d037a03f1290&filename=__out-0&node-kind=output>

<details id="com.nextjournal.article">
<summary>This notebook was exported from <a href="https://nextjournal.com/a/LJz5oibaF236dKpGXF5BA?change-id=CVLGSK86cR6hzTZfMvEEVJ">https://nextjournal.com/a/LJz5oibaF236dKpGXF5BA?change-id=CVLGSK86cR6hzTZfMvEEVJ</a></summary>

```edn nextjournal-metadata
{:article
 {:settings {:authors? true, :subtitle? true},
  :nodes
  {"0152f4e7-2e1e-4167-b562-8267e39e2426"
   {:compute-ref #uuid "9f2967fb-abc6-4ce1-8ebc-6b865d4da57b",
    :exec-duration 564,
    :id "0152f4e7-2e1e-4167-b562-8267e39e2426",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "055e2f80-0282-46f6-9d12-f79d4503038c"
   {:compute-ref #uuid "0d5286be-0a22-484f-97e8-7910fec8786b",
    :exec-duration 983,
    :id "055e2f80-0282-46f6-9d12-f79d4503038c",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "09c35830-3eae-4c8b-a561-3a7a07325e1c"
   {:compute-ref #uuid "d1fe22a6-4236-4573-a086-6b341c746cf2",
    :exec-duration 678,
    :id "09c35830-3eae-4c8b-a561-3a7a07325e1c",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "1309da62-146a-4279-9f7a-d7fe8144e298"
   {:compute-ref #uuid "3987f8e0-adf8-4f68-ad6a-b33742192f56",
    :exec-duration 1406,
    :id "1309da62-146a-4279-9f7a-d7fe8144e298",
    :kind "code",
    :output-log-lines {:stdout 20},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "13ef8a7c-b67a-4d86-a9ee-3e8e0660f516"
   {:compute-ref #uuid "a0471df7-c842-42d4-bd51-fc1fdf2ce8c9",
    :exec-duration 631,
    :id "13ef8a7c-b67a-4d86-a9ee-3e8e0660f516",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "14063e8c-6d26-48b5-bc23-d037a03f1290"
   {:compute-ref #uuid "7b40829f-d897-429d-a073-122b8aa626e9",
    :exec-duration 675,
    :id "14063e8c-6d26-48b5-bc23-d037a03f1290",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "245ec53b-f1be-44b0-8f4e-9361a21963fd"
   {:compute-ref #uuid "d8f74422-1e04-4252-95b3-d2f11aade99b",
    :exec-duration 1354,
    :id "245ec53b-f1be-44b0-8f4e-9361a21963fd",
    :kind "code",
    :locked? true,
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "27d01133-edf9-46d7-b383-8a60e2e25c26"
   {:compute-ref #uuid "cf651016-a99a-4789-b614-334117cd2cf2",
    :exec-duration 1344,
    :id "27d01133-edf9-46d7-b383-8a60e2e25c26",
    :kind "code",
    :output-log-lines {:stdout 4},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "2c21f7c8-07e3-45eb-be38-ccf7d69cf04f"
   {:compute-ref #uuid "0b0e4486-7cc9-4beb-846e-73cd8a6a8a0c",
    :exec-duration 742,
    :id "2c21f7c8-07e3-45eb-be38-ccf7d69cf04f",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "301f2420-acf7-405a-b971-dc5102143ebf"
   {:compute-ref #uuid "520283ff-fd98-44fe-a287-be7ccef4c054",
    :exec-duration 645,
    :id "301f2420-acf7-405a-b971-dc5102143ebf",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "31a53de1-8b52-4e61-a19f-35043b50c5be"
   {:compute-ref #uuid "0865faf1-a1b7-4803-9f9c-bb970ddd9485",
    :exec-duration 629,
    :id "31a53de1-8b52-4e61-a19f-35043b50c5be",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "3ebb6dbe-6da3-46d2-806a-baa7e51bb2bd"
   {:compute-ref #uuid "4fb7a30c-acfc-45b1-93cd-7c076e23beed",
    :exec-duration 746,
    :id "3ebb6dbe-6da3-46d2-806a-baa7e51bb2bd",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "41d716e7-6216-497b-94ec-6d838fe84dce"
   {:compute-ref #uuid "33787e46-ff36-4d76-9f04-ab23f555b020",
    :exec-duration 654,
    :id "41d716e7-6216-497b-94ec-6d838fe84dce",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "485e1f69-20de-4421-8d91-9e13bc555c8a"
   {:compute-ref #uuid "5bce9555-d948-491b-9e9d-c085c04b14bf",
    :exec-duration 884,
    :id "485e1f69-20de-4421-8d91-9e13bc555c8a",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "493966e3-c829-46ba-9615-58571ef44711"
   {:compute-ref #uuid "400abb4f-f57b-45bc-8241-c0a14aabef5b",
    :exec-duration 576,
    :id "493966e3-c829-46ba-9615-58571ef44711",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "4bdb5870-5187-47a8-ab84-c8426a731e17"
   {:compute-ref #uuid "6e5d8f51-ca84-4693-b691-7f98ed347554",
    :exec-duration 572,
    :id "4bdb5870-5187-47a8-ab84-c8426a731e17",
    :kind "code",
    :locked? true,
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "5c90104a-7127-41e9-848c-780d87851deb"
   {:compute-ref #uuid "66d81970-6561-472b-8536-dc1cea25b735",
    :exec-duration 922,
    :id "5c90104a-7127-41e9-848c-780d87851deb",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "63c79e7b-d4a8-4f34-a61a-5f9e1bd92bcc"
   {:compute-ref #uuid "2ffd76b8-a7ce-44e1-a4f7-8eca68660fd2",
    :exec-duration 625,
    :id "63c79e7b-d4a8-4f34-a61a-5f9e1bd92bcc",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "69f94f73-7688-4a80-b964-ad9d0da42ada"
   {:compute-ref #uuid "28264122-1a7a-4f9e-9a11-eb50b7847ec1",
    :exec-duration 1147,
    :id "69f94f73-7688-4a80-b964-ad9d0da42ada",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "6c36df22-92a4-460b-8254-ce45fb5fe6e5"
   {:environment
    [:environment
     {:article/nextjournal.id
      #uuid "02b1dc95-25ed-46d7-afec-b1ea208bcce0",
      :change/nextjournal.id
      #uuid "5d08d3f1-ce5a-4f18-9d35-87307f1046e5",
      :node/id "505a2ee8-f928-47ce-94bd-6fc78f74be3a"}],
    :id "6c36df22-92a4-460b-8254-ce45fb5fe6e5",
    :kind "runtime",
    :language "r",
    :type :nextjournal,
    :runtime/checkpoint? true},
   "769d6d65-6da6-43a5-ac6f-50b13f0f55b5"
   {:compute-ref #uuid "e8b2eb18-b309-45cc-bcb1-cdb3bb4c5478",
    :exec-duration 3807,
    :id "769d6d65-6da6-43a5-ac6f-50b13f0f55b5",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "841ce35e-3bb7-4afa-8456-ee1651aa1d92"
   {:compute-ref #uuid "d1b18951-168a-48cf-8daf-e6601e368fe0",
    :exec-duration 712,
    :id "841ce35e-3bb7-4afa-8456-ee1651aa1d92",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "88695a89-75ca-4475-8bf2-82f29530265b"
   {:compute-ref #uuid "7efbcefb-a984-4655-8b6d-47c7eef5192a",
    :exec-duration 603,
    :id "88695a89-75ca-4475-8bf2-82f29530265b",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "8b745dbf-dd4a-4e55-88ff-e0092344c922"
   {:compute-ref #uuid "d000891a-991b-4dba-b257-0036c3bd9221",
    :exec-duration 1173,
    :id "8b745dbf-dd4a-4e55-88ff-e0092344c922",
    :kind "code",
    :locked? true,
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "8bfe096b-aadb-45c3-ab99-3aa6e5546b1a"
   {:compute-ref #uuid "2274e11a-b86d-479e-b0a8-b6d156da866d",
    :exec-duration 1206,
    :id "8bfe096b-aadb-45c3-ab99-3aa6e5546b1a",
    :kind "code",
    :output-log-lines {:stdout 20},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "8dc1918b-0604-4df5-8c0a-e82f8098398b"
   {:compute-ref #uuid "d6f50901-f217-434b-8bd8-5ead754a063b",
    :exec-duration 632,
    :id "8dc1918b-0604-4df5-8c0a-e82f8098398b",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "8e81ef83-7e75-49d5-a7b4-8a0618fea60f"
   {:compute-ref #uuid "5405b5e8-e8bc-462e-84bf-3a87cf2848ec",
    :exec-duration 714,
    :id "8e81ef83-7e75-49d5-a7b4-8a0618fea60f",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "8f34d550-b5f1-4714-897b-bea40be282d1"
   {:compute-ref #uuid "0e8ca5c6-1c81-41d4-ab69-644e2cf7fb01",
    :exec-duration 582,
    :id "8f34d550-b5f1-4714-897b-bea40be282d1",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "a222099e-ea7c-4c3f-9653-725fbd9b21f9"
   {:compute-ref #uuid "d0e1dec6-b5c1-44a3-a4f4-f80c68983687",
    :exec-duration 683,
    :id "a222099e-ea7c-4c3f-9653-725fbd9b21f9",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "aa125966-03eb-4e80-ae86-e6818ebcdc17"
   {:compute-ref #uuid "cdc13b41-53ff-4a3a-8d0f-4e3bd2252a1e",
    :exec-duration 1244,
    :id "aa125966-03eb-4e80-ae86-e6818ebcdc17",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "aa5fc0f3-738e-4de0-9aaf-ae9fc6baa00e"
   {:compute-ref #uuid "77124fa9-be95-47e6-9f8b-50e30e529c57",
    :exec-duration 1364,
    :id "aa5fc0f3-738e-4de0-9aaf-ae9fc6baa00e",
    :kind "code",
    :output-log-lines {:stdout 26},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "b8acdd3d-3851-485c-8e9b-495366ccd7f2"
   {:compute-ref #uuid "a075df29-06af-40cd-9ecf-54103d32a10a",
    :exec-duration 677,
    :id "b8acdd3d-3851-485c-8e9b-495366ccd7f2",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "bdbc00be-b887-4a47-ba1b-e6b69fce4280"
   {:compute-ref #uuid "21abacb1-f8ad-46e3-b7e1-ee6425079b70",
    :exec-duration 753,
    :id "bdbc00be-b887-4a47-ba1b-e6b69fce4280",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "c496c6b1-0cb9-4ab2-a74a-fb523471aa7c"
   {:compute-ref #uuid "e20f3eff-cb0b-48c4-a1e0-c088fe67f4b8",
    :exec-duration 604,
    :id "c496c6b1-0cb9-4ab2-a74a-fb523471aa7c",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "ca8543ca-4ec0-44e8-a12f-7678c05de9bf"
   {:compute-ref #uuid "31151fba-1939-41d3-82d1-76dae67815ca",
    :exec-duration 636,
    :id "ca8543ca-4ec0-44e8-a12f-7678c05de9bf",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "d2e9b42b-efd4-43cc-b168-1b20ea6532e8"
   {:compute-ref #uuid "145293a8-f1e4-4b8e-93cf-471b5522b945",
    :exec-duration 605,
    :id "d2e9b42b-efd4-43cc-b168-1b20ea6532e8",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "d4f39b27-3337-497b-8ff8-644e73a54d12"
   {:compute-ref #uuid "5479337c-0262-4551-be0d-44a4ffdb8c1d",
    :exec-duration 700,
    :id "d4f39b27-3337-497b-8ff8-644e73a54d12",
    :kind "code",
    :locked? true,
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "d586d1e8-fcb9-49d7-a444-c8bc369c9c8b"
   {:compute-ref #uuid "c9a0fbfd-7f6c-4c59-9449-0cc51b6c4d55",
    :exec-duration 1011,
    :id "d586d1e8-fcb9-49d7-a444-c8bc369c9c8b",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "da9ce426-335a-4ba6-9e38-12d0c7d3c9a0"
   {:compute-ref #uuid "2acb7ffe-e2f7-4fd1-b914-9ffd2d36237e",
    :exec-duration 613,
    :id "da9ce426-335a-4ba6-9e38-12d0c7d3c9a0",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "de97d7f8-ad8c-4e0e-af05-766d73ac6905"
   {:compute-ref #uuid "0c9e23fc-ef6b-4ee2-a742-5d1e734b6b0a",
    :exec-duration 643,
    :id "de97d7f8-ad8c-4e0e-af05-766d73ac6905",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "df23f5a5-ba32-497d-a13f-f1f04230cf75"
   {:compute-ref #uuid "00d27fd1-91ce-4ee3-8c65-a8bbf2477606",
    :exec-duration 1467,
    :id "df23f5a5-ba32-497d-a13f-f1f04230cf75",
    :kind "code",
    :output-log-lines {:stdout 35},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "e0997a7e-dfec-4346-ba44-4a54b6dfe419"
   {:compute-ref #uuid "36aaef08-0a13-4567-814a-1e3aef45ded7",
    :exec-duration 936,
    :id "e0997a7e-dfec-4346-ba44-4a54b6dfe419",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "e188eb12-f8c8-4d82-80f2-3cd75d288caf"
   {:compute-ref #uuid "647352d4-c8cc-404b-ae5a-db35c3002016",
    :exec-duration 917,
    :id "e188eb12-f8c8-4d82-80f2-3cd75d288caf",
    :kind "code",
    :locked? true,
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "f0dae82c-ac58-4ff9-8c26-21716fa78d93"
   {:compute-ref #uuid "4f66f852-5e62-4ce4-b1e6-568e999c4b18",
    :exec-duration 632,
    :id "f0dae82c-ac58-4ff9-8c26-21716fa78d93",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "f28dc23d-1cda-4e4e-966c-9d2cabef57cf"
   {:compute-ref #uuid "18038d24-3973-4460-94e6-5b803cb29193",
    :exec-duration 585,
    :id "f28dc23d-1cda-4e4e-966c-9d2cabef57cf",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "f36f1c8d-4e0b-4abd-a248-41d0ec7b5fef"
   {:compute-ref #uuid "89b576f9-2a53-4abf-861e-769db3baa91e",
    :exec-duration 567,
    :id "f36f1c8d-4e0b-4abd-a248-41d0ec7b5fef",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]},
   "fc85cc1b-d812-4cec-b116-bfc2f2eb414f"
   {:compute-ref #uuid "013f342c-43f2-4a6d-aadc-88f1f96dca0d",
    :exec-duration 885,
    :id "fc85cc1b-d812-4cec-b116-bfc2f2eb414f",
    :kind "code",
    :output-log-lines {:stdout 0},
    :runtime [:runtime "6c36df22-92a4-460b-8254-ce45fb5fe6e5"]}},
  :nextjournal/id #uuid "02b236fb-c35b-4849-95d2-4b21ec344a0d",
  :article/change
  {:nextjournal/id #uuid "5d093958-4ec3-4069-b0c4-39c660c27395"}}}

```
</details>