Read in some data

The source data is some Strave-generated GPX files relating to various hikes - I have downloaded several of these into a directory called gpx inside my working directory. The files are:

1. Berenthanti_Ghorepani_Ghandruk_Loop_Hike_Day_1_of_3_.gpx
2. Berenthanti_Ghorepani_Ghandruk_Loop_Hike_Day_2_of_3_.gpx
3. Berenthanti_Ghorepani_Ghandruk_Loop_Hike_Day_3_of_3_.gpx
4. Colorado_Belford_Oxford_and_Missouri_Mountains_Hike.gpx
5. Colorado_Longs_Peak_and_Chasm_Lake_Hike.gpx

The code below reads in all of the files in the gpx directory to give a list of data frames, one for each file. Each data frame is then given an id number, and then all of the data frames are merged. The time variable is transformed to the standard lubridate time format, and the elevation is converted to numeric (for some reason it reads in as character). Finally this data frame is converted into a spatial sf object.

files <- dir('gpx') %>% paste0("gpx/",.)
tracks <- map(files,~{readGPX(.)$tracks[[1]][[1]]}) %>% 
  imap(~ .x %>% mutate(id=.y)) %>%
  bind_rows() %>% 
  mutate(time=ymd_hms(time),ele=as.numeric(ele)) %>%
  st_as_sf(coords=c('lon','lat'),crs=4326)

The Colorado data is in track ids 4 and 5. Pull out track 4 to start with. To prepare for the analysis, set the projection to EPSG 26954 (appropriate for Colorado) and set the time zone to Mountain Time (as a default, GPX files record time as GMT). Finally, extract the coordinates and add them as ordinary numerical columns called x and y - these will be needed to compute the slope and walking rate.

colorado <- tracks %>% filter(id == 4) %>% st_transform(26954)
xy <- st_coordinates(colorado)
colorado <- colorado %>% 
  mutate(x=xy[,1],y=xy[,2],time=with_tz(time,"America/Denver")) 

Now compute the slope and speed for each point on the track:

colorado <- colorado %>% mutate(tlag = as.numeric(as.period(time - lag(time)),'seconds'),
                    dist = sqrt((x - lag(x))^2 + (y - lag(y))^2),
                    rate = 3.6 * dist/tlag,
                    elag = ele - lag(ele),
                    slope = elag/dist)

We can now plot the slope against the walking rate. Here a spline smooth is also fitted. Although there is a lot of scatter around the trend, Tobler’s general idea appears sound - on average the rate peaks a a small downward slope, and falls off gradually for increasing or decreasing slope with respect to this reference slope.

colorado %>% ggplot(aes(x=slope,y=rate)) + geom_point(alpha=0.4) + geom_smooth()

We can also fit parametric curves - here a Tobler-style hiking function of the form

\[ r = a \exp\left( -b\left|S + c\right| \right) \]

is fitted where:

• \(r,S\) are as in the original Hiking function
• \(a,b,c\) are parameters in the hiking function to be estimated.

In Toblers original model, we have \(a=6,b=3.5,c=0.05\) but here we estimate these parameters in terms of least squares fit, via the R function nls (non-linear least squares).

mod1 <- nls(rate ~ a*exp(-b*abs(slope+c)),
            data=as.tibble(colorado),
            start=list(a=3.6,b=3,c=0.2),
            control=nls.control(tol=1.0e-4))
summary(mod1)

## 
## Formula: rate ~ a * exp(-b * abs(slope + c))
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 3.565782   0.044255   80.57   <2e-16 ***
## b 2.032773   0.068779   29.55   <2e-16 ***
## c 0.133212   0.004466   29.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.16 on 3117 degrees of freedom
## 
## Number of iterations to convergence: 10 
## Achieved convergence tolerance: 9.42e-05
##   (1 observation deleted due to missingness)

A curve can then be created and superimposed on the emprical data to see the quality of the fit.

fitcurve1 <- data_frame(slope=seq(-0.75,0.75,l=101)) %>% 
  {mutate(.,rate=predict(mod1,.))}
ggplot(colorado,aes(x=slope,y=rate)) + 
  geom_point(alpha=0.03) + 
  geom_line(data=fitcurve1,col='dodgerblue')

The cusp in this curve is very much an artifact of the model specification. The derivative of the curve, \(\frac{dr}{dS}\) is greatest in magnitude as \(S\) approaches \(-c\), and changes sign when \(S=-c\). Although the non-parametric trend plotted initially is forced to have a continuous first derivative, and so cannot exhibit a cusp, it does show a tendency to have an increasing \(\frac{dr}{dS}\) near to its peak value - suggesting some evidence that a Tobler-style double exponential model is reasonable.

However, an alternative formulation, where \(\frac{dr}{dS}\) flattens out near the point where \(S=-c\) may also be plausible. One such formulation would be to use a Gaussian-like function instead of the double exponential -

\[ r = a\exp\left(-b(S+c)^2\right) \] This is fitted here, again using nls.

mod2 <- nls(rate ~ a*exp(-b*(slope+c)^2),
            data=as.tibble(colorado),
            start=list(a=3,b=2,c=0.2), 
            control=nls.control(tol=1.0e-4))
summary(mod2)

## 
## Formula: rate ~ a * exp(-b * (slope + c)^2)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 3.088581   0.032346   95.48   <2e-16 ***
## b 3.640182   0.202784   17.95   <2e-16 ***
## c 0.168333   0.009563   17.60   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.189 on 3117 degrees of freedom
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 8.357e-05
##   (1 observation deleted due to missingness)

fitcurve2 <- data_frame(slope=seq(-0.75,0.75,l=101)) %>% 
  {mutate(.,rate=predict(mod2,.))}
ggplot(colorado,aes(x=slope,y=rate)) + 
  geom_point(alpha=0.03) + 
  geom_line(data=fitcurve2,col='dodgerblue')

We can look at more sophisticated models, where the falloff in walking speed differs for each side of the optimum slope. For this, we use a more advanced model:

\[ r = a\exp\left(-b\left|S-c\right| - b'\left[(S-c)\right]_+\right) \]

where

\[ [x]_+ = \left\{ \begin{array}{ll} x & \textrm{ if } x > 0 \\ 0& \textrm{ otherwise} \end{array} \right. \] The function is defined in R as follows:

diode <- function(x) ifelse(x>0, x, 0)

The name diode is used as this emulates the action of an idealised electronic component of that name. The effect of this model is similar to the original Tobler form, but the coefficient on the exponential term when \(x>c\) is \(b + b'\) whereas for \(x \le c\) it is just \(b\). When \(b'=0\) the model is the same as the standard Tobler form.

As before, nls may be used to fit the model to the data.

mod3 <- nls(rate ~ a*exp(-abs(b*slope+c)-diode(b1*slope+c)),
            data=as.tibble(colorado),
            start=list(a=1,b=2,b1=0.4,c=0.2),
            control=nls.control(tol=1.0e-4))
summary(mod3)

## 
## Formula: rate ~ a * exp(-abs(b * slope + c) - diode(b1 * slope + c))
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## a   4.18683    0.14100   29.69  < 2e-16 ***
## b   1.84836    0.08448   21.88  < 2e-16 ***
## b1  0.32261    0.08468    3.81 0.000142 ***
## c   0.21181    0.01794   11.81  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.157 on 3116 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 5.503e-05
##   (1 observation deleted due to missingness)

fitcurve3 <- data_frame(slope=seq(-0.75,0.75,l=101)) %>% 
  {mutate(.,rate=predict(mod3,.))}
ggplot(colorado,aes(x=slope,y=rate)) + 
  geom_point(alpha=0.03) + 
  geom_line(data=fitcurve3,col='dodgerblue')

It might be helpful to investigate the slope as an angle (in degrees), instead of a gradient. If \(s\) is the gradient, then and \(\theta\) is the angle of incline in degrees, then \[ \theta = \frac{180}{\pi} \tan^{-1}(s) \] Here, a helper function is created:

atandeg <- function(s) 180*atan(s)/pi

Then, all models may be viewed as functions of angle of slope:

colorado <- mutate(colorado,angle=atandeg(slope))
fitcurve1 <- mutate(fitcurve1,angle=atandeg(slope),model='Abs/Exp')
fitcurve2 <- mutate(fitcurve2,angle=atandeg(slope),model='Gaussian')
fitcurve3 <- mutate(fitcurve3,angle=atandeg(slope),model='Abs/Exp/Asym')
fitcurves <- bind_rows(fitcurve1,fitcurve2,fitcurve3)
ggplot(colorado,aes(x=angle,y=rate)) + 
  geom_point(alpha=0.03) + 
  geom_line(aes(x=angle,y=rate,col=model),data=fitcurves)

Bayesian Analysis via Laplace’s Approximation

Laplace’s approximation is a very useful tool for carrying out Bayesian inference based on standard regression analysis techniques. An overview of the approach, due to Rasmus Bååth is very helpful in explaining the key idea:

“Have you noticed that posterior distributions often are heap shaped and symmetric? What other thing, that statisticians love, is heap shaped and symmetric? The normal distribution of course! Turns out that, under most conditions, the posterior is asymptotically normally distributed as the number of data points goes to infinity . Hopefully we would then not be too off if we approximated the posterior using a (possibly multivariate) normal distribution. Laplace approximation is a method that does exactly this by first locating the mode of the posterior, taking this as the mean of the normal approximation, and then calculating the variance of the normal by “looking at” the curvature of of the posterior at the mode. That was the handwaving, if you want math, check wikipedia.”

It also turns out that, with a non-informative prior distribution, the ‘curvature’ is essentially the variance-covariance function of the estimated model parameters in nls. Thus, if we assume that \(r\) can be approximated as a normally distributed observation, with \[ \textrm{E}(r_i) = a\exp\left(-b\left|S_i-c\right| - b'\left[(S_i-c)\right]_+\right) \] and \[ \textrm{Var}(r_i) = \sigma^2 \]

This makes it possible to similate draws from the posterior distribution of \((a,b,b',c) |\{r_i,S_i, \textrm{ for } i \in 1 \cdots n\}\), via the mvtnorm package. Here, a sample of 1000 draws from this distribution is generated.

post_sample <- mvtnorm::rmvnorm(10000,coef(mod3),vcov(mod3)) %>% as_data_frame

Noting that \(-c\) is the slope at which maximum hiking speed occurs, a posterior estimate of the angle corresponding to this may be useful.

post_sample %>%
  transmute(angle=-atandeg(c)) %>% 
  ggplot(aes(x=angle)) + geom_density(fill='indianred',alpha=0.5)

A posterior mean estimate with a standard deviation are also possible, as are 95% credibility intervals.

post_sample %>% 
  transmute(angle=-atandeg(c)) %>%
  summarise(Mean=mean(angle),
            Sd=sd(angle),
            Lower_CI=quantile(angle,0.025),
            Upper_CI=quantile(angle,0.975))

## # A tibble: 1 x 4
##    Mean    Sd Lower_CI Upper_CI
##   <dbl> <dbl>    <dbl>    <dbl>
## 1 -12.0 0.979    -13.9    -10.1

Thus, the maximum speed is achieved (on average) at a downhill inclination of around \(12{}^{\circ}\).

It is also possible to obtain a range of simulated fitted \(r\) values for a range of given \(S\) gradients.

sim_hfun <- 
  post_sample %>% 
  mutate(slope=list(seq(-0.75,0.75,l=101)),sim=row_number()) %>% 
  unnest() %>%  
  mutate(r = a*exp(-abs(b*slope+c)-diode(b1*slope+c)))

These can then be used for the basis of a series of superimposed fitted curves, based on the posterior distribution of the parameters.

ggplot(sim_hfun,aes(x=slope,y=r,group=sim)) + geom_line(alpha=1/250)

Grouping by the distinct \(S\) values, posterior means and quantiles may be calculated for each distinct slope value, and used in a ribbon plot, as is done frequently using geom_smooth.

hf_ests <- sim_hfun %>% 
  group_by(slope) %>% 
  summarise(r_l=quantile(r,0.025),r_m=mean(r),r_h=quantile(r,0.975))
ggplot(hf_ests,aes(x=slope,y=r_m,ymin=r_l,ymax=r_h)) + 
  geom_line(col='indianred') + 
  geom_ribbon(alpha=0.2,fill='indianred')