Calculate cycling uptake for 'Go Dutch' scenario : uptake_pct_godutch()

[Yuhei110]

I try to calculate cycling uptake for 'Go Dutch' scenario with sample values. The function, uptake_pct_godutch() and the formula give different results.

The formula that I used is:

Equation 1B:
logit(pcycle) = Equation 1A + Dutch parameters = -3.959 + (-0.5963 * distance) + (1.866 * distancesqrt) + (0.008050 * distancesq) + (-0.2710 * gradient) + (0.009394 * distance*gradient) + (-0.05135 * distancesqrt *gradient) + (2.523 *dutch) + (-0.07626 * dutch * distance)

which is obtained from:
APPENDIX 1 : The National Propensity to Cycle Tool for England ( PCT-England ), version 1 : Description of methodology for scenario-building (https://www.jtlu.org/index.php/jtlu/article/downloadSuppFile/862/360)
Thank you in advance, @Robinlovelace

Distance = 15
Gradient = 2

logit =
    -3.959 + 
  (-0.5963 * Distance) + 
  (1.866 * sqrt(Distance)) + 
  (0.008050 * Distance^2) +
  (-0.2710 * Gradient) + 
  (0.009394 * Distance*Gradient) + 
  (-0.05135 * sqrt(Distance) *Gradient) +
  2.523 +
  (-0.07626*Distance)
# Result: -3.144098

pcycle = exp(logit) / (1 + exp(logit))
# Result: 0.04132445

uptake_pct_godutch(Distance, Gradient)
# Result: 0.03159427

[Robinlovelace]

Many thanks for posting this @Yuhei110 I will take a look next week.

[Yuhei110]

I'm just writing to quickly follow up on the above question. Any clarification you could provide would be greatly appreciated. Thank you.

[Robinlovelace]

Hi @Yuhei110 I'm afraid I've not had a chance to look into it but will aim to do so next week. In the meantime if you make any progress let me know. See here for the recently updated manual that contains the up-to-date formula for cycling potential: https://github.com/npct/pct-shiny/blob/a59ebd1619af4400eeb7ffb2a8ecdd8ce4c3753d/regions_www/www/static/03a_manual/pct-bike-eng-user-manual-c1.pdf

[Yuhei110]

Thank you for your reply @Robinlovelace, I will look into the pdf.

[Robinlovelace]

Hi @Yuhei110 I can reproduce your issue for the Government Target. Does this make sense to you, can you see any issues with my code and have you made any further progress on this matter? FYI the parameter values changed from 2016 when the PCT paper was published, with parameters that are in the package, to 2020 as documented in that manual and in the reproducible code below. Heads-up also @mem48 👍

# Aim: test the pct pkg can reproduce data on www.pct.bike

# see https://github.com/ITSLeeds/pct/issues/65#issuecomment-674237160

# Uptake formula from
# https://raw.githubusercontent.com/npct/pct-shiny/a59ebd1619af4400eeb7ffb2a8ecdd8ce4c3753d/regions_www/www/static/03a_manual/pct-bike-eng-user-manual-c1.pdf
#
# logit (pcycle)= -4.018 +  (-0.6369 *  distance)  +
#   (1.988  * distancesqrt)  +  (0.008775* distancesq) +
#   (-0.2555* gradient) + (0.02006* distance*gradient) +
#   (-0.1234* distancesqrt*gradient)
library(pct)

l_pct_2020 = get_pct_lines(region = "isle-of-wight")
pcycle_pct_2019 = uptake_pct_govtarget(distance = l_pct_2020$rf_dist_km, gradient = l_pct_2020$rf_avslope_perc)
summary(pcycle_pct_2019)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> 0.004712 0.007869 0.011871 0.014671 0.020958 0.043210
# Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
# 0.004712 0.007869 0.011871 0.014671 0.020958 0.043210
govtarget_slc_2019 = round(pcycle_pct_2019 * l_pct_2020$all + l_pct_2020$bicycle)

cor(l_pct_2020$govtarget_slc, govtarget_slc_2019) # 0.9953986
#> [1] 0.9953986
round(1 - mean(govtarget_slc_2019) / mean(l_pct_2020$govtarget_slc), digits = 5) * 100
#> [1] 4.774
# [1] 4.774 # values in pct package 2019 are 5% lower

# previous (2019) values:
alpha = -3.959
d1 = -0.5963
d2 = 1.866
d3 = 0.008050
h1 = -0.2710
i1 = 0.009394
i2 = -0.05135

# New (2020) values:
alpha = -4.018
d1 = -0.6369
d2 = 1.988
d3 = 0.008775
h1 = -0.2555
i1 = 0.02006
i2 = -0.1234

# try to reproduce 2020 values:
pcycle_pct_2020 = uptake_pct_govtarget(
  distance = l_pct_2020$rf_dist_km,
  gradient = l_pct_2020$rf_avslope_perc,
  alpha = alpha,
  d1 = d1,
  d2 = d2,
  d3 = d3,
  h1 = h1,
  i1 = i1,
  i2 = i2
)

govtarget_slc_2020 = round(pcycle_pct_2020 * l_pct_2020$all + l_pct_2020$bicycle)
cor(l_pct_2020$govtarget_slc, govtarget_slc_2020) # 0.9953986
#> [1] 0.9947073
round(1 - mean(govtarget_slc_2020) / mean(l_pct_2020$govtarget_slc), digits = 5) * 100
#> [1] 11.599

^{Created on 2020-08-19 by the reprex package (v0.3.0)}

[Robinlovelace]

Update, I think this is close to a solution in terms of reproducibility, could be useful for the dissertation @Yuhei110 👍

l_pct_2020 = pct::get_pct_lines(region = "isle-of-wight")
l_pct_2020$rf_avslope_perc
#>   [1] 0.33 0.34 0.35 0.63 0.63 0.47 1.45 0.37 0.63 1.43 0.38 0.65 0.67 0.51 1.26
#>  [16] 1.49 3.63 3.61 2.72 2.26 2.64 2.81 1.83 1.79 1.56 2.82 2.57 2.54 2.11 3.53
#>  [31] 2.88 2.48 2.81 2.97 2.19 2.29 1.91 2.84 2.73 2.69 2.27 3.60 3.35 3.55 3.72
#>  [46] 1.85 2.09 3.51 1.83 2.42 2.43 2.39 2.35 2.15 3.78 2.36 3.82 3.25 3.25 3.55
#>  [61] 3.04 3.11 2.66 2.64 2.57 2.79 3.40 3.87 4.15 2.86 2.94 3.60 2.60 2.89 3.05
#>  [76] 3.47 2.94 2.80 3.49 3.94 3.28 3.27 3.75 3.05 3.12 2.41 2.45 2.39 3.00 3.37
#>  [91] 3.85 3.44 3.63 3.23 3.23 2.62 2.60 2.53 2.81 3.42 3.29 3.47 1.96 2.70 2.89
#> [106] 2.97 2.92 2.80 2.34 3.10 3.50 2.10 2.69 2.88 3.08 3.03 2.91 2.25 2.81 3.59
#> [121] 2.53 2.54 2.47 3.34 2.79 3.19 2.07 2.57 2.84 2.61 3.10 2.52 1.01 2.80 3.06
#> [136] 2.53 2.31 2.62 2.96 1.96 2.67 3.13 3.27 3.80 4.35
l_pct_2020$rf_dist_km
#>   [1] 15.00 12.40 11.50 15.40 12.80 12.50 18.60 11.70 14.60 19.80 11.20 14.00
#>  [13] 12.50 12.20 19.60 18.30  1.10  1.94 18.00 13.00 17.10 18.30  7.48  7.26
#>  [25]  7.89 12.90 18.80 18.80 17.40  3.06 17.70 12.70 16.80 18.00  7.16  6.81
#>  [37]  7.71 11.90 18.50 18.50 17.10 11.10  6.14 10.20 11.40  6.75  8.13 18.90
#>  [49]  8.49 14.70 16.90 16.90 18.80 17.60  5.09  1.82  1.07 10.70 12.10  7.57
#>  [61] 12.50 18.70  9.15 12.20 12.90 19.50 18.40  4.39  5.59  5.93  7.32 13.10
#>  [73]  7.68 13.90 12.00 13.30 14.50 15.40 17.70  1.98 10.10 11.50  8.54 11.80
#>  [85] 18.00  9.12 12.20 12.80 18.60 17.90 12.30 12.60  6.97 13.00 19.20  7.95
#>  [97] 11.00 11.70 19.40 18.40  1.85 17.10  2.35 19.50  8.55 11.60 11.60 13.50
#> [109] 10.90 16.00 18.70  1.62 18.30  7.33 13.30 13.30 15.10 11.70 17.10 18.50
#> [121]  8.30 11.30 12.00 19.20 17.60  6.65 14.10 13.00 14.90 10.20 15.60 11.00
#> [133] 11.50 19.20 16.20  3.12  3.81 13.00 11.50  2.04 10.40  8.83  9.84  7.93
#> [145]  7.38
uptake_pct_govtarget_2020 = function(
  distance,
  gradient,
  alpha = -4.018,
  d1 = -0.6369,
  d2 = 1.988,
  d3 = 0.008775,
  h1 = -0.2555,
  i1 = 0.02006,
  i2 = -0.1234
) {
  ned_rf_avslope_perc = gradient - 0.78
  distancesqrt = sqrt(distance)
  distancesq = distance^2
  logit_pcycle = alpha + (d1 *  distance)  +
    (d2 * distancesqrt)  +  (d3 * distancesq) +
    (h1 * ned_rf_avslope_perc) + (i1 * distance*ned_rf_avslope_perc) +
    (i2* distancesqrt*ned_rf_avslope_perc)
  boot::inv.logit(logit_pcycle)
}


pcycle_pct_2020 = uptake_pct_govtarget_2020(
  distance = l_pct_2020$rf_dist_km,
  gradient = l_pct_2020$rf_avslope_perc
  )


govtarget_slc_2020 = pcycle_pct_2020 * l_pct_2020$all + l_pct_2020$bicycle
plot(l_pct_2020$govtarget_slc, govtarget_slc_2020)

cor(l_pct_2020$govtarget_slc, govtarget_slc_2020) # 0.9928463
#> [1] 0.993046
round(1 - mean(govtarget_slc_2020) / mean(l_pct_2020$govtarget_slc), digits = 5) * 100
#> [1] 0.254
# New estimate is less than 0.3% out

^{Created on 2020-08-19 by the reprex package (v0.3.0)}

[Robinlovelace]

Hi @Yuhei110 I think I've solved the problem. Can you reproduce the code below?

# install updated version of package:
remotes::install_github("itsleeds/pct")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'pct' from a github remote, the SHA1 (f4484521) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(pct)
# https://www.jtlu.org/index.php/jtlu/article/downloadSuppFile/862/360)
# Equation 1B:
distance = 15
gradient = 2
logit = -3.959 + 2.523 +
 ((-0.5963 - 0.07626) * distance) +
 (1.866 * sqrt(distance)) +
 (0.008050 * distance^2) +
 (-0.2710 * gradient) +
 (0.009394 * distance*gradient) +
 (-0.05135 * sqrt(distance) *gradient)
logit
#> [1] -3.144098
# Result: -3.144098
pcycle = exp(logit) / (1 + exp(logit))
# Result: 0.04132445
boot::inv.logit(logit)
#> [1] 0.04132445
uptake_pct_godutch(distance, gradient, alpha = -3.959 + 2.523, d1 = -0.5963 - 0.07626,
d2 = 1.866, d3 = 0.008050, h1 = -0.2710, i1 = 0.009394, i2 = -0.05135
)
#> [1] 0.04132445
# these are the default values
uptake_pct_godutch(distance, gradient)
#> [1] 0.04132445
l = routes_fast_leeds
pcycle_scenario = uptake_pct_godutch(l$length, l$av_incline)
#> Distance assumed in m, switching to km
plot(l$length, pcycle_scenario)

^{Created on 2020-08-21 by the reprex package (v0.3.0)}

I am very grateful to you for testing this issue, it has led to an improvement in the package and an increase in the level of reproducibility of the PCT. If you find any other issues or have any other comments/suggestions, don't hesitate to ask or even put in a Pull Request if you want to contribute directly to the source code of the package.

[Robinlovelace]

For purposes of testing...