Today we are talking about structural equation models (SEM). There are lots of synonyms, sub-categories, and adjacent techniques that you may have heard before–covariance structure analysis, linear structural relations, path analysis, latent variable modeling, causal modeling, confirmatory factor analysis, exploratory factor analysis, latent growth modeling, mixture/multigroup/heirarchical/multilevel structural modeling, construct analysis, etc., etc…

So many names exist because there are a LOT of things that you can do with these types of models. It turns out that allowing for and testing structure in modeling can help solve lots of problems in research.

I think that the best way to understand SEMs is to start with simple regressions. Let’s consider a regression predicting a rooster crowing from the rising of the sun:

y = b₀ + b₁x + ε

y ~ Rooster Crowing,
x ~ Sun Rising,
b₀ = 0,
ε ~ N(0, 1)

This is a silly regression because we all understand the relationship: the rooster goes cock-a-doodle-do when the sun crests the horizon. But in the language of mathematics, there is no reason I can’t rewrite this equation as:

x = ( y – ε ) b₁^-1

This formulation make no sense. Basically, we are saying that the rising of the sun is a function of the rooster crowing! Even though this is totally mathematically viable, it defies our common sense of causation.

The language of structural equation modeling allows one way to impose some directional structure on mathematical equations. We usually visualize that language as directed graphical models like the one below.

In a graphical model, the observed variables are displayed as boxes Unobserved or latent variables are displayed as circles. Constants are triangles. The functional relationship between variables is displayed as directional arrows. Non-directional or double-headed arrows indicate a variance or covariance.

This graphical model above is the same model as the regression in the equation above. Our independent variable x has a linear relationship with the dependent variable y with the slope parameter b₁. y has a constant intercept b₀ of 0. Finally, the residual variation in y not caused by x is assumed to come from some other unobserved cause with it’s own variance. Rather than thinking of variables as independent or dependent, we used the terms exogenous or endogenous. Exogenous variables are those (like x) that have no incoming paths. These are the most ‘upstream’ variables in the causal paths. Endogenous variables are those that receive causal paths. We’ll see later that some endogenous variables can also be causes of other endogenous variables, but they are still considered endogenous.

In practice, we generally ignore the peripherals and intercepts in structural models, yielding a simplified graph:

Now that we have a shared language. Let’s take a look at a toy example to understand why SEMs can be so useful in research.

Motivating example

You will need a handful of packages for this tutorial:

packages_for_sem_workshop <-
  c(
    'tidyverse', # basic data wrangling
    'tidygraph', # graph visualization
    'ggraph', # graph visualization
    'lavaan', # sem tools
    'piecewiseSEM', # sem tools
    'mgcv', # nonlinear modeling
    'lme4', # random effect modeling
    'cvsem' # cross-validating sems
    )

install_and_load_packages <-
  function(x){
    for( i in x ){
      if( ! require( i , character.only = TRUE ) ){
        install.packages( i , dependencies = TRUE )
        library( i , character.only = TRUE )
      }
    }
  }

install_and_load_packages(packages_for_sem_workshop)

Fitting a simple regression

I find that, for most applied researchers, the language of code is more intuitive than the language of math. So, let’s simulate a toy dataset and see how we can fit the same system as a linear model or an SEM.

simple_ex <-
  data.frame(
    x = runif(n = 100, min = 0, max = 10),
    e = rnorm(n = 100, 0, 1)
  )

simple_ex <- simple_ex %>%
  mutate(
    y = 1 + 2.5*x + e
  )

Now, let’s fit a good old-fashioned linear regression model:

fit_simple_ex_lm <- 
  lm(y ~ x, data = simple_ex)

summary(fit_simple_ex_lm)

> summary(fit_simple_ex_lm)

Call:
lm(formula = y ~ x, data = simple_ex)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.61653 -0.49110 -0.01622  0.51680  2.76976 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.08641    0.17916   6.064 2.49e-08 ***
x            2.52441    0.02871  87.913  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8972 on 98 degrees of freedom
Multiple R-squared:  0.9875,	Adjusted R-squared:  0.9874 
F-statistic:  7729 on 1 and 98 DF,  p-value: < 2.2e-16

We can fit the same model using the lavaan package. The syntax for lavaan is to pass the list of models as a simple character string separated by next lines.

simple_ex_sem <-
  '
  y ~ x
  y ~ 1
'

fit_simple_ex_sem <-
  sem(model = simple_ex_sem,
      data = simple_ex)

summary(fit_simple_ex_sem)

> summary(fit_simple_ex_sem)
lavaan 0.6-18 ended normally after 1 iteration

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         3

  Number of observations                           100

Model Test User Model:
                                                      
  Test statistic                                 0.000
  Degrees of freedom                                 0

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  y ~                                                 
    x                 2.524    0.028   88.805    0.000

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .y                 1.086    0.177    6.125    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .y                 0.789    0.112    7.071    0.000

Lavaan give us an informative output. We see that the parameters are estimated with ‘ML’, maximum likelihood, which differs from ‘lm’ which estimates via ordinary least squares (OLS). This means that estimates might differ slightly. We will trust the default optimizer ‘NLMINB‘.

We see that the model is estimating 3 parameters–the regression slope, the intercept, and the residual variance. Note that we typically do not estimate the intercepts in SEMs (I’m doing so here for continuity). To fit the SEM without the intercept, simply remove the 'y ~ 1' from the model list.

We will pass over the Model Tests and focus on the parameter estimates. The slope estimate of 2.524 and the intercept estimate of 1.086 perfectly match our lm estimates. We don’t get an estimate of the residual variance from lm, but we can extract the residuals from the lm model and calculate the variance ourselves.

fit_simple_ex_lm %>%
  resid() %>%
  var()

> fit_simple_ex_lm %>%
+   resid() %>%
+   var()
[1] 0.7968917

That values is a bit off due to the difference in estimation techniques, I believe.

Fitting complex systems

Now that we’ve seen how linear regression models can be fit as a special case of SEM, let’s take a look at an example that shows where SEM surpasses linear regression.

Again, we will use simulated data, but we’ll download the data so that we can’t see the underlying structure used to generate it. Instead, we’ll use SEMs to determine the structure.

source("https://raw.githubusercontent.com/andisa01/202407_SEM_turorial/main/scripts/SEM_tutorial_example_source.R")

The traditional approach

First, let’s analyze this data using linear regression modeling with stepwise variable selection like almost (the technique use by almost half of all ecological and animal behavioral researcher (including my past self)). In this methodology, we start with a ‘full’ model including all variables and their interactions. Then, we drop non-significant parameters starting with the least parsimonious and refit the model until only significant parameters remain.

# Fit the full or 'global model'
mod_ex01_full <- lm(Y ~ X1 + X2 + X1:X2, data = example01_data_anon)
summary(mod_ex01_full)

> summary(mod_ex01_full)

Call:
lm(formula = Y ~ X1 + X2 + X1:X2, data = example01_data_anon)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.5935  -3.7755   0.1861   3.6929  13.7593 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.163110   7.944663  -0.146    0.884
X1          -0.264413   0.387296  -0.683    0.496
X2           0.155742   0.109397   1.424    0.156
X1:X2        0.002022   0.003950   0.512    0.609

Residual standard error: 5.066 on 196 degrees of freedom
Multiple R-squared:  0.144,	Adjusted R-squared:  0.1309 
F-statistic: 10.99 on 3 and 196 DF,  p-value: 1.06e-06

None of our parameters are significant, so we drop the interaction term.

# Drop the interaction
lm(Y ~ X1 + X2, data = example01_data_anon) %>% 
  summary()

> lm(Y ~ X1 + X2, data = example01_data_anon) %>% 
+   summary()

Call:
lm(formula = Y ~ X1 + X2, data = example01_data_anon)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.5962  -3.8684   0.1564   3.6213  13.8291 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -4.73378    3.79270  -1.248  0.21346   
X1          -0.08449    0.16216  -0.521  0.60295   
X2           0.19734    0.07308   2.701  0.00753 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.057 on 197 degrees of freedom
Multiple R-squared:  0.1428,	Adjusted R-squared:  0.1341 
F-statistic: 16.41 on 2 and 197 DF,  p-value: 2.555e-07

Now, x₂  is significant, but x₁ is not. So, we drop x₁.

lm(Y ~ X2, data = example01_data_anon) %>% 
  summary()

> lm(Y ~ X2, data = example01_data_anon) %>% 
+   summary() # Drop X1

Call:
lm(formula = Y ~ X2, data = example01_data_anon)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.4915  -3.8082  -0.0172   3.6272  13.5921 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.28441    2.57317  -1.276    0.203    
X2           0.16227    0.02839   5.716 3.97e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.047 on 198 degrees of freedom
Multiple R-squared:  0.1416,	Adjusted R-squared:  0.1373 
F-statistic: 32.67 on 1 and 198 DF,  p-value: 3.966e-08

And, now we have a minimum adequate model to explain the system! We can explain 14% of the variation of y with x₂. For every 1 unit increase in x₂, we expect y to increase by 0.16. And look at that p-value! Let’s go publish a paper!!!

In fact, this method of stepwise variable selection is so common that you can do it all in one command with dredge() from the MuMIn package.

MuMIn::dredge(mod_ex01_full)

> MuMIn::dredge(mod_ex01_full)
Fixed term is "(Intercept)"
Global model call: lm(formula = Y ~ X1 + X2 + X1:X2, data = example01_data_anon)
---
Model selection table 
   (Int)       X1     X2    X1:X2 df   logLik   AICc delta weight
3 -3.284          0.1623           3 -606.551 1219.2  0.00  0.626
4 -4.734 -0.08449 0.1973           4 -606.413 1221.0  1.81  0.254
8 -1.163 -0.26440 0.1557 0.002022  5 -606.279 1222.9  3.64  0.101
2  4.870  0.31890                  3 -610.048 1226.2  6.99  0.019
1 11.280                           2 -621.824 1247.7 28.49  0.000
Models ranked by AICc(x) 

The complexity of complex systems

The problem with the linear regression approach is that there is an inherent assumption that the independent variables are, well… independent. But, in natural systems, there are almost always relationships between the independent variables. For instance, here are just a few structures that could underly the relationship between y, x₁, and x₂.

The model structure we derived in our stepwise regression model is the single effect graph in the top left. The model implied by the multiple regression model is a ‘common effect’ structure where the exogenous variables are uncorrelated–in other words x₁ and x₂ have independent effects on y.

But other structures could exist that we cannot easily capture in regression. For instance, For instance,  x₁ and x₂ have independent effects on y but remain correlated with each other (common cause with correlation (bottom right)). Or x₁ may have no direct effect on y, but may affect x₂ which affects y in turn (this is called a chain or fully mediated effect (top middle)). Or, x₁ might directly affect y in addition to the indirect effect (partial mediator (top right)).

x₂ may not have any effect on y at all, but may still covary because both are directly affected by x₁ (common cause (bottom left)).

Using SEMs to compare structural hypotheses

Given all the possible structures, how can we ever know which governs our study system? Well, the first pass is to use common sense. If x₁ is height and x₂ is weight, the principle of allometry should have us exclude any model without a relationship between them.

The second pass is to use knowledge of your system from the literature. For example, if x₁ is maternal genotype and x₂ is F1 genotype (of the progeny), there can be no direct effect of the maternal effect on y if the species has an annual lifecycle, but might be partially mediated for perennial species.

Once you have a set of plausible structural hypotheses, we can use the mechanics of SEM to ask which structure best fits the data.

For now, we will assume that all 6 hypotheses above are plausible. We’ll fit each in turn. To do so, I’ll introduce a new operator ~~ which relates the undirected covariation between two variables (i.e. covariance) or with itself (i.e. variance).

# single effect
ex01_formula_x2effect <- '
Y ~ X2
X1 ~~ X1
'

ex01_sem_x2effect <- sem(ex01_formula_x2effect, data = example01_data_anon)

summary(ex01_sem_x2effect)

> summary(ex01_sem_x2effect)
lavaan 0.6-18 ended normally after 8 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         3

  Number of observations                           200

Model Test User Model:
                                                      
  Test statistic                               377.733
  Degrees of freedom                                 2
  P-value (Chi-square)                           0.000

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  Y ~                                                 
    X2                0.162    0.028    5.745    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
    X1               32.093    3.209   10.000    0.000
   .Y                25.220    2.522   10.000    0.000

We won’t use the summary output for now, so I will exclude it when fitting the rest of the models.

# chain
ex01_formula_chain <- '
X2 ~ X1
Y ~ X2
'

ex01_sem_chain <- sem(ex01_formula_chain, data = example01_data_anon)

# partially mediated
ex01_formula_mediator <- '
X2 ~ X1
Y ~ X2
Y ~ X1
'

ex01_sem_mediator <- sem(ex01_formula_mediator, data = example01_data_anon)

# common cause
ex01_formula_commoncause <- '
X2 ~ X1
Y ~ X1
'

ex01_sem_commoncause <- sem(ex01_formula_commoncause, data = example01_data_anon)

# common effect (uncorrelated)
ex01_formula_commoneffect <- '
Y ~ X1
Y ~ X2
'

ex01_sem_commoneffect <- sem(ex01_formula_commoneffect, data = example01_data_anon)

# common effect (correlated)
ex01_formula_commoneffect2 <- '
Y ~ X1
Y ~ X2
X1 ~~ X2
'

ex01_sem_commoneffect2 <- sem(ex01_formula_commoneffect2, data = example01_data_anon)

We’ll use the anova() command to get some summary statistics on all of the models.

anova(
  ex01_sem_chain, 
  ex01_sem_commoneffect, 
  ex01_sem_commoncause, 
  ex01_sem_x2effect, 
  ex01_sem_commoneffect2, 
  ex01_sem_mediator
)

Chi-Squared Difference Test

                       Df    AIC    BIC    Chisq Chisq diff RMSEA Df diff Pr(>Chisq)    
ex01_sem_commoneffect   0 1218.8 1228.7   0.0000                                        
ex01_sem_commoncause    0 2425.5 2442.0   0.0000       0.00 0.000       0               
ex01_sem_commoneffect2  0 3688.8 3708.6   0.0000       0.00 0.000       0               
ex01_sem_mediator       0 2425.5 2442.0   0.0000       0.00 0.000       0               
ex01_sem_chain          1 2423.8 2437.0   0.2754       0.28 0.000       1     0.5997    
ex01_sem_x2effect       2 2480.4 2490.3 377.7330     377.46 1.372       1     <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Warning message: lavaan->lavTestLRT():  
   some models have the same degrees of freedom 

The models are not mutually nested, so we can’t interpret the chi-squared test. However, we can use the two information criteria (AIC and BIC) since those are agnostic to the particular models in comparison. Both AIC and BIC indicate that the uncorrelated common effect structure is the most likely model.

That’s not the end of the story, though…

Cross-validating SEMs

As an applied data scientist, I’ve become very skeptical of parametric goodness-of-fit metrics. After all, I can fit a neural net that perfectly fits any dataset by memorizing it. What I am really interested is in how well a model can make accurate predictions.

Ideally, we would fit our model on one dataset and then collect separate data to validate its predictive accuracy. This is called cross-validation. Given the limitations of field research, however, a more common approach is to split your single dataset into two groups. One for fitting and one for validation.

The simple form of cross-validation only gives us one chance to measure accuracy (or lack thereof, i.e. error). In order to maximize the use of our data and get better estimates, it is a standard practice to split the dataset into a given number of evenly distributed, randomly sampled splits (i.e. folds). Then we can perform the cross-validation as many times as we have folds. This is called k-fold cross-validation.

We can use the library cvsem to easily conduct a 10-fold cross validation on all of our models.

# Cross validating
models_to_cv <- 
  cvgather(
    ex01_sem_chain, 
    ex01_sem_commoneffect, 
    ex01_sem_commoncause, 
    ex01_sem_x2effect, 
    ex01_sem_commoneffect2, 
    ex01_sem_mediator
    )

cvsem(
  data = example01_data_anon,
  Models = models_to_cv,
  k = 10
)

Cross-Validation Results of 6 models 
based on  k =  10 folds. 

                   Model E(KL-D)   SE
1         ex01_sem_chain    0.32 0.16
3   ex01_sem_commoncause    0.33 0.16
6      ex01_sem_mediator    0.33 0.16
4      ex01_sem_x2effect    5.90 1.32
5 ex01_sem_commoneffect2    7.16 1.23
2  ex01_sem_commoneffect    7.20 1.23

Both of the ‘common cause’ structures exhibit the worst predictive performance. They have the highest errors (by default, the error is calculated via KL-Divergence–the difference between the covariance matrix of the validation set and the covariance matric implied by the fitted model). To be fair, a very quick look at the relationship between x₁ and x₂ could have told us that any model that fails to account for the relationship of these two models is not realistic.

The next highest error comes from the model suggested by our initial stepwise regression method. Thus, based on predictive performance, we can conclusively exclude these hypothesized structures.

The chain, common cause, and mediator structures all have similar predictive accuracy. How can we decide which is the true structure of our system?  Unfortunately, there is no statistical method to differentiate with our data. However, we can still use the structural model to help us develop further experiments to select the final model. For example, we could envision conducting an experiment where we held x₂ constant in one set of plots and let it vary in another.

If the system were a common cause structure, holding x₂ would cause no difference in the y values between plots. However, if the system were a fully mediated chain through x₂ , holding x₂ would completely decouple the association between x₁ and y, whereas controlling x₂ in a partially mediated system would only attenuate the covariation between x₁ and y.

In our toy example, the true causal structure will be easy to diagnose when I tell you what the variables really are:

y ~ Shark attacks,
x₁ ~ Weather,
x₂ ~ Ice cream sales  

Estimating SEM coefficients

You might be tempted to take the SEM for the common cause structure and use that in your paper. I’d advise against it. Since we used all of our data in estimating the structure of the system, we don’t want to reuse the same data to estimate parameter coefficients. That would be double-dipping on our data. Having our cake and eating it too.

Instead, we need to either collect additional data (ideal) or reserve a split of our original data to use as the final fitting data. Since this is a simulated dataset, we can ‘collect’ all of the new data we want.

set.seed(666)

n <- 200

weather <- rnorm(n, mean = 20, sd = 5)

ice_cream_sales <- 50 + 2 * weather + rnorm(n, mean = 0, sd = 5)

shark_attacks <- 5 + 0.3 * weather + rnorm(n, mean = 0, sd = 5)

example01_newdata <- data.frame(shark_attacks, weather, ice_cream_sales)

# common cause
ex01_new_commoncause <- '
ice_cream_sales ~ weather
shark_attacks ~ weather
# shark_attacks ~ 1
'

ex01_sem_new_commoncause <- sem(ex01_new_commoncause, data = example01_newdata)

summary(ex01_sem_new_commoncause)

> summary(ex01_sem_new_commoncause)
lavaan 0.6-18 ended normally after 10 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         5

  Number of observations                           200

Model Test User Model:
                                                      
  Test statistic                                 0.000
  Degrees of freedom                                 0

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Regressions:
                    Estimate  Std.Err  z-value  P(>|z|)
  ice_cream_sales ~                                    
    weather            2.131    0.065   32.853    0.000
  shark_attacks ~                                      
    weather            0.419    0.061    6.852    0.000

Covariances:
                     Estimate  Std.Err  z-value  P(>|z|)
 .ice_cream_sales ~~                                    
   .shark_attacks       3.026    1.644    1.841    0.066

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .ice_cream_sals   24.452    2.445   10.000    0.000
   .shark_attacks    21.727    2.173   10.000    0.000

Our coefficient estimates are close, but not spot on the true parameter values we simulated. The reason is that the simulated errors are pretty high (sd = 5). Rather than fitting the final model once, we can borrow another technique from machine learning and repeatedly fit bootstrap replicates. This will both stabilize our estimates and provide a convenient, non-parametric confidence interval.

Boostrapping SEM coefficients

Content to come!

Extensions to SEM

Content to come!

Laten variables

Content to come!

Composite variables

Content to come!

Non-linear relationships

Content to come!

Hierarchical relationships

Content to come!

Background

This tutorial was originally developed for a graduate workshop at Yale University on 2024 07 25.

In an effort to make life easiest, I spent some time building out the pipeline entirely in R, including a sourceable function for converting spherical panos to hemispherical images (all available on this repo).

The easiest way to convert all of your spherical panos to hemispherical projections is to source the function from my github:

source("https://raw.githubusercontent.com/andisa01/Spherical-Pano-UPDATE/main/Spheres_to_Hemis.R")

When you source the script, it will install and load all necessary packages. It also downloads the masking file that we will use to black out the periphery of the images.

The script contains the function convert_spheres_to_hemis, which does exactly what is says. You’ll need to put all of your raw spherical panos into a subdirectory within your working directory. We can then pass the path to the directory as an argument to the function.

convert_spheres_to_hemis(focal_path = "./raw_panos/")

This function will loop through all of your raw panos, convert them to masked, north-oriented upward-facing hemispherical images and put them all in a folder called “masked_hemispheres” in your working directory. It will also output a csv file called “canopy_output.csv” that contains information about the image.

Below, I will walk through the steps of the workflow that happened in the convert_spheres_to_hemis function. If you just want to use the function, you can skip to the analysis. I’ve also written a script to do all of the conversion AND analysis in batch in the repo titled “SphericalCanopyPanoProcessing.R”.

library(tidyverse) # For data manipulation
library(exifr) # For extracting metadata

R is not really the best tool for working with image data. So, we’ll use the magick package to call the ImageMagick program from within R.

library(magick) # For image manipulation
# Check to ensure that ImageMagick was installed.
magick_config()$version

You should see a numeric version code like ‘6.9.12.3’ if ImageMagick was properly installed.

# ImageR also requires ImageMagick
library(imager) # For image display

For binarizing and calculating some canopy metrics, we will use Chiannuci’s hemispheR package, which we need to install from the development version.

library(devtools)
devtools::install_git("https://gitlab.com/fchianucci/hemispheR")
library(hemispheR) # For binarization and estimating canopy measures

To get started, you’ll need to have all of your raw equirectangular panoramas in a folder. I like to keep my raw images in a subdirectory called ‘raw_panos’ within my working directory. Regardless of where you store your files, set the directory path as the focal_path variable. For this tutorial we’ll process a single image, but I’ve included scripts for batch processing in the repo. Set the name of the image to process as the focal_image variable.

focal_path <- "./raw_panos/"
focal_image <- "PXL_20230519_164804198.PHOTOSPHERE_small.jpg"

focal_image_path <- paste0(focal_path, focal_image)
focal_image_name <- sub("\\.[^.]+$", "", basename(focal_image_path))

Let’s take a look at the equirectangular image.

pano <- image_read(focal_image_path)
pano # Visualize the pano

Note: One advantage of spherical panos is that they are large, and therefore, high resolution. The images from my Google Pixel 4a are 38 mega pixels. For this tutorial, I downsized the example pano to 10% resolution to make processing and visualizing easier. For your analysis, I’d recommend using full resolution images.

Spherical panoramas contain far more metadata than an average image. We can take a look at all of this additional information with the read_exif function.

read_exif(focal_image_path) %>%
  glimpse()

We’ll extract some of this information about the image, the date it was created, the georeference, and altitude, etc. to output alongside our canopy metrics. You can decide which elements are most important or useful for you.

xmp_data <- read_exif(focal_image_path) %>%
  select(
    SourceFile,
    Make,
    Model,
    FullPanoWidthPixels,
    FullPanoHeightPixels,
    SourcePhotosCount,
    Megapixels,
    LastPhotoDate,
    GPSLatitude,
    GPSLongitude,
    GPSAltitude,
    PoseHeadingDegrees
  )

The first image processing step is to convert the equirectangular panorama into a hemispherical image. We’ll need to store the image width dimension and the heading for processing. The pose heading is a particularly important feature that is a unique advantage of spherical panoramas. Since the camera automatically stores the compass heading of the first image of the panorama, we can use that information to automatically orient all of our hemispherical images such that true north is the top of the image. This is critical for analyses of understory light which requires plotting the sunpath onto the hemisphere.

# Store the pano width to use in scaling and cropping the image
pano_width <- image_info(pano)$width
image_heading <- read_exif(focal_image_path)$PoseHeadingDegrees

The steps to reproject the spherical panorama into an upward-looking hemispherical image go like this:

Crop the upper hemisphere (this is easy with smartphone spheres because the phone’s gyro ensures that the horizon line is always the midpoint of the y-axis).

Rescale the cropped image into a square to retain the correct scaling when reprojected into polar coordinate space.

Project into polar coordinate space and flip the perspective so that it is upward-looking.

Rotate the image so that the top of the image points true north and crop the image so that the diameter of the circle fills the frame.

We can accomplish all of those steps with the code below.

pano_hemisphere <- pano %>%
  # Crop to retain the upper hemisphere
  image_crop(geometry_size_percent(100, 50)) %>%
  # Rescale into a square to keep correct scale when projecting in to polar coordinate space
  image_resize(geometry_size_percent(100, 400)) %>%
  # Remap the pixels into polar projection
  image_distort("Polar",
                c(0),
                bestfit = TRUE) %>%
  image_flip() %>%
  # Rotate the image to orient true north to the top of the image
  image_rotate(image_heading) %>%
  # Rotating expands the canvas, so we crop back to the dimensions of the hemisphere's diameter
  image_crop(paste0(pano_width, "x", pano_width, "-", pano_width/2, "-", pano_width/2))

The resulting image looks funny because the outer pixels are extended by interpolation and we’ve rotated the image which leaves white space at the corners. Most analyses define a bounding perimeter to exclude any pixels outside of the circular hemisphere; so, the weird border shouldn’t matter. But, we can add a black mask to make the images look better.

I’ve included a vector file for a black mask to lay over the image in the repo.

# Get the image mask vector file
image_mask <- image_read("./HemiPhotoMask.svg") %>%
  image_transparent("white") %>%
  image_resize(geometry_size_pixels(width = pano_width, height = pano_width)) %>%
  image_convert("png")

masked_hemisphere <- image_mosaic(c(pano_hemisphere, image_mask))

masked_hemisphere

We’ll store the masked hemispheres in their own subdirectory. This script makes that directory, if it doesn’t already exists and writes our file into it.

if(dir.exists("./masked_hemispheres/") == FALSE){
  dir.create("./masked_hemispheres/")
} # If the subdirectory doesn't exist, we create it.

masked_hemisphere_path <- paste0("./masked_hemispheres/", focal_image_name, "hemi_masked.jpg") # Set the filepath for the new image

image_write(masked_hemisphere, masked_hemisphere_path) # Save the masked hemispherical image

At this point, you can use the hemispherical image in any program you’d like either in R or other software. For this example, I’m going to use Chiannuci’s hemispheR package in order to keep this entire pipeline in R.

The next step is to import the image. hemispheR allows for lots of fine-tuning. Check out the docs to learn what all of the options are. These settings most closely replicate the processing I used in my 2021 paper.

fisheye <- import_fisheye(masked_hemisphere_path,
                          channel = '2BG',
                          circ.mask = list(xc = pano_width/2, yc = pano_width/2, rc = pano_width/2),
                          gamma = 2.2,
                          stretch = FALSE,
                          display = TRUE,
                          message = TRUE)

Now, we need to binarize the images, converting all sky pixels to white and everything else to black (at least as close as possible). Again, there are lots of options available in hemispheR. You can decides which settings are right for you. However, I would suggest keeping the zonal argument set to FALSE. The documentation describes this argument as:

zonal: if set to TRUE, it divides the image in four sectors
(NE,SE,SW,NW directions) and applies an automated classification
separatedly to each region; useful in case of uneven light conditions in
the image

Because spherical panoramas are exposing each of the 36 images separately, there is no need to use this correction.

I also suggest keeping the export argument set to TRUE so that the binarized images will be automatically saved into a subdirectory named ‘results’.

binimage <- binarize_fisheye(fisheye,
                 method = 'Otsu',
                 # We do NOT want to use zonal threshold estimation since this is done by the camera
                 zonal = FALSE,
                 manual = NULL,
                 display = TRUE,
                 export = TRUE)

Unfortunately, hemispheR does not allow for estimation of understory light metrics like through-canopy radiation or Global Site Factors. If you need light estimates, you’ll have to take the binarized images and follow my instructions and code for implementing Gap Light Analyzer.

Assuming all you need is canopy metrics, we can continue with hemispheR and finalize the whole pipeline in R. We estimate canopy metrics with the gapfrac_fisheye() function.

gapfrac <- gapfrac_fisheye(
  binimage,
  maxVZA = 90,
  # Spherical panoramas are equidistant perforce
  lens = "equidistant",
  startVZA = 0,
  endVZA = 90,
  nrings = 5,
  nseg = 8,
  display = TRUE,
  message = TRUE
)

Finally, we can estimate the canopy metrics with the canopy_fisheye() function, join those to the metadata from our image, and output our report.

canopy_report <- canopy_fisheye(
  gapfrac
)

output_report <- xmp_data %>%
  bind_cols(
    canopy_report
  ) %>%
  rename(
    GF = x,
    HemiFile = id
  )

glimpse(output_report)

Rows: 1
Columns: 32
$ SourceFile            "./raw_panos/PXL_20230519_164804198.PHOTOSPHERE_small.jpg"
$ Make                  "Google"
$ Model                 "Pixel 4a"
$ FullPanoWidthPixels   8704
$ FullPanoHeightPixels  4352
$ SourcePhotosCount     36
$ Megapixels            0.37845
$ LastPhotoDate         "2023:05:19 16:49:57.671Z"
$ GPSLatitude           41.33512
$ GPSLongitude          -72.91103
$ GPSAltitude           -23.1
$ PoseHeadingDegrees    86
$ HemiFile              "PXL_20230519_164804198.PHOTOSPHERE_smallhemi_masked.jpg"
$ Le                    2.44
$ L                     3.37
$ LX                    0.72
$ LXG1                  0.67
$ LXG2                  0.55
$ DIFN                  9.981
$ MTA.ell               19
$ GF                    4.15
$ VZA                   "9_27_45_63_81"
$ rings                 5
$ azimuths              8
$ mask                  "435_435_434"
$ lens                  "equidistant"
$ channel               "2BG"
$ stretch               FALSE
$ gamma                 2.2
$ zonal                 FALSE
$ method                "Otsu"
$ thd                   116

Be sure to check out my prior posts on working with hemispherical images.

In this post, I will show some methods of displaying mixed effect regression models and associated uncertainty using non-parametric bootstrapping. This is kind of a follow-up to my previous post on visualizing custom main effect models.

Introduction

Mixed models have quickly become the model du jour in many observation-oriented fields because these models obviate many of the issues of pseudoreplication inherent in blocked or repeated measures experiments and structured data.

They do this by treating the levels of categorical variables not as unique instances to be parameterized individually, but as random samples from an infinite distribution of levels. Instead of wasting degrees of freedom estimating parameter values for each level, we only need to estimate that global distribution (which requires only a handful of parameters) and instead focus our statistical power on the variables of interest.

Thanks to packages like nlme and lme4, mixed models are simple to implement. For all of their virtues, mixed models can also be a pain to visualize and interpret. Although linear mixed models are conceptually similar to the plain old ordinary least-squares regression we know and love, they harbor a lot more math under the hood, which can be intimidating.

One of the reasons mixed models are difficult to intuitively visualize is because they allow us to manage many levels of uncertainty. Depending on the focus of our analyses, we usually want to focus on certain aspects of the trends and associated uncertainty. For instance, an ecologist might be interested in the effect of nutrient input across many plots, but not interested in the difference between plots (i.e. traditional random effect). Or, an educator might be interested in the effect of different curricula, but not the difference between specific classes within specific schools (i.e. nested random effects). Or, a physician might be interested in the effect of a long-term treatment on a patient after accounting for baseline difference between patients (i.e. repeated measures).

In this tutorial, I’m going to focus on how to visualize the results of mixed effect models from lme4 using ggplot2. You can also clone the annotated code from my Github.

First, load in the necessary libraries.

library(tidyverse)
library(lme4)
library(ggsci)
library(see)
library(cowplot)
theme_set(theme_classic()) # This sets the default ggplot theme

To begin, I am going to simulate an experiment with 10 experimental units each containing 100 observations. These could be 10 plots with 100 random samples, or 10 schools with 100 student test scores, or the records of 10 patients from each of 100 follow-up visits. Each of the experimental units will ultimately get its own intercept and slope effect coefficient. The rand_eff data frame is essentially a Z matrix in classic mixed model notation. For this example, I’ll assume that the intercepts come from a distribution with standard deviation of 20 and the slopes from a distribution with standard deviation of 0.5. The random effects define the variation of the experimental unit around the main effect, so the mean of these distributions is necessarily 0.

set.seed(666)
rand_eff <- data.frame(unit = as.factor(seq(1:10)),
            b0 = rnorm(10, 0, 20),
            b1 = rnorm(10, 0, 0.5))

We can now join our random effect matrix to the full dataset and define our y values as y_i = B0_i + b0_j + B1x_i + b1x_j + ε.

X <- expand.grid(unit = as.factor(seq(1:10)), obs = as.factor(seq(1:100))) %>%
  left_join(rand_eff,
            by = "unit") %>%
  mutate(x = runif(n = nrow(.), 0, 10),
         B0 = 20,
         B1 = 2,
         E = rnorm(n = nrow(.), 0, 10)) %>%
  mutate(y = B0 + b0 + x * (B1 + b1) + E)

Here’s a look at the data.

X %>%
  ggplot(aes(x = x, y = y, col = unit)) +
  geom_point() +
  facet_wrap(vars(unit))

Random intercept model

For demonstration, let’s first assume that we are primarily interested in the overall slope of the relationship. For instance, if these are 10 field plots, we might want to know the effect of adding 1 unit of nutrient fertilizer, regardless of the baseline level of nutrients in a given plot.

We can do this by fitting a random intercept model and then looking at the summary of the resulting model.

lmer1 <- lmer(y ~ x + (1|unit), data = X)

summary(lmer1)

> summary(lmer1)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (1 | unit)
   Data: X

REML criterion at convergence: 7449.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.88156 -0.68745  0.01641  0.71022  2.84532 

Random effects:
 Groups   Name        Variance Std.Dev.
 unit     (Intercept) 811.59   28.488  
 Residual              94.63    9.728  
Number of obs: 1000, groups:  unit, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  18.2652     9.0301   2.023
x             2.0091     0.1077  18.651

Correlation of Fixed Effects:
  (Intr)
x -0.060

We can see that the fitted model does a good job estimating the fixed effect slope (B1), which we simulated with a coefficient of 2 as 2.0091. However, the model is underestimating the fixed effect intercept (B0) as 18.3 and overestimating the standard deviation of the random effect slopes (b1) as 28.5, when we simulated those values as 20 and 20.

If we think we could live with that fit, how would we go about visualizing our model?

Here is a look at our data with a linear regression fit to each experimental unit. It is clear that there is a wide spread in the intercepts, but the slopes are similar.

X %>%
ggplot(aes(x = x, y = y, col = unit)) +
geom_point() +
geom_smooth(method = 'lm', se = F)

Marginal (Fixed effect) versus Conditional (Fixed + Random effect)

We might be tempted to use this built-in regression by group from ggplot as a visualization of the mixed model. However, this would be WRONG!!! GGplot is fitting an ordinary least squares regression without accounting for the random effect. That means that the estimates and the confidence intervals do not reflect our model. In this case, the estimates might be pretty close since our samples sizes across species are pretty even, but this could be wildly off, or even opposite, of mixed model slope estimate.

In a prior post, I showed how we can use the predict function to display our custom models in ggplot. In the case of mixed effect models, you can predict both the marginal and conditional values. The marginal value is the fixed effect. The conditional value is the mixed effect of the fixed and random effects.

In other words, the marginal effect is asking “What would I expect y to be for a given x without knowing which experimental unit it came from?” whereas the conditional effect is asking “What would I expect y to be for a given x from a given experimental unit?”

We can specify which prediction we want with the random effect formula argument re.form:

X <- X %>% 
  mutate(fit.m = predict(lmer1, re.form = NA),
         fit.c = predict(lmer1, re.form = NULL))

The simplest visualization would be to display the marginal fit on the raw values.

X %>%
  ggplot(aes(x = x, y = y)) +
    geom_point(pch = 16, col = "grey") +
    geom_line(aes(y = fit.m), col = 1, size = 2) +
    coord_cartesian(ylim = c(-40, 100))

However, this is a bit misleading because it underrepresents our confidence in the slope by making it look like the residuals are huge.

But the residuals, and our confidence in the fit, is based on the conditional residual variance, which is much tighter. We can see that easily when we look at the conditional fits. This is one option for visualization, but it highlights the wrong element if our primary interest is the overall slope trend.

X %>%
  ggplot(aes(x = x, y = y, col = unit)) +
  geom_point(pch = 16) +
  geom_line(aes(y = fit.c, col = unit), size = 2) +
  facet_wrap(vars(unit)) +
  coord_cartesian(ylim = c(-40, 100))

Displaying the conditional fits on the same facet helps. Now we can see the variation in the conditional intercepts, but as a tradeoff it makes it difficult to get a sense of the residual variance because there are too many points.

X %>%
  ggplot(aes(x = x, y = y, col = unit)) +
  geom_point(pch = 16) +
  geom_line(aes(y = fit.c, col = unit), size = 2)  +
  coord_cartesian(ylim = c(-40, 100))

Instead, I think it makes more sense to display the conditional residuals around the marginal effect. You can kind of think of this as collapsing all of the conditional fits from the previous plot into the single marginal fit. We can do this by extracting the residuals (which are the conditional residuals) and then displaying the points as the marginal fit plus the residuals.

X <- X %>%
  mutate(resid = resid(lmer1))

X %>%
  ggplot(aes(x = x, y = fit.m + resid, col = unit)) +
  geom_point(pch = 16) +
  geom_line(aes(y = fit.m), col = 1, size = 2)  +
  coord_cartesian(ylim = c(-40, 100))

In some cases, we might also want to give the reader a sense of the variation in the conditional intercepts. For instance, the fact that the slope is so consistent across a wide range of baselines might actually increase our confidence in the relationship even further.

There are a couple of ways to simultaneously display both our confidence in the fit of the marginal trend and the variance in the conditional fits.

Depending on the number of conditional units, one option is to display the conditional fits below the scatter plot of the conditional residuals.

X %>%
  ggplot(aes(x = x, y = fit.m + resid)) +
  geom_line(aes(y = fit.c, col = unit), size = 1) +
  geom_point(pch = 16, col = "grey") +
  geom_line(aes(y = fit.m), col = 1, size = 2) +
  coord_cartesian(ylim = c(-40, 100))

Another option is to display a density plot or histogram of the estimated conditional intercepts (also known as the Best Linear Unbiased Predictors or BLUPs). In a random effect framework, we are assuming that the conditional intercepts are samples of some infinite distribution of intercepts, so this histogram from the BLUPs of our model is essentially an empirical representation of that idealized distribution. (Alternatively, we could also simply plot the idealized distribution as a normal distribution from the estimated variance of the random effect, but I like the empirical density plot because it also give a sense of when our conditionals do NOT conform to the assumption of being samples from a normal distribution).

We can extract the BLUPs from the model object (b0_hat) and add those to the model estimate of the marginal intercept (B0_hat) to get the estimated conditional intercepts. This is our data frame of conditional estimates.

Cond_DF <- as.data.frame(ranef(lmer1)) %>% transmute(unit = grp, b0_hat = condval) %>% mutate(Intercept_cond = b0_hat + summary(lmer1)$coef[1,1])

X %>%
  ggplot(aes(x = x, y = fit.m + resid)) +
  geom_point(pch = 16, col = "grey") +
  geom_violinhalf(data = Cond_DF, aes(x = 0, y = Intercept_cond), trim = FALSE, width = 3, fill = NA) +
  geom_line(aes(y = fit.m), col = 1, size = 2) +
  coord_cartesian(ylim = c(-40, 100))

Random slope and intercept model

Now let’s imagine that we are not satisfied with the random intercept model and also want to fit a random slope parameter. In this case, we want to estimate the distribution of slopes for all experimental units across the values of x.

lmer2 <- lmer(y ~ x + (x|unit), data = X)

We can see that the values from the model are getting much closer to the known values that we simulated.

summary(lmer2)

> summary(lmer2)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (x | unit)
   Data: X

REML criterion at convergence: 7442.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2071 -0.6957  0.0268  0.7067  2.8482 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 unit     (Intercept) 846.3021 29.0913       
          x             0.2151  0.4638  -0.27
 Residual              93.0853  9.6481       
Number of obs: 1000, groups:  unit, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  18.3474     9.2202    1.99
x             2.0075     0.1816   11.05

Correlation of Fixed Effects:
  (Intr)
x -0.250

Correlation of random effects

One important addition from the random intercept-only model is the estimate for the correlation between the distribution of the random slopes and random intercepts (which the model estimates as -0.268, see output below). Because we simulated these data, we know that there is no true correlation between the unit slopes and intercepts. But, because we have a small number of units, we just happened to have an emergent correlation.

summary(lmer2)$varcor

cor.test(rand_eff$b0, rand_eff$b1)

> summary(lmer2)$varcor
 Groups   Name        Std.Dev. Corr  
 unit     (Intercept) 29.09127       
          x            0.46378 -0.268
 Residual              9.64807       
> 
> cor.test(rand_eff$b0, rand_eff$b1)

	Pearson's product-moment correlation

data:  rand_eff$b0 and rand_eff$b1
t = -1.2158, df = 8, p-value = 0.2587
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.8205150  0.3123993
sample estimates:
       cor 
-0.3949022 

If you had a reason to assume NO correlation between your random effects, you could specify that as (x||unit) in this way:

lmer3 <- lmer(y ~ x + (x||unit), data = X)

summary(lmer3)

summary(lmer3)$varcor

> summary(lmer3)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + ((1 | unit) + (0 + x | unit))
   Data: X

REML criterion at convergence: 7443.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1812 -0.6951  0.0223  0.7074  2.8523 

Random effects:
 Groups   Name        Variance Std.Dev.
 unit     (Intercept) 833.2681 28.8664 
 unit.1   x             0.2092  0.4574 
 Residual              93.1011  9.6489 
Number of obs: 1000, groups:  unit, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)   18.325      9.149   2.003
x              2.010      0.180  11.165

Correlation of Fixed Effects:
  (Intr)
x -0.035
> 
> summary(lmer3)$varcor
 Groups   Name        Std.Dev.
 unit     (Intercept) 28.86638
 unit.1   x            0.45738
 Residual              9.64889

As with the random intercept model, we can use the predict function to get expected values of y based on the marginal or conditional estimates. Note that re.form = NULL is the same as re.form = ~ (x|unit).

X <- X %>% 
  mutate(fit2.m = predict(lmer2, re.form = NA),
         fit2.c = predict(lmer2, re.form = NULL),
         resid2 = resid(lmer2))

As with the random intercept model, one way to visualize the model is to show the conditional intercept/slopes as fitted lines and the conditional residuals as points.

pmain_lmer2 <- X %>%
  ggplot(aes(x = x, y = fit2.m + resid2)) +
  geom_line(aes(y = fit2.c, col = unit), size = 1) +
  geom_point(pch = 16, col = "grey") +
  geom_line(aes(y = fit2.m), col = 1, size = 2) +
  coord_cartesian(ylim = c(-40, 100))
pmain_lmer2

Visualizing the random effect variance gets a bit more difficult with two random parameters. One strategy I like is to include an additional plot of the correlation and distribution of the random effects.

Basically, we make a scatter plot of the BLUPs of slopes and intercepts, Then we make histograms for each and add those to the margins of the plot. (I’m relying heavily on Claus Wilke’s post to code the histograms at the margins of the plot). Finally, we patch it all together with cowplot.

Cond_DF2 <- as.data.frame(ranef(lmer2)) %>% 
  transmute(unit = grp,
            term = case_when(term == "(Intercept)" ~ "b0_hat",
                             term == "x" ~ "b1_hat"),
            value = condval) %>%
  pivot_wider(id_cols = "unit", names_from = "term", values_from = "value") %>%
  mutate(Intercept_cond = b0_hat + summary(lmer2)$coef[1,1],
         Slope_cond = b1_hat + summary(lmer2)$coef[2,1])

pmain <- Cond_DF2 %>%
  ggplot(aes(x = Intercept_cond, y = Slope_cond)) +
    geom_point(aes(col = unit), size = 3) +
    geom_density2d(bins = 4, col = "grey", adjust = 3)

xdens <- axis_canvas(pmain, axis = "x") +
  geom_density(data = Cond_DF2, aes(x = Intercept_cond), fill = "grey", col = NA, trim = FALSE, adjust = 2)

ydens <- axis_canvas(pmain, axis = "y", coord_flip = TRUE) +
  geom_density(data = Cond_DF2, aes(x = Slope_cond), fill = "grey", col = NA, trim = FALSE, adjust = 2) +
  coord_flip()

p1 <- insert_xaxis_grob(pmain, xdens, grid::unit(.2, "null"), position = "top")
p2 <- insert_yaxis_grob(p1, ydens, grid::unit(.2, "null"), position = "right")

pinsert_lmer2 <- ggdraw(p2)

plot_grid(
  pmain_lmer2,
  pinsert_lmer2,
  nrow = 1
)

Bootstrapping uncertainty

For many researchers, one of the most frustrating aspects of mixed models is that estimating confidence intervals and testing the significance of parameters is not straight forward. I highly encourage that folks take a look at Ben Bolker’s thorough considerations on the topic. Dr. Bolker suggests many different problems and solutions depending on the structure of your model.

I think that the most generalized solution is to use non-parametric bootstrapping. This method essentially asks the question, “How would our model fit change if we could go back in time, select different samples, and then rerun our analysis?”

We can’t go back in time, but maybe we CAN assume that our original samples were representative of the population. If so, instead of resampling the actual population, we could resample our original observations WITH REPLACEMENT to approximate resamples of the population. If we do this many time, we can then make intuitive statements like, “95% of the time, if we reran this experiment, we’d expect the main effect to be between X and X value.”

It is important to stop an consider what re-doing your data generation process would look like. For instance, imagine our mock data had come from 1000 independent random observation that we then categorized into “units” to control for autocorrelation after the fact. If we re-ran the process to generate a new dataset, we may not always get the same number of observations in each “unit”.

However, if our mock data came from an experiment where we planted 100 trees in each of 10 “units”, then when we re-ran the experiment, we could control the number of individuals per unit. We would also need to consider if we would always choose to plant in the same 10 units, or if we would also choose units at random.

The structure of the data generation process can guide our bootstrap resampling strategy. In the first example, we could simply bootstrap all individual observations (although we may need to worry about non-convergence and small sample sizes). In the second example, where unit choice is constrained, we might decide to bootstrap within units. If, in the second example, we could also randomize units, we should probably take a hierarchical approach, first bootstrapping the units and then bootstrapping the observations within each unit.

NOTE: The problem with non-parametric boostrapping of this kind is that it can be computationally expensive. One trick is to parallelize the bootstraps across all of your computer’s processors. By default, R uses one processor, so it will fit one bootstrap iteration at a time, sequentially. But the bootstraps are independent and the order doesn’t matter. If you computer has 8 cores, there is no reason not to fit 8 models simultaneously on 8 processors in 1/8 of the time. Unfortunately, setting up for parallel processing can be an adventure of its own. I won’t detail it here, but will try to dedicate a post on it in the future. If you have a very large dataset and can’t run bootstraps in parallel, you might consider some of the other methods suggested by Dr. Bolker.

Since our dataset is fairly small and simple, I’ll demonstrate how we can use bootstrapping to simultaneously estimate confidence intervals of our model parameters and visualize error bands.

Just keep in mind that if you are fitting models with complex random effect designs, you’ll have to think critically about which elements and levels of variance are most important for your data story. Hopefully, these examples will at least get you started and inspired!

The bootstrapping process begins by initiating an empty dataframes to accept the parameter estimates for the fixed effect coefficients and random effect variances. The number of rows in the dataframe will be the same as the number of bootstrap iterations, so we set that first. The number of iterations is the number of times we want to simulate re-doing our data generation. The convention is 1000, but the more the merrier!

Number_of_boots <- 1000

The number of columns for the dataframes will equal the number of fixed effect coefficients and random effect variances. We can extract these from the initial model. First we extract the coefficients, then transpose the table into wide format.

# Extract the fixed effect coefficients.
FE_df <- fixef(lmer2) %>% 
  t() %>%
  as.data.frame()

# Extract the random effects variance and residual variance
RE_df <- VarCorr(lmer2) %>%
  as.data.frame() %>%
  unite("Level", -c(vcov, sdcor)) %>%
  select(-vcov) %>%
  t() %>%
  as.data.frame()

Next, we create empty dataframes to take our bootstraps.

BS_params <- data.frame(matrix(nrow = Number_of_boots, ncol = ncol(FE_df)))
colnames(BS_params) <- colnames(FE_df)

BS_var <- data.frame(matrix(nrow = Number_of_boots, ncol = ncol(RE_df)))
colnames(BS_var) <- RE_df["Level",]

In addition, we will be predicting marginal values from each model. So, we need to create a prediction dataframe with an empty column to store the predicted values. For this example, we only need to predict ŷ values for a handful of x values that represent the range of xs. I chose to use 10-quantiles because I want to be able to fit a non-linear confidence band later. If this was a non-linear fit, we might want even more prediction values.

BS_pred <- expand.grid(x = quantile(X$x, probs = seq(0, 1, length.out = 10)),
                       iterration = 1:Number_of_boots,
                       pred = NA)

Finally, we can write a loop that creates a resampled dataset (with replacement) and fits the original model formula to the new dataset. From the new model, we can then extract the fixed and random effects and predict ŷs for the subset of x values. All of these get stored in their respective dataframes, indexed by the iteration number.

for(i in 1:Number_of_boots){
BS_X <- slice_sample(X, prop = 1, replace = TRUE)
BS_lmer <- lmer(formula = lmer2@call$formula,
data = BS_X)

BS_params[i,] <- BS_lmer %>%
fixef() %>%
t() %>%
as.data.frame()

BS_var[i,] <- BS_lmer %>%
VarCorr() %>%
as.data.frame() %>%
.$sdcor

BS_pred[which(BS_pred$iterration == i),]$pred <- predict(BS_lmer,
newdata = BS_pred[which(BS_pred$iterration == i),],
re.form = ~0)
}

Now we have a dataframe of the marginal (i.e. fixed effect) intercept and slope parameter estimates from 1000 models fit to bootstraps of our original data.

head(BS_params)

> head(BS_params)
  (Intercept)        x
1    18.06942 2.096240
2    18.18005 2.070043
3    18.13506 2.093110
4    18.77862 1.928048
5    18.47963 2.013831
6    18.28875 2.005947

One way to get a sense of our confidence in these parameter estimates is to take a look at their distributions.

BS_hist_x <- BS_params %>%
  ggplot(aes(x = x)) +
  geom_histogram()

BS_hist_intercept <- BS_params %>%
  ggplot(aes(x = `(Intercept)`)) +
  geom_histogram()

BS_hists <- plot_grid(
  BS_hist_intercept,
  BS_hist_x,
  nrow = 1)

BS_hists

These histograms tell us so much more than a typical confidence interval because we can see the full distribution. We can see that the baseline effect on y given x = 0 is around 18.5 and we are very confident it is between 17 and 20. We can also see that the effect of 1 unit change in x is expected to yield about a 2.0 unit change in y and we are extremely confident that the slope is positive and greater than 1.7 but probably less than 2.3.

Although we won’t directly visualize the random effect variance, we can see the estimates in the BS_var dataframe.

BS_var %>% head()

> BS_var %>% head()
  unit_(Intercept)_NA unit_x_NA unit_(Intercept)_x Residual_NA_NA
1            28.30440 0.4427450         0.01009018       9.786112
2            30.31353 0.6265168        -0.42315981       9.365630
3            30.33369 0.5261389        -0.78128072       9.320844
4            28.51896 0.6079722        -0.08575830       9.739881
5            30.50690 0.5835157        -0.64068868       9.718633
6            29.17012 0.5182395        -0.28244154       9.930305

Here, the first column is the estimated variance attributed to the random intercepts, the second column is the variance estimate of the random slopes, and the third column is the correlation between the random effects. The fourth and final column is the residual variation after accounting for the random effects (i.e. the conditional residual variation). This is information you’d want to include in a table of the model output.

We also have a dataframe of 10 predicted ŷs for each iteration.

head(BS_pred, n = 20)

> head(BS_pred, n = 20)
             x iterration     pred
1  0.008197445          1 18.08660
2  1.098954093          1 20.37309
3  2.231045775          1 22.74622
4  3.266109165          1 24.91596
5  4.323574828          1 27.13267
6  5.667717615          1 29.95031
7  6.672807375          1 32.05722
8  7.748743587          1 34.31264
9  8.793447507          1 36.50259
10 9.968969584          1 38.96677
11 0.008197445          2 18.19702
12 1.098954093          2 20.45493
13 2.231045775          2 22.79841
14 3.266109165          2 24.94104
15 4.323574828          2 27.13003
16 5.667717615          2 29.91247
17 6.672807375          2 31.99305
18 7.748743587          2 34.22028
19 8.793447507          2 36.38286
20 9.968969584          2 38.81624

Rather than using a traditional confidence band (which basically reduces the distribution of our bootstraps down to two points: high and low), I prefer to actually show all of the iterations and let the density of the lines make a kind of optical illusion of a confidence band.

Since our estimates are pretty tightly distributed, I’m using cowplot to show what this looks like while also zooming in to a section.

plot_grid(
  BS_pred %>%
    ggplot(aes(x = x, y = pred)) +
    geom_line(aes(group = iterration), alpha = 0.1, col = "grey50") +
    geom_line(data = X,
              aes(x = x, y = fit2.m)) +
    geom_rect(aes(ymin = 35, ymax = 40,
                  xmin = 8, xmax = 10),
              col = "firebrick",
              fill = NA,
              size = 2),
  BS_pred %>%
    ggplot(aes(x = x, y = pred)) +
    geom_line(aes(group = iterration), alpha = 0.1, col = "black") +
    geom_line(data = X,
              aes(x = x, y = fit2.m),
              col = "grey", 
              size = 2) +
    coord_cartesian(xlim = c(8, 10),
                    ylim = c(35, 40)) +
    geom_rect(aes(ymin = 35, ymax = 40,
                  xmin = 8, xmax = 10),
              col = "firebrick",
              fill = NA,
              size = 2) +
    theme(axis.line = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank()) +
    labs(x = "", y = ""),
  nrow = 1
)

Of course, we can also create a more traditional 90% confidence band by summarizing the 5th and 95th percentiles of ŷ across each x value.

BS_pred %>%
  group_by(x) %>%
  summarise(hi = quantile(pred, 0.95),
            lo = quantile(pred, 0.05)) %>%
  ggplot(aes(x = x)) +
    geom_ribbon(aes(ymin = lo, ymax = hi),
                fill = "grey50",
                alpha = 0.3) +
    geom_line(data = X,
              aes(x = x, y = fit2.m))

Putting it all together

Putting it all together, here is my preferred visualization of a mixed effect model with random intercepts and slopes, using bootstrapping to display uncertainty.

BS_ci_lines <- BS_pred %>%
  ggplot(aes(x = x, y = pred)) +
    geom_line(aes(group = iterration), alpha = 0.1, col = "grey") +
    geom_line(data = X,
              aes(x = x, y = fit2.m)) +
    geom_point(data = X,
               aes(x = x, y = fit2.m + resid2, col = unit),
               alpha = 0.3,
               pch = 16)

plot_grid(
  BS_ci_lines,
  plot_grid(
    pinsert_lmer2,
    plot_grid(
      BS_hist_intercept,
      BS_hist_x,
      nrow = 1,
      labels = c("C", "D")),
    nrow = 2,
    rel_heights = c(1, 0.7),
    labels = c("B", NA)
  ),
  nrow = 1,
  rel_widths = c(1, 0.7),
  labels = c("A", NA)
)plot_grid(
  BS_ci_lines,
  plot_grid(
    pinsert_lmer2,
    plot_grid(
      BS_hist_intercept,
      BS_hist_x,
      nrow = 1,
      labels = c("C", "D")),
    nrow = 2,
    rel_heights = c(1, 0.7),
    labels = c("B", NA)
  ),
  nrow = 1,
  rel_widths = c(1, 0.7),
  labels = c("A", NA)
)

Or, if we want to get really fancy with it, we could inset everything into one plot panel.

BS_ci_lines +
  coord_cartesian(ylim = c(-30, 100)) +
  annotation_custom(ggplotGrob(BS_hists),
                    xmin = 5,
                    xmax = 10,
                    ymin = -30,
                    ymax = 5
  ) +
  annotation_custom(ggplotGrob(pinsert_lmer2),
                    xmin = 0,
                    xmax = 4,
                    ymin = 50,
                    ymax = 110
  )

I am deep in the Alaskan Arctic,  300 miles from the nearest road system, attempting to conduct the kind of science that usually requires a specialized laboratory. We rowed 30 miles of meandering flatwater today, bringing our total to 200 river miles in 12 days since we landed at a lonely gravel bar on the headwaters of Ambler River in Gates of the Arctic National Park.

Mosquitoes spangle the tent canopy arching over me. Backlit by summer solstice sun, the silhouettes of the insects make an inverted night sky of shifting constellations. The sun never sets on the banks of the Kobuk River this time of year. It hangs high above the horizon even now at 11 p.m., transforming my tent into a solar oven as I, ironically, work to uncover the secrets of a frog that can turn into ice.

… Read the rest of the article here.

The other day, I wanted to tweet out a map showing the distribution of some wood frog tissue samples compared to the entire range of the species. I’m not much for GIS and I didn’t need anything complicated so I wanted to plot in R. If you find yourself in a similar situation needing to plot a simple range map with GGplot, then this post is for you.

This post WILL NOT explain how to estimate a species range. There’s a lot more involved in that process and there are dedicated packages like rangemappr to help.

In this post, we will take a pre-made polygon of a species range and slap it onto a basemap.

The first task is to find a polygon of the range of your favorite species. Unless you already have a shapefile in your possession, one of the easiest places to acquire one is through the IUCN database. You can follow the links to go down the taxonomic rabbit hole for specific species. For this tutorial I’ll be using the red salamander (Pseudotriton ruber), just because I think that they are gorgeous.

Here’s the  IUCN page for the red salamander. In the top right of the page, you’ll see a drop down button labelled, “Download”. You will want to select and download the “Range data – Polygons (SHP)” option. You will need to agree to some conditions in order to download the file, but it is quick and painless. Once you download the zipped folder, extract it into your project directory.

Here are the packages we will be using.

library(tidyverse)
library(rgdal)
library(ggspatial)

First, we need to do a bit of work to get the shapefiles into R. First, we will use rgdal package to read in the shapefiles. Then, we will coerce that object into a dataframe that can play nicely with GGplot.

# Read in the shapefile
PSRU_shape <- readOGR(dsn = "./PSRU_range", layer = "data_0")
# Coerce into a dataframe to play nicely with GGplot
PSRU_shape_df <- fortify(PSRU_shape)

Now we can start mapping. We will use the ggmap package in GGplot that all comes bundled in the tidyverse. First, we define regions for our basemap. In this case, red salamanders range is entirely within the lower 48 United States. Then we plot the basemap and layer on the range map on top. It’s really that simple.

Low48_map <- map_data("state") # Create basemap form GGplot

ggplot(Low48_map, aes(x = long, y = lat, group = group)) +
geom_polygon() +
geom_polygon(data = PSRU_shape_df, fill = "orangered")

But we can make this look a lot nicer. First off, we can ditch the ugly GGplot default theme and add our own color scheme. We can also force GGplot to use a polyconic projection (you’ll need to make sure that you also have the mapproj package installed).

ggplot(Low48_map, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = PSRU_shape_df, fill = "darkorange3", alpha = 0.8) +
theme_void() +
theme(panel.background = element_rect(fill = "cornsilk")) +
coord_map(projection = "polyconic")

Beyond the basics:

For species with larger ranges, you might want to incorporate more of the continent. We can easily do that by defining larger regions for our basemap. Here, I am still plotting the Lower 48 states basemap to get the state borders.

NAm_map <- map_data("world", region = c("Mexico", "Canada")) # Create basemap form GGplot

ggplot(NAm_map, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = Low48_map, fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = PSRU_shape_df, fill = "darkorange3", alpha = 0.8) +
theme_void() +
theme(panel.background = element_rect(fill = "cornsilk")) +
coord_map(projection = "gilbert")

One issue is that it is extremely difficult to crop in if, for instance, we don’t need all of Canada’s norther islands. Using xlim and ylim truncates the polygons in weird ways. And you should NEVER crop with xlim/ylim, anyway.

ggplot(NAm_map, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = Low48_map, fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = PSRU_shape_df, fill = "darkorange3", alpha = 0.8) +
theme_void() +
theme(panel.background = element_rect(fill = "cornsilk")) +
coord_map(projection = "gilbert") +
xlim(-125, -60) +
ylim(25, 64)

The better way to crop is to set the xlim and ylim of the clipping mask with coord_cartesian(). But, then we lose the projection.

ggplot(NAm_map, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = Low48_map, fill = "cornsilk4", col = "cornsilk") +
geom_polygon(data = PSRU_shape_df, fill = "darkorange3", alpha = 0.8) +
theme_void() +
theme(panel.background = element_rect(fill = "cornsilk")) +
coord_map(projection = "mercator") +
coord_cartesian(xlim = c(-125, -60), ylim = c(25, 52))

There are a handful of work-arounds to the problem.

First, there is a code workaround here, but it is a bit complicated.

Second, you can simply export the image at an appropriate size to get a resulting ratio that works. One potential issue with this is that GGplot exports ALL of the information in the PDF, even the polygons outside of the image. In this case, that means a lot of data points to outlines all those tiny Arctic Islands. This can make for some really large PDF files. A very slick solution to this is to rasterize the maps with ggraster before exporting the images. You can still export the image at twice or three times the size to retain high resolution, but the file size will be much smaller.

Third, you could export the uncropped version as a pdf (which is a vector graphic) and open it in a vector graphics program like Adobe Illustrator or Inkscape. From there, you can crop in where ever you want.

Take me to the photos!

Step 1: Come up with a hair-brained scheme.

My labmate Yara and I had been dreaming up the idea studying wood frog genomes from across the species’ range since she started her PhD. Wood frogs have the largest range of any North American amphibian. They also happen to be the only North American amphibian that can survive North of the Arctic circle.

Dr. Julie Lee-Yaw had done a similar study back in 2008. She embarked on a road trip from Quebec all the way up to Alaska to collect wood frog tissue. So, out first step was to ask Dr. Lee-Yaw if she would collaborate and share her samples.

Those samples gave us a solid backbone across the wood frog range, but we were missing population in expansive regions north and west of the road systems. We worked with the Peabody Museum to search for tissue samples that were already housed in natural history collections around the world. We filled a few gaps, but huge portions of the range were still missing.

We knew that there must be samples out there sitting in freezers and labrooms that were not catalogued in museum databases. So, our next step was to begin sleuthing. We looked up author lists from papers and cold-called leads. I even reached out to friends on Facebook (…which actually turned out to be a big success. The aunt of a friend from undergrad happens to do herpetology research in Galena, Alaska and was able to collect fresh samples for us this year!). This effort greatly expanded our sample coverage with new connections (and friends) from Inuvik and Norman Wells in the Northwest Territories, Churchill on the Hudson Bay, and the Stikine River Delta in Southeast Alaska.

But as the points accumulated on the map, we noticed some glaring holes in our coverage. Most importantly, we had no samples from Northwestern Alaska. Populations in this region are the most distant from the ancestral origin of all wood frogs in the southern Great Lakes. If we wanted a truly “range-wide” representation of wood frog samples, we needed tissue from that blank spot on the map!

Step 2: Convince your advisor and funders it’s a good idea.

This might be the hardest step. In our case, Yara and I were lucky that our advisor, Dave, was immediately supportive of the project. After we made the case for the importance of these samples, funders came around to the idea as well.

Step 3: Make a plan …then remake it …then make a new plan yet again.

Once we knew where we required samples from, we needed to figure out how to get there. Alaska in general is remote, but northwestern Alaska is REALLY remote. The road system doesn’t stretch farther than the middle of the state. All of the communities–mainly small villages–are only accessible by plane, and most of them only have runways for tiny prop planes. Travelling out from the villages into the bush is another layer of difficulty. Most people here either travel by boat on the river or by snowmachine during the winter. Traveling on land, over the soggy and brush-choked permafrost, is brutal and most locals only do it when necessary, if at all.

Prior to academia, I made a career of organizing expeditions to the most remote places in the rugged southeastern archipelago of Alaska. Despite my background, the logistic in the Arctic were even inscrutable to me. Fortunately, I had a couple of friends, Nick Jans and Seth Kantner, who know the area well. In fact, Seth grew up in a cabin out on the Kobuk. (Seth and Nick are both talented authors. I suggest checking out Ordinary Wolves by Seth and The Last Light Breaking by Nick). With their help, I was able to piece together the skeleton of a trip.

After many logistic iterations, Yara and I decided to follow in the footsteps of local hunters who, for generations, have used the rivers as conduits into the heart of the wilderness. Our plan was to travel down one of the major arterial rivers and hike inland to search for frog as we went.

Our original itinerary was to raft the 100 mile section of the Kobuk River from just north of Ambler village to the village of Kiana. But at the last minute (literally), our plans changed. As we were loading up the plane, the pilot told us that he couldn’t fly into our planned starting point. Instead, he suggested that we fly into a gravel bar 30 miles up river in Gate of the Arctic. Those “30 miles” turn out to be AIR MILES. Following the river, it ended up adding over 60 miles to our trip.

We packed two inflatable oar rafts, almost 150 pounds of food, and another 300 pounds of camping, rescue, and science gear, into the balloon-wheeled plane. For the next two weeks, we rowed down the swift Ambler River from the headwaters to the confluence of the Kobuk. Then, we rowed down the massively wide and meandering Kobuk River, eventually extending our trip by an additional 30 miles, by-passing Kiana, and continuing to Noorvik, the last village on the river.

Step 4: Recruit a crew.

Despite being the worlds first and only Saudi Arabian Arctic Ecologist with limited camping experience, I knew Yara would be a stellar field partner. But I never like traveling in brown bear country with fewer than four people. Plus, expedition research involves too many daily chores for the two of us to manage alone. So, we recruited a team.

Sam Jordan is a dry land ecologist, but he had been willing to help me with my dissertation fieldwork in wetlands before, so I knew he would be willing to defect for a good adventure. Sam is also an exceptional whitewater paddler and all-around outdoor guru. Plus, he’s just a great guy (when he leaves his banjo at home). He and I spend two weeks floating the Grand Canyon in the dead of winter and there are few people I would want along on a remote river trip.

Kaylyn Messer and I guided sea kayak expeditions in Southeast Alaska back in our youth. I am a bit particular about how I manage my camp system (read: “extremely picky and fastidious to a fault”) on big trips. Kaylyn is one of the few people as scrupulous as me, but she’s also a super amenable Midwesterner at heart. I knew she’d be a huge help out in the field.

We fell into an effective rhythm on the trip.  Each morning we woke, made breakfast, broke camp, packed the boats, and launched early in the day. While one person on each boat rowed, the other person checked the maps for frog surveying spots, fished, or photographed. We stopped along the way to bushwhack back into wetlands we’d identified from satellite images. We typically arrived at camp late. Yara and I would set up one tent to process the specimens from the day while Same and Kay made camp and cooked dinner. One of the hidden disadvantages of 24-hour Arctic sunlight is that it is easy to overwork. Most nights we only managed to get sampled finished, dinner cleaned up, and camp bearproofed with enough time to crawl into tents with just eight hours till beginning again the next day.

Step 5: Do the science.

Doing science in the field is difficult. Tedious dissections seem impossible while baking in the omnipresent sun and being alternately hounded by hundreds of mosquitoes or blasted by windblown sand. Trading lab coats for rain jackets and benchtops for sleeping pads covered in trashbags compounds the trouble. Not to mention, keeping tissues safe and cool. Organization and adaptability go a long way.

On remote, self-supported trips, it is inevitable that equipment fails or is lost. On one of the first days, we discovered that our formalin jar was leaking—and formalin is not something you want sloshing around! We cleaned the boats and found a creative solution to replace the offending container: a 750ml Jack Daniel’s bottle!

Planning ahead and engineering backup plans also helps. One of our main struggles was figuring out how to preserve specimens and get them home. It is illegal to ship alcohol by mail and you can’t fly with the high-proof alcohol needed for genetic samples. You can ship formalin, but it is difficult to fly with. To make matters worse, we were flying in and out of “dry” or “damp” villages where alcohol is strictly regulated or forbidden. Also, we happened to be flying out on a Sunday, making it impossible to mail samples home. The solution we arrived at was to ship RNAlater and formaldehyde to our hotel room ahead of time. Tissue would remain stable in RNAlater for a couple of weeks and we could make formalin to fix the specimens. After fixing, we cycled the specimens through water to leach out the formalin. This made it possible for me to fly with all of the tissue tubes and damp specimens in my carry on. Other than a few concerned looks from the TSA folks, all of the samples made it back without issue!

Step 6: Enjoy the adventure.

Despite the hard work, there was a lot to appreciate about the Arctic. We witnessed major changes in ecology as we travelled from the steep headwater streams in the mountains to the gigantic Kobuk. Every day was an entirely new scene.

Here’s an illustrated thread on our collecting trip to the Alaskan Arctic for the @Yale @yalepeabody this past summer. 1/22 pic.twitter.com/Jh3VCqVEgx

— A. Z. Andis Arietta (@azandisarietta) August 24, 2021

Step 7: Forget the hardships

Looking back, it is really easy to forget the sweltering heat, swarms of mosquitoes, inescapable sun, and freak lightning storms. And, it’s probably better to forget those anyway!

Arctic CircleSam illustrates that we are, in fact, above the Arctic Circle.

GearGear (minus food and rafts) ready for weigh-in before loading the plane.

Upper Ambler RiverOur drop-off point on the braided upper reach of the Ambler River in Gate of the Arctic National Park.

Gates of the ArcticGetting dropped off on a gravel bar in Gates of the Arctic National Park.

The first dayUnloading the plane on the first day of the trip.

Camp, night 1Campsite along the Kobuk River.

Rigging the raftsPacking up the two-person rafts on the first day.

Kaylyn Messer at the helm.In the upper rivers, near the headwaters, the current was pleasantly swift.

Upper riverSam and Yara rowing down the Kobuk River.

Early morning rowingRowing the swiftwater on the headwaters of the Ambler River in Gates of the Arctic National Park.

Searching for frogs

First frogs!Yara and I after finding the first frogs of the trip!

Science in the fieldYara preparing to do the sciencing with a garbage back as a makeshift lab bench.

Dinner. night 6Dinner and a view on the banks of the Kobuk River.

Midnight sunTrying to avoid getting a sunburn while sleeping in the Arctic.

Midnight sunThe sun baking our tent at 1130pm.

Rainy rowingSam and Yara floating downstream.

Arctic frogYara proudly presents her first wood frog collection.

Kobuk viewsSmoke and clouds mute the midnight sun.

Sam JordanSam Jordan in his happy place.

Chocolate con culicidae.With hundreds of bugs circling your head at all times, it is impossible not to end up with a few in every meal.

CutbanksMassive cutbanks made much of the shore inaccessible for sampling.

StormA freak lightning storm.

SunsetBeautiful light from the midnight sun.

Field dissection Taking tissue samples for sequencing.

TadpolesTadpoles from "Seth's Pond".

Prepping collectionsYara and I taking tissue samples and prepping specimens.

Arctic Wood FrogBeautiful arctic wood frog.

Fixing specimensHere, I am fixing specimens with formalin to preserve them for the Peabody Museum.

Bushwhacking The view while bushwhacking in to search for frogs.

BugsA random sample of friendly mosquitoes.

BearGrizzly bear prints on the beach.

FishingSam with his first of many greyling.

The Kobuk River valleyThe Kobuk Valley watershed is just an endless expanse of wetlands.

The end of the tripFlying back to Kotzebue in a packed plane at the end of the trip.

Heading homeCatching a ride out to the airstrip with our truckload of gear on the last day of the trip.

Rotate around the camera

Perfecting panoramas took some practice for me. The first key is to avoid rotating the camera around your body, but rather, to rotate your body around the camera. I found this post helpful when I was learning. The principle is to keep the entrance pupil of the camera lens as the exact center of the sphere (as much as possible).

Calibrate your compass and gyro

Another tip is to calibrate your compass and gyro regularly (here’s how for Android users). If the compass is off, the phone won’t always know when you’ve made a full rotation.

Take your time

It also helps to go slowly. The camera is constantly processing and mapping the scene to guide the next image, which takes a few seconds.

Be consistent

Even though you can technically take images in any order, I always rotate the same direction around the horizon and then the upper two circles. The lower sphere doesn’t really matter, so I don’t spend much time on it. The key here is that consistency, both for you and the camera, helps.

Use your body as a tripod

You can also try to use your body as a kind of tripod. I used a specific method when I was first learning to take spherical images: I hold my right elbow to my hip bone directly over my right knee and pinch the phone at the top in line with the camera. Then I keep my right foot planted and walk my left foot around the right, using my left hand to tip the phone as need. After a few dozen spheres, I don’t really need to do that any more. I usually know if the sphere will turn out bad and just restart.

Review your images

It is important to remember that you can (and should) review your images in the field so that you can retake any with issues. You can also get a Google Cardboard viewer for $13 which lets you review the sphere in VR in the field to more clearly see the images.

Over time, you get the hang of capturing spheres to minimize problems. It takes me much less time to get a good photosphere than it does to level and set my exposure with my DSLR, now that I have practiced.

Nevertheless, there are always some artifacts. After looking at lots of images side-by-side, my sense is that the artifactual errors in the panos, even when fairly large, contribute less error than the total errors of incorrect orientation, unleveled camera, poor lens calibrate, exposure misspecification, and much lower pixel count from other methods.

Be sure to check out my other posts on this topic:

Smartphone hemispherical photography

Analyzing hemispherical photos

Taking hemispherical canopy photos

Hardware for hemispherical photos

Hemispherical light estimates

The study:

One of the most consistent findings arising from 20 years of study in our lab is that wood frogs seem to adapt to life in cold, dark ponds. In general, cold-blooded animals like reptiles and amphibians are not suited for the cold and function much better in warmer conditions. So, wood frogs that live in colder ponds should have a harder time competing against their neighbors in warmer ponds.

In response, cold-pond wood frogs seem to have developed adaptations that level the playing field. In separate experiments, we’ve found that wood frog tadpoles in cold-ponds tend to seek out warmer water (like in sunflecks) and have lower tolerance to extremely warm temperatures. Most importantly, they can mature faster as eggs and larvae.

But I’ve always struggled with a lingering question: if cold pond frogs have adapted these beneficial adaptation to compete with warm-pond frogs, what is keeping those genes out of the warm-ponds? Shouldn’t cold-pond genes in a warm pond mean double the benefits? One would expect the extrinsic environmental influence and the intrinsically elevated growth rates to produce super tadpoles that metamorphose and leave the ponds long before all the others.

Kaija, who was an undergrad in the lab at the time, decide to tackle that question for her senior thesis.

We hypothesized that there might be a cost to developing too quickly. Studies in fish suggested that the trade-off could be between development and performance. The idea is that, like building Ikea furniture, if you build the tissue of a tadpole too quickly, the price is loss of performance.

So we collected eggs from 10 frog ponds that spanned the gradient from sunny and warm to dark and cold. We split clutches across two incubators that we set to bracket the warmest and coldest of the ponds.

Then we played parents to 400 tadpoles, feeding and changing water in 400 jars two to three times a week.

We reared the tadpoles to an appropriate age (Gosner stage 35ish). Those in the warm incubator developed about 68% faster than those in the cold incubator. In addition to our lab-reared tadpoles, we also captured tadpoles from the same ponds in the wild as a comparison. Development rates in the lab perfectly bounded those in the wild.

Once they reached an appropriate size, we put them to the test. We simulated a predator attack by a dragon fly naiad by poking them in the tale. Dragonfly naiads are fast, fierce, tadpole-eating machines and a tadpole’s fast-twitch flight response is a good indicator of their chance of evading their insect hunters. It’s a measure of performance that directly relates to a tadpole’s fitness.

Above the test arenas, we positioned highspeed cameras to capture the tadpoles’ burst responses. We recorded 1245 trials, to be exact—way more than we ever wanted to track by hand. Fortunately, Kaija is a wiz at coding; and with a bit of help, she was able to write a Matlab script that could identify the centroid of a tadpole and record its position 60 times per second.

We measured the tadpoles’ speed during the first half second of their burst response and looked for an association with their developmental rates. One complicating factor is that a tadpole’s fin and body shape can influence burst speeds. So, a weak tadpole with a giant fin might have a similar burst speed to a super fit tadpole with a small fin. To account for this, we took photos of each tadpole and ran a separate analysis mapping their morphometry and included body shape into our models.

As we had hypothesized, tadpoles reared at warmer temperatures show much slower burst speed than their genetic half-sibling reared in the cold incubator. We even saw a similar, but weaker effect for the tadpoles that were allowed to develop in their natal ponds. It seems that developing too fast reduces performance.

Thus, it certainly seems that the counter-gradient pattern we see of faster development in cold-pond populations, but not in warm-pond populations, is at least partially driven by the trade-off between development rate and performance.

In fact, it may even be the case that we’ve been viewing the pattern backwards all along. Perhaps instead we should consider if warm-pond populations have developed adaptively slower development rates to avoid the performance cost. This especially makes sense given the range of wood frogs. Our populations are at the warm, southern end of the range. Maybe this tradeoff is also a factor constraining wood frogs to the cold north of the continent?

If warm weather and faster development are a real liability for wood frogs, it is only going to get worse in the future. We know from another recent study that our ponds have been warming quickly, especially during the late spring and early summer months. But climate change is also causing snow to fall later in the winter forcing frogs to breed later. The net result is that wood frogs may be forced to develop fast intrinsic developmental rates in response to a contracting developmental window, while at the same time, extrinsic forces drive development even faster. That’s a double whammy in the trade-off with performance. And might lead to too many “Ikea furniture mistakes” at the cellular level.

As a separate part of this study, we also measured metabolic rates in out tadpoles in hopes of understanding the relationship between developmental rates, performance, and cellular respiration. I’m in the process of analyzing those data, so stay tuned for more!

Now, five years, four blog posts, and uncountable hours later, I can say that measuring canopy structure with hemispherical photos is surprisingly difficult.

One of the biggest hurdles is understanding the equipment and deploying it properly in the field. For me, nothing is more tedious than standing waste deep in a mucky, mosquito-infested pond while I fiddle around with camera exposure settings and fine-tuning my leveling device. Add to that the constant fear of dropping a couple thousands of dollars of camera and specialized lens into the water, and you get a good sense of my summers.

So, it is with great pleasure that I offer an alternative method for capturing canopy photos that requires nothing but a cell phone, yet produces higher quality images than a traditional DSLR setup. This new method exploits the spherical panorama function available on most cameras (or in the free Google Street View app). Hemispherical photos can then be easily extracted a remapped from the spheres. You can check out the manuscript at Forestry here (a PDF is available on my Publications page) or continue reading while I walk through the paper below.

The old way

The most common way to measure canopy structure these days is with the use of hemispherical photographs. These images capture the entire canopy and sky from horizon to horizon. Assuming proper exposure, we can categorize individual pixels as either sky or canopy and run simple statistics to count the amount of sky pixels versus canopy or the number and size of the gaps between canopy pixels. We can also plot a sun path onto the image and estimate the amount of direct and indirect light that penetrated through the canopy. (You can follow my analysis pipeline in this post).

All of this analysis relies on good hemispherical images. But the problem is that there are many things that can go wrong when taking canopy photos, including poor lighting conditions, bad exposure settings, improperly oriented camera, etc. Another problem is that capturing images of high-enough quality requires a camera with a large sensor, typically a DSLR, a specialized lens, and a leveling device, which can cost a lot of money. Most labs only have one hemispherical photography setup (if any), which means that we sacrifice the number of photos we can take in favor of high-quality images.

The new way

In the past few years, researchers have tried to figure out ways to get around this equipment barrier. Folks have tried eschewing the leveling device, using clip-on hemispherical lenses for smartphones, or using non-hemispherical smartphone images. I even tested using a hemispherical lens attachment on a GoPro.

But, none of these methods really produce images that are comparable to the images from DSLRs, for three reasons:

1. Smartphone sensors are tiny compared to DSLR, so there is a huge reduction in quality.
2. Clip-on smartphone lenses are tiny compared to DSLR, so again, there is a huge reduction in optical quality.
3. Canopy estimates are sensitive to exposure settings and DLSRs allow for more control over exposure.

The method I developed gets around all of these issues by using multiple, individual cell phone images to stitch together a single hemispherical image. Thus, instead of relying on one tiny cell phone sensor, we are effectively using many tiny cell phone sensors to make up the difference.

Another advantage of creating a hemispherical image out of many images is that each individual image only has to be exposed for a portion of the sky. This avoids the problems of glare and variable sky conditions that plague traditional systems. An added benefit is that, smartphone cameras operate in a completely different way than DSLRs, so they are much less sensitive to exposure issues in general.

Smartphones are less sensitive to exposure issues because, unlike DSLRs that capture a single instance on the sensor when you hit the shutter button, smartphone cameras use computational photography techniques that blend the best parts of many images taken in short succession. You may not realize it, but your smartphone is constantly taking photos as soon as you turn it on (which makes sense since you can see the scene from the camera on your screen). The phone stores about 15 images at a time, constantly dumping the older versions out of temporary memory as updated images pour in. When you hit the button to take a picture, your phone then automatically blends the last few images with the next few images. The phone’s software selects the sharpest pixels with the most even contrast and color from each image and then composites those into the picture presented to you. With every new software update, the algorithms for processing images get better and better. That’s why modern cell phones are able to take photos that can compete with mid-range DSLRs despite the limitations of their tiny sensors.

So, if each phone photos is essentially a composite of 15 images, and then we take 18 of those composite images and stitch them into a hemispherical image, we are effectively comparing a sensor the size of 270 individual phone camera sensors to the DSLR sensor.

The best part is that there is already software that can do this for us via the spherical panorama feature included with most contemporary smartphone cameras. This feature was introduced in the Google Camera app back in 2012 and iOS users can access the feature via the Google StreetView app. It is incredibly simple to use.

Update: Check out my post on tips for taking spherical panoramas

Once you’ve taken a spherical panorama, it is stored in your phone as a 2D JPEG in equirectangular format. The best part about the photo sphere software is that it utilizes your phone’s gyroscope and spatial mapping abilities to automatically level the horizon. This is advantageous for two reasons. First, it means we can ditch the tedious leveling devices. Second, it means that the equirectangular image can be perfectly split between the upper and lower hemisphere. We simply have to crop the top half of the rectangular image and remap it to polar coordinates to get a circular hemispherical image.

How to extract hemispherical images from spherical panoramas

UPDATE: Please see my latest post to process spherical images with R.

Command line instructions

If you are proficient with the command line, the easiest way to extract hemispherical images from photo spheres is to use ImageMagick. After you download and install the program you can run the script below to convert all of your images with just a couple lines of code.

cd "YOUR_IMAGE_DIR"

magick mogrify -level 2%,98% -crop 8704x2176-0-0 -resize "8704x8704!" -virtual-pixel horizontal-tile -background black +distort Polar 0 -flop -flip *jpg

You may need to make a few modifications to the script for your own images. The -crop 8704x2176-0-0 flag crops the top half of the image (i.e. upper hemisphere). Be sure to adjust this to 1.00×0.25 the dimensions of your panorama dimensions. The -resize "8704x8704!" flag resizes the image into a square in order to apply a polar transformation. Be sure to adjust this to 1.00×1.00 the width of your panorama

Note that the code above will convert and overwrite all of the .jpg files in your folder to hemispherical images. I suggest that you practice on a folder of test images or a folder of duplicates to avoid any mishaps.

GUI instructions

If you are intimidated by the command line, extracting hemispherical images from photo spheres is also easy with GIMP (I used GIMP because it is free, but you can follow the same steps in Photoshop).

Update: You can also try out this cool web app developed by researchers in Helsinki which allows you to upload spherical panoramas from your computer or phone and automatically converts them to hemispherical images that you can download. However, I would not suggest using this tool for research purposes because the app fixes the output resolution at 1000p, so you lose all of the benefits of high-resolution spherical images.

First, crop the top half of the rectangular photo sphere.

Second, scale the image into a square. I do this by stretching the image so that the height is the same size as the width. I go into why I do this below.

Third, remap the image to a polar projection. Go to Filter > Distort > Polar Coordinates

Fourth, I found that increasing the contrast slightly helps the binarization algorithms find the correct threshold.

All of these steps can be automated in batch with BIMP plugin (a BIMP recipe is available in the supplemental files of the paper). This can also be automated from the command line with ImageMagick (see scripts above and in the supplemental materials of the paper).

The result is a large image with a diameter equal to the width of the equirectangular sphere. Because we are essentially taking columns of pixels from the rectangular image and mapping them into “wedges” of the circular image, we will always need to down sample pixels toward the center of the circular image. Remember that each step out from the center of the image is the same as each step down the rows of the rectangular image. So, the circumference of every ring of the circular image is generated from a row of pixels that is the width of the rectangular image.

With a bit of geometry, we can see that, the circumference matches the width of our rectangular image (i.e. 1:1 resolution) at zenith 57.3 degrees. Zenith rings below 57.3 will be downsampled and those above will be scaled up and new pixels will be interpolated into the gaps. Conveniently, 57.3 degrees is 1 radian. The area within 1 rad, from zenith 0° to 57°, is important for canopy estimates as gap fraction measurements in this portion of the hemisphere are insensitive to leaf inclination angle, allowing for estimated of leaf area index without accounting for leaf orientation.

Thus, we retain most of our original pixel information within this critical portion of the image, but it does mean that we are expanding the pixels (increasing the resolution) closer to the horizon. I tested the impact of resolution directly in my paper and found almost no difference in canopy estimates; so, it is probably okay to downscale images for ease of processing if high resolution is not needed.

The hemispherical image produced can be now be analyzed in any pipeline used to analyze DSLR hemispherical images. You can see the pipeline I uses in this post.

How do images from smartphone panoramas compare to DSLR

In my paper, I compared hemispherical photos taken with a DSLR against those extracted from a spherical panorama. I took consecutive photos at 72 sites. Overall, I found close concordance between measures of canopy structure (canopy openness) and light transmittance (global site factor) between the methods (R² > 0.9). However, the smartphone images were of much greater clarity and therefore retained more detailed canopy structure that was lost in the DSLR images.

Although the stitching process occasionally produces artifacts in the image, the benefits of this method far outweigh the minor problems. Care when taking the panorama images, as well as ever-improving software will help to minimize imperfect stitching.

Overall, this method is not only a good alternative, it is probably even more accurate than traditional methods because of the greater clarity and robustness to variable exposure. My hope is that this paper will help drive more studies in the use of smartphone spheres for forest research. For instance, 360 horizontal panoramas could be extracted for basal measurement or entire spheres could be used to spatially map tree stands. The lower hemisphere could also be extracted and used to assess understory plant communities or leaf litter composition. Researchers could even enter the sphere with a virtual reality headset in order to identify tree species at their field sites from the comfort of their home.

Mostly, I’m hopeful that the ease of this method will allow more citizen scientists and non-experts to collect data for large-scale projects. After all, this method requires no exposure settings, no additional lenses, and is automatically levelled. The geolocation and compass heading can even be extracted from the image metadata to automatically orient the hemispherical image and set the location parameters in analysis software. Really, anyone with a cell phone can capture research-grade spherical images!

Be sure to check out my other posts about canopy research that cover the theory, hardware, field sampling, and analysis pipelines for hemispherical photos, and my tips for taking spherical panoramas.

Update (Nov. 2020): And, I’m especially thrilled that our piece will be the cover article for the journal, featuring a pair of breeding wood frogs from our population! My hands nearly froze trying to get this underwater shot, so I’m glad it was worth the effort!

Over the past 20+ years, our lab has been monitoring over 50 populations of wood frogs at Yale Myers Forest. Each year in early spring, we listen for the duck-like clucks of the male frogs which means that they have emerged from under the snow and moved into the breeding ponds. Shortly afterward, we head out into the freezing ponds to count the egg masses as a way to monitor population density over time.

In this study, we looked at how the oviposition date (the day on which frogs deposited eggs) has changed over time. As climates warm, we usually expect for the timing life-history events (like oviposition, emergence from hibernation, flowering time, etc.) called ‘phenology’ to advance in the year as winters get shorter. That’s just what most species do. And the trend of advancing phenology is strongest for amphibians.

Given that annual temperatures at our field site have increased by almost 0.6 C in the past two decades, we expected frogs to breed and lay eggs earlier. If our frogs were like other amphibians, we might expect oviposition to come around 6 days earlier.

Surprisingly, we found the opposite. Our frogs seem to be breeding 3 days LATER.

To figure out what might be going on with our frogs, we decided to look more closely at climate across the season, not just annual averages. It turns out that most of the increase in annual temperatures are felt later in the summer, but relatively less when frogs are breeding. Snowpack, on the other hand, is actually accumulating later and lasting longer. In the figure (Fig. 3 from the paper) below, you can see these trends. On the left are the comparisons between temperature, precipitation, and snowpack between 1980 and 2018. On the right, we plot only the difference in trends over time. At the top-right, we plot the oviposition dates to show how seasonal changes in climate line up with frog breeding.

We also looked at how the timing of oviposition correlated with climate across the season. We found that breeding occurs later when there is more snow at the beginning of the breeding window. Also colder temperatures just before breeding correlate with delayed oviposition (which makes sense if colder temps mean more snow and less melting).

So, we think that frogs may be kind of stuck. Persistent snowpack might be keeping them from breeding earlier. But at the same time, warmer summer temperatures might be drying up their ponds faster. If so, this could be a big problem for tadpoles that need to maximize their time for development. The figure below shows that frogs tend to breed earlier when winter and early spring air temperature are high. As we’d expect, more snowpack correlates with later breeding. High precipitation during the spring delays breeding (probably because it is falling as snow).

Twenty years is a long time to be collecting ecological data, but it is a pretty short window into the evolutionary history of wood frogs. And, we don’t know how long snow and temperatures may have been working against these frogs. So, as a final piece of our analysis, we used a machine learning technique called a random forest to predict oviposition dates backwards in time an additional 20 years. It doesn’t seem like much has changed over the past half-century or so. In one way, that could be good news in that at least things don’t seem to be getting any worse.

The big question is, how will frogs cope with these climate changes? If tadpoles are faced with an ever-shrinking window of time to develop into frogs, will they be able to keep up? Or, will they lose the race and end up as tadpole-shaped raisins in our ponds?

I won’t give away any spoilers, but I’m looking at our long-term larval datasets to ask that question next.