31:T945,md`# Comparing two methods for determining an efficient frontier In this notebook, I explore two methods for computing an [efficient frontier](https://en.wikipedia.org/wiki/Efficient_frontier). The traditional method I use is as follows: 1. Compute a [convex hull](https://en.wikipedia.org/wiki/Convex_hull) from the points with [d3-polygonHull](https://github.com/d3/d3-polygon#polygonHull)¹ 2. Locate the point on the hull with the smallest cost 3. Follow the points on the hull until we reach the highest benefit In their [paper](https://www.ncbi.nlm.nih.gov/pubmed/25926282) entitled _An Efficient, Noniterative Method of Identifying the Cost-Effectiveness Frontier_ authors Sze-chuan Suen and Jeremy D. Goldhaber-Fiebert claim that their method is 'is conceptually simple, easier to implement and potentially faster for situations that challenge the conventional approach.' They refer to a process that uses the 'net monetary benefit (NMB) defined as a strategy’s benefits, converted to monetary units for a given WTP, minus its costs' (p2 of the article). I've loosely termed their frontier method as the NMB method. I was curious how much of an improvement we'd see with this newer method. I implemented a javascript version of both methods and compared the run time. At a first glance it seems that the traditional hull-based method is actually quite faster. On this browser it was about **${format( NMB_duration / hull_duration )} faster** than the NMB method. We could optimize the NMB implementation somewhat, but not enough to compete with the hull-based method. It's quite possible I'm not fully grasping the nuances here (and welcome comments), but I'm guessing that in this case the iterative approach would seem more efficient than the so-called single pass method. Duration for ${iterations} iterations: | Hull method | NMB method | Ratio | | -------------- | -------------- | -------------- | | ${format(hull_duration)} ms| ${format(NMB_duration)} ms| ${format( NMB_duration / hull_duration )} | ¹Mike Bostock's [d3-polygonHull](https://github.com/d3/d3-polygon#polygonHull) uses A. M. Andrew's [monotone chain method](https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain) I plot two versions of the same frontier for a visual inspection of the results.`32:T4ad,// NMB method to determine frontier roughly adapted from Suen and Goldhaber-Fiebert's article // Admitly this may not be an optimal implementation // I've preserved the variable names to match the original R source code for easy comparison function frontierByNMB(data, maxWTP) { const CEmat = data; let WTPtestPoints = [0]; let ICERmat = CEmat.forEach((d, i) => CEmat.forEach((e, j) => { if (j >= i) return; // avoid redundant calcs let wtp = (d.cost - e.cost) / (d.benefit - e.benefit); if ( wtp > 0 && isFinite(wtp) && WTPtestPoints.indexOf(wtp) < 0 && // unique wtp <= maxWTP ) { WTPtestPoints.push(wtp); } }) ); WTPtestPoints = WTPtestPoints.sort((a, b) => a - b); let NMBmat = []; CEmat.forEach(d => NMBmat.push(WTPtestPoints.map(e => e * d.benefit - d.cost)) ); let NMBmatTransposed = d3.transpose(NMBmat); let maxVals = NMBmatTransposed.map(arr => d3.max(arr)); let indiciesOfMax = []; NMBmatTransposed.forEach((arr, i) => { arr.forEach((d, j) => { if (d == maxVals[i] && indiciesOfMax.indexOf(j) < 0) indiciesOfMax.push(j); }); }); return indiciesOfMax; }33:T525,// convex hull method to determine frontier function frontierByHull(data, maxWTP) { const points = data.map((d, i) => [d.cost, d.benefit, d.index]); const initial_hull = d3.polygonHull(points); // rehydrate point info into an object let hull = initial_hull.map((d, i) => { return _.assign({ cost: d[0], benefit: d[1], original_index: d[2], hull_index: i, on_frontier: false }); }); // Slice the list of points so that we start with the point with the smallest cost // e.g. if array = [1,2,3,4,5,6] and min_cost_index is 3 --> [3,4,5,6,1,2] const min_cost_index = _.minBy(hull, 'cost').hull_index; hull = hull.slice(min_cost_index).concat(hull.slice(0, min_cost_index)); // Iterate through the points on the hull in order // Stop if wtp is > maxWTP or if benefit is less than benefit of previous point let previous_benefit = -Infinity; let previous_cost = 0; hull.forEach(d => { let wtp = (d.cost - previous_cost) / (d.benefit - previous_benefit); if (d.benefit >= previous_benefit && wtp <= maxWTP) { d.on_frontier = true; previous_cost = d.cost; previous_benefit = d.benefit; } }); // return original point index for each point on the frontier return hull.filter(d => d.on_frontier).map(d => d.original_index); }34:[] 2f:{"title":"Comparing two methods for determining an efficient frontier","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":235,"value":"{\n const line1 = vl\n .markLine()\n .data(data.filter(d => _.includes(frontierIndicesNMB, d.index)))\n .encode(vl.x().fieldQ('benefit'), vl.y().fieldQ('cost'))\n .title('Efficient frontier (cost/benefit)');\n\n const point1 = vl\n .markCircle()\n .data(data)\n .encode(\n vl.x().fieldQ('benefit'),\n vl.y().fieldQ('cost'),\n vl.tooltip(['name', 'cost', 'benefit'])\n );\n\n const line2 = vl\n .markLine()\n .data(data.filter(d => _.includes(frontierIndicesNMB, d.index)))\n .encode(vl.y().fieldQ('benefit'), vl.x().fieldQ('cost'))\n .title('Efficient frontier (benefit/cost)');\n\n const point2 = vl\n .markCircle()\n .data(data)\n .encode(\n vl.y().fieldQ('benefit'),\n vl.x().fieldQ('cost'),\n vl.tooltip(['name', 'cost', 'benefit'])\n );\n\n return vl.hconcat(vl.layer(line1, point1), vl.layer(line2, point2)).render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":906,"value":"viewof maxWTP = slider({\n min: 0,\n max: 210000,\n step: 500,\n value: 210000,\n format: \",\",\n description: \"Maximum willingness to pay (WTP)\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":351,"value":"frontierIndicesNMB = frontierByNMB(data, maxWTP)","pinned":false,"hidden":false,"mode":"ojs"},{"id":787,"value":"frontierIndicesHull = frontierByHull(data, maxWTP)","pinned":false,"hidden":false,"mode":"ojs"},{"id":693,"value":"iterations = 100","pinned":false,"hidden":false,"mode":"ojs"},{"id":2,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":366,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":612,"value":"NMB_duration = {\n let time = Date.now();\n for (let i = 0; i < iterations; ++i) frontierByNMB(data, 210000);\n return Date.now() - time;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":575,"value":"hull_duration = {\n let time = Date.now();\n for (let i = 0; i < iterations; ++i) frontierByHull(data, 210000);\n return Date.now() - time;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":9,"value":"// Dataset from Suen and Goldhaber-Fiebert's paper\ndata = {\n const data = await d3.tsv(\n 'https://gist.githubusercontent.com/john-clarke/b01c28eea5f264570eff9c3b6ca6b195/raw/68db83369612d51be45305055b50e13a24e337cd/CostQualyInput.csv',\n (d, index) => ({\n index,\n name: +d.Strategy,\n cost: +d.Cost,\n benefit: +d.QALYs\n })\n );\n // Remove duplicates\n return _.uniqWith(data, _.isEqual);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":711,"value":"format = d3.format(',.0f')","pinned":false,"hidden":false,"mode":"ojs"},{"id":236,"value":"import { vl } from '@vega/vega-lite-api'","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"d3 = require('d3')","pinned":false,"hidden":false,"mode":"ojs"},{"id":807,"value":"_ = require('lodash')","pinned":false,"hidden":false,"mode":"ojs"},{"id":908,"value":"import { slider } from \"@tmcw/inputs/2\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}