#### Modeling in the tidyverse, an example predicting tree biomass

##### March 24, 2020

`statistics`

`Regression`

`tidyverse`

`broom`

`modelr`

One of the most revolutionary advances in the last 10 years in the R program has been the tidyverse. I recently watched Hadley Wickham’s talk on the tidyverse from rstudio::conf 2020 and it was a great synopsis of the current status and future of the packages.

The growth of the tidyverse packages is astounding. As an example, the cumulative number of downloads for the `dplyr`

package has reached over 10 million. Incredible!

The tidyverse offers a suite of packages and functions to import, tidy, transform, visualize, and model data. I’ve known well how to use the data wrangling functions from packages like `dplyr`

and data visualization from `ggplot2`

. It’s been that last core function of the tidyverse, modeling data, that I’ve never really realized.

The `modelr`

and `broom`

packages handle modeling in the tidyverse. In this post, I’ll also use the `knitr`

and `kableExtra`

to design some of the tables:

```
library(tidyverse)
library(modelr)
library(knitr)
library(kableExtra)
```

## Tree biomass data

My goal in this post is to fit a series of models that determine the aboveground biomass of trees using tree diameter as a predictor variable. I’ve gathered data from LegacyTreeData, an online repository of individual tree measurements such as volume, weight, and wood density. I queried the database to provide all tree measurements from the US State of Georgia.

After the query, there are 608 observations from seven species that contain a value for the tree’s diameter at breast height, measured in inches (`ST_OB_D_BH`

) and its aboveground dry weight in pounds (`AG_DW`

).

Here’s a “hexagonal heat map” produced in `ggplot`

that divides the x- and y-axes into 20 hexagons. The color of each hexagon reflects the number of observations in each area of the figure. In this data set, most trees are small in diameter and do not weigh a lot:

```
ggplot(tree, aes(ST_OB_D_BH, AG_DW)) +
geom_hex(bins = 20) +
labs(x = "Diameter at breast height (inches)",
y = "Aboveground dry weight (pounds)") +
theme(panel.background = element_rect(fill = "NA"),
axis.line = element_line(color = "black"))
```

The data contain relatively small trees, with DBH ranging values from 0.7 to 8.5 inches. There are at least 40 observations for seven primary species common to Georgia. These are mostly different kinds of pine trees in addition to one hardwood tree (sweetgum):

Species | Num trees | Mean DBH | Max DBH | Min DBH | Mean weight | Max weight | Min weight |
---|---|---|---|---|---|---|---|

Loblolly pine | 186 | 3.1 | 8.5 | 0.7 | 37.1 | 422.8 | 1.7 |

Shortleaf pine | 100 | 3.0 | 4.9 | 1.1 | 26.3 | 78.5 | 1.7 |

Longleaf pine | 80 | 3.0 | 4.9 | 1.2 | 36.5 | 119.1 | 1.9 |

Slash pine | 80 | 3.0 | 4.9 | 1.2 | 30.6 | 105.8 | 1.9 |

Virginia pine | 80 | 3.0 | 4.9 | 1.0 | 37.7 | 128.1 | 2.1 |

Sweetgum | 42 | 1.8 | 3.3 | 1.0 | 7.8 | 26.5 | 1.5 |

Eastern white pine | 40 | 3.0 | 4.9 | 1.0 | 27.5 | 67.7 | 2.1 |

## Models of aboveground tree biomass

To start our modeling analysis, we can fit a model predicting aboveground dry weight using tree diameter. From above, we can see a clear nonlinear trend in the data, indicating that we may want to log-transform the variables. We can do that with the `mutate`

function.

Then, we can start by fitting a linear model using `lm()`

to all observations in the data set. Printing the summary of `tree_mod`

indicates that the model fits the data well, with small standard errors of the coefficients and an R-squared value of 0.9251:

```
tree <- tree %>%
mutate(l_ST_OB_D_BH = log2(ST_OB_D_BH),
l_AG_DW = log2(AG_DW))
tree_mod <- lm(l_AG_DW ~ l_ST_OB_D_BH, data = tree)
summary(tree_mod)
```

```
##
## Call:
## lm(formula = l_AG_DW ~ l_ST_OB_D_BH, data = tree)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21620 -0.27276 -0.03442 0.22539 2.60255
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.02268 0.04085 25.04 <2e-16 ***
## l_ST_OB_D_BH 2.26555 0.02619 86.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4428 on 606 degrees of freedom
## Multiple R-squared: 0.9251, Adjusted R-squared: 0.925
## F-statistic: 7485 on 1 and 606 DF, p-value: < 2.2e-16
```

Chapter 24 in the popular book R for Data Science discusses making a grid of data to investigate model predictions. In this step, I’ll use the `data_grid()`

function to generate a grid of data points. Then, I’ll use the `add_predictions()`

function to add the model predictions from `tree_mod`

to complete our data grid. Finally, I’ll back-transform the model predictions to get them into units we understand well.

The model appears to perform well, albeit underpredicting the biomass of a few larger-diameter loblolly pine trees:

```
grid <- tree %>%
data_grid(ST_OB_D_BH = seq_range(ST_OB_D_BH, 20)) %>%
mutate(l_ST_OB_D_BH = log2(ST_OB_D_BH)) %>%
add_predictions(tree_mod, "l_AG_DW") %>%
mutate(AG_DW = 2 ^ l_AG_DW)
ggplot(tree, aes(ST_OB_D_BH, AG_DW)) +
geom_hex(bins = 20) +
geom_line(data = grid, color = "red", size = 1) +
labs(x = "Diameter at breast height (inches)",
y = "Aboveground dry weight (pounds)")
```

A data analyst’s best friend in the tidyverse is the `group_by`

statement. We can fit the same model as we did earlier, but this time I’ll specify it for each of the seven species using `group_by`

. The `tidy()`

function available in the broom package provides a set of functions that put model output into data frames.

Here, we can see that the species have a different set of coefficients and other attributes like p-values:

```
tree_coef<- tree %>%
group_by(Species) %>%
do(broom::tidy(lm(l_AG_DW ~ l_ST_OB_D_BH, .)))
tree_coef
```

```
## # A tibble: 14 x 6
## # Groups: Species [7]
## Species term estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Eastern white pine (Intercept) 1.02 0.109 9.40 1.88e-11
## 2 Eastern white pine l_ST_OB_D_BH 2.21 0.0693 31.9 4.86e-29
## 3 Loblolly pine (Intercept) 1.53 0.0775 19.8 2.45e-47
## 4 Loblolly pine l_ST_OB_D_BH 1.96 0.0475 41.3 5.82e-95
## 5 Longleaf pine (Intercept) 0.445 0.103 4.33 4.36e- 5
## 6 Longleaf pine l_ST_OB_D_BH 2.73 0.0652 41.8 3.46e-55
## 7 Shortleaf pine (Intercept) 0.600 0.0683 8.80 4.95e-14
## 8 Shortleaf pine l_ST_OB_D_BH 2.40 0.0432 55.4 8.62e-76
## 9 Slash pine (Intercept) 0.532 0.0962 5.54 4.00e- 7
## 10 Slash pine l_ST_OB_D_BH 2.54 0.0610 41.6 5.33e-55
## 11 Sweetgum (Intercept) 0.749 0.0740 10.1 1.34e-12
## 12 Sweetgum l_ST_OB_D_BH 2.38 0.0834 28.5 3.42e-28
## 13 Virginia pine (Intercept) 1.12 0.0766 14.6 5.04e-24
## 14 Virginia pine l_ST_OB_D_BH 2.38 0.0480 49.6 9.77e-61
```

A way to visualize the species differences is to plot the intercept and slope coefficients with standard errors. Here we can see that all errors bars do not overlap with zero, indicating they’re good models:

```
ggplot(tree_coef, aes(estimate, 1)) +
geom_point() +
geom_errorbarh(aes(xmin = estimate - (2*std.error),
xmax = estimate + (2*std.error),
height = 0.5)) +
scale_y_continuous(limits = c(0,2))+
facet_grid(term~Species)+
labs(x = "Coefficient", y = " ")+
ggtitle("Coefficients for determining aboveground dry weight \n(+/- two standard errors) for seven tree species")+
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank())
```

## Analysis of model predictions

Aside from coefficients, we might be interested in species-specific predictions from a model. The `nest()`

function creates a list of data frames containing all the nested variables in an object. I think of a nested data frame as a “data frame containing many data frames”.

The `by_spp`

object will store this list of data frames for each species so that we can work with them:

```
by_spp <- tree %>%
group_by(Species) %>%
nest()
species_model <- function(df){
lm(l_AG_DW ~ l_ST_OB_D_BH, data = df)
}
models <- map(by_spp$data, species_model)
by_spp <- by_spp %>%
mutate(model = map(data, species_model))
by_spp
```

```
## # A tibble: 7 x 3
## # Groups: Species [7]
## Species data model
## <chr> <list<df[,7]>> <list>
## 1 Sweetgum [42 x 7] <lm>
## 2 Loblolly pine [186 x 7] <lm>
## 3 Shortleaf pine [100 x 7] <lm>
## 4 Virginia pine [80 x 7] <lm>
## 5 Eastern white pine [40 x 7] <lm>
## 6 Slash pine [80 x 7] <lm>
## 7 Longleaf pine [80 x 7] <lm>
```

Similar to what we did above to the all-species equation, we can map the model predictions to the nested object, adding another variable called `preds`

:

```
by_spp <- by_spp %>%
mutate(preds = map2(data, model, add_predictions))
by_spp
```

```
## # A tibble: 7 x 4
## # Groups: Species [7]
## Species data model preds
## <chr> <list<df[,7]>> <list> <list>
## 1 Sweetgum [42 x 7] <lm> <tibble [42 x 8]>
## 2 Loblolly pine [186 x 7] <lm> <tibble [186 x 8]>
## 3 Shortleaf pine [100 x 7] <lm> <tibble [100 x 8]>
## 4 Virginia pine [80 x 7] <lm> <tibble [80 x 8]>
## 5 Eastern white pine [40 x 7] <lm> <tibble [40 x 8]>
## 6 Slash pine [80 x 7] <lm> <tibble [80 x 8]>
## 7 Longleaf pine [80 x 7] <lm> <tibble [80 x 8]>
```

Then, we can unnest the model predictions. Unnesting is the opposite of what we’ve done in the previous step. This time we’re taking the nested data frame and turning it into a “regular” data frame.

The model predictions can now be visualized by species to better understand differences in aboveground biomass:

```
preds <- unnest(by_spp, preds)
preds %>%
mutate(tpred = 2 ^ pred) %>%
ggplot(aes(ST_OB_D_BH, tpred)) +
geom_line(aes(group = Species), size = 1) +
labs(x = "Diameter at breast height (inches)",
y = "Predicted aboveground\ndry weight (pounds)")+
facet_wrap(~ Species)
```

## Analysis of model residuals

Any good data analysis involving modeling will also analyze the residuals. Just like we added model predictions with the `add_predictions`

function and nesting/unnesting, we can add residuals with the `add_residuals`

statement.

We see that the model residuals are generally centered around zero, with those four bigger loblolly pine trees still giving us some issues:

```
by_spp <- by_spp %>%
mutate(resids = map2(data, model, add_residuals))
by_spp
```

```
## # A tibble: 7 x 5
## # Groups: Species [7]
## Species data model preds resids
## <chr> <list<df[,7]>> <list> <list> <list>
## 1 Sweetgum [42 x 7] <lm> <tibble [42 x 8~ <tibble [42 x 8~
## 2 Loblolly pine [186 x 7] <lm> <tibble [186 x ~ <tibble [186 x ~
## 3 Shortleaf pine [100 x 7] <lm> <tibble [100 x ~ <tibble [100 x ~
## 4 Virginia pine [80 x 7] <lm> <tibble [80 x 8~ <tibble [80 x 8~
## 5 Eastern white pi~ [40 x 7] <lm> <tibble [40 x 8~ <tibble [40 x 8~
## 6 Slash pine [80 x 7] <lm> <tibble [80 x 8~ <tibble [80 x 8~
## 7 Longleaf pine [80 x 7] <lm> <tibble [80 x 8~ <tibble [80 x 8~
```

```
resids <- unnest(by_spp, resids)
resids %>%
ggplot(aes(ST_OB_D_BH, resid)) +
geom_point(aes(group = Species), alpha = 1/3) +
geom_smooth(se = F) +
labs(x = "Diameter at breast height (inches)",
y = "Residual")+
geom_abline(intercept = 0, slope = 0, color = "pink")+
facet_wrap(~Species)
```

## Assessing model quality

We also might be interested in other model attributes such as the R-squared, log-likelihood, AIC, and BIC values. These indicators of model fit can be useful to compare with different model forms. The `glance()`

function from the `broom`

package allows you to obtain a data frame with a single row.

```
glance <- by_spp %>%
mutate(glance = map(model, broom::glance)) %>%
unnest(glance, .drop = T)
glance %>%
arrange(desc(r.squared))
```

```
## # A tibble: 7 x 16
## # Groups: Species [7]
## Species data model preds resids r.squared adj.r.squared sigma
## <chr> <list<df> <lis> <lis> <list> <dbl> <dbl> <dbl>
## 1 Virgin~ [80 x 7] <lm> <tib~ <tibb~ 0.969 0.969 0.268
## 2 Shortl~ [100 x 7] <lm> <tib~ <tibb~ 0.969 0.969 0.266
## 3 Easter~ [40 x 7] <lm> <tib~ <tibb~ 0.964 0.963 0.288
## 4 Longle~ [80 x 7] <lm> <tib~ <tibb~ 0.957 0.957 0.350
## 5 Slash ~ [80 x 7] <lm> <tib~ <tibb~ 0.957 0.956 0.327
## 6 Sweetg~ [42 x 7] <lm> <tib~ <tibb~ 0.953 0.952 0.245
## 7 Loblol~ [186 x 7] <lm> <tib~ <tibb~ 0.903 0.902 0.509
## # ... with 8 more variables: statistic <dbl>, p.value <dbl>, df <int>,
## # logLik <dbl>, AIC <dbl>, BIC <dbl>, deviance <dbl>, df.residual <int>
```

We can visualize the R-squared values from the seven species and see that Virginia pine and shortleaf pine show the highest R-squared, and loblolly pine the lowest:

```
glance %>%
ggplot(aes(reorder(Species, r.squared), r.squared)) +
geom_bar(stat = "identity")+
coord_flip()+
geom_smooth(se = F) +
labs(x = " ",
y = "R-squared")+
theme(legend.position = "none")
```

## Conclusion

Like other data analysts, I’ve spent a lot of time with the tidyverse doing data wrangling and data visualization. But don’t underestimate the power of the modeling functions available in the tidyverse. The suite of functions available in the `modelr`

and `broom`

packages make modeling easier, particularly if the same model forms need to be fit iteratively to different levels (in our case, tree species).

*By Matt Russell. Leave a comment below or email Matt with any questions or comments.*