R versus python

[mpadge]

To avoid junking up #12, I've created a new document at RehabCatchments/rvspy (in the README so can be directly viewed on github) that attempts to examine reasons for the different results generated by @gboeing in python and myself in R. To repeat the table in #12, these differences were in terms of final estimated case loads on each rehab centre (here just the relative percentages):

Destination	R	python
CaseyHospital	19.4	12.8
DandenongHospital	29.4	37.9
KingstonHospital	51.2	49.4

I had hypothesised there that differences could be due to

• Thanks to input of @richardbeare, the R code uses the PSMA::fetch_postcode() function which generates random samples of actual street addresses within a postcode, where the python code uses a simple sample of network nodes

• I've been uncertain for a while whether osmx.network_from_xxx -> networkx.shortest_path_length actually does the same thing as dodgr_dists(), and suspect in fact not. @gboeing Your wisdom greatly appreciated here, but in my potentially uniformed view, osmx.nework_from_xxx extracts the specified part of the network (here, network_type = "drive"), but does not actually weight that to generate weighted route preferences. The networkx.shortest_path_length calculation is then in linear km (or whatever), but absent a preferential weighting scheme. In contrast, dodgr calculates a dual graph with a weighted version used for preferential routing, and an unweighted version used for distance calculations. That could be the source of some discrepancy here?

The new document here compares both of those differences, yet generates quite robust estimates reflecting the previous R values, yet still failing to recreate the python values. Any insights, help, solutions appreciated here!

[richardbeare]

Snap - I've just been looking at the same thing. The issue is that different sets of postcodes get selected - check the displays in the existing documents. There's a very different set in the vicinity of Casey. In python the 10km radius is to the centroid, while the R version appears to be to the nearest point of the postcode. Either is fine, but lets go with the version that is easier to reproduce on both systems.

[mpadge]

Oh yeah, they really are quite different. I'll incorporate that difference as well and update the code.

---

EDIT: Bingo! That's the origin of the difference. Code push impending ...

[mpadge]

hmm... that commit reduced the load to Casey to the desired 13% or so, but also reduced the load to Dandenong so that Kingston is then around 60%. That means there are still other differences at play here.

[mpadge]

The postcodes in the 2 sets weren't quite exactly matched, so the above commit fixed that. Results still don't agree too strongly, something like

Destination	R	python
CaseyHospital	13.8	12.8
DandenongHospital	31.2	37.9
KingstonHospital	55.0	49.4

I suppose that's not too bad, but it would still be nice to have them somewhat closer ...

[mpadge]

So the python code of @gboeing first estimates the stroke incidence per postcode based on the basic demographic data, which the R code does not do. The latest commit modifies the sampling scheme so that numbers of cases are scaled to the estimate incidence rates per postcode, resulting in ...

Destination	R	python
CaseyHospital	11.5	12.8
DandenongHospital	31.6	37.9
KingstonHospital	56.9	49.4

[mpadge]

That commit gets it to something like:

Destination	R	python
CaseyHospital	15.9 +/- 0.9	12.8
DandenongHospital	37.6 +/- 1.1	37.9
KingstonHospital	46.5 +/- 1.3	49.4

Those values for R come from:

1. Using an unweighted street network retaining only components useable for vehicular routing;
2. Sampling a fixed number of random points directly from this network, and allocating those to postcode areas;
3. Weighting the final case load by per-postcode estimates of stroke incidence.

The error estimates come from the code at the end of the README. Importantly, these error estimates themselves are not very reproducible, indicating that some portion of the R-vs-py differences must be presumed to arise from sampling effects alone. Because the python code was based on 1,000 random points in total, while the R code used that number per postcode (and there are 57 of those), the latter must be presumed more accurate here.

---

In contrast to the above values, a more realistic estimate from the R code is likely to come from using the weighted street network, and sampling addresses within each postcode, which corresponds to these values:

Destination	load
CaseyHospital	13.8
DandenongHospital	39.2
KingstonHospital	46.9

[mpadge]

Future ref for @gboeing: I'll dig more deeply in osmnx-vs-osmdata - something I've been wanting to do for a while. Thanks @richardbeare for the incentive here.