2d:null 2e:["$","title",null,{"children":"COVID SuperEmitters and the logNormal distribution | chonghorizons | Observable"}] 33:Tdb4,{ const chart1 = vl.markCircle() .title(`Infectivity as LogNormal Distribution, center=10^${logCenter.toPrecision(2)}, log StdDev=center=10^${logStdDev.toPrecision(2)}`) .data(data_percentiles) .encode( vl.x().fieldQ("p").title('percentile'), vl.y().fieldQ("x").scale({type: 'log'}).title('infectivity [GE/vol]').axis({format: ".2e"}), vl.tooltip(["p","x"]) ) const chart2 = vl.markLine({color: 'firebrick'}) .data(data_percentiles) .encode( vl.x().fieldQ("p").title('percentile'), vl.y().fieldQ("averageInfectivity").scale({type: 'log'}).title(''), vl.tooltip({format: ".2e"}).fieldN("averageInfectivity") ) const chart3 = vl.markLine({color: 'yellow'}) .data(data_percentiles) .encode( vl.x().fieldQ("p").title('percentile'), vl.y().fieldQ("median").scale({type: 'log'}).title(''), vl.tooltip({format: ".2e"}).fieldN("median") ) return md`This is the companion calculator to a "forthcoming" medrvix post on the logNormal distribution of COVID RNA concentrations. The logNormal distribution underlies many viral load distributions. ## SUMMARY This logNormal RNA distribution has implications on risk and infectivity. When calibrated to observed COVID-19 RNA concentrations: ${await vl.layer(chart1, chart2, chart3).width(300).height(200).render()} The infectivity (as proxied by viral concentration) is highly skewed. * The top 5% have an infectivity share of ${data_percentiles.reduce( (acc,row) => { if (row.p<0.95) return acc; return acc+ row.infectivityshare; },0).toLocaleString(undefined, { style: "percent", maximumSignificantDigits: 3 })} * The top 10% have an infectivity share of ${data_percentiles.reduce( (acc,row) => { if (row.p<0.90) return acc; return acc+ row.infectivityshare; },0).toLocaleString(undefined, { style: "percent", maximumSignificantDigits: 3 })} In other words, given a random sample of 100 COVID positive patients, if they are all emitting quanta (per the Wells-Riley model), the top 10% of people are responsible for ${data_percentiles.reduce( (acc,row) => { if (row.p<0.90) return acc; return acc+ row.infectivityshare; },0).toLocaleString(undefined, { style: "percent", maximumSignificantDigits: 3})} of the virus. (!!!) --- Another way to look at it is by comparing the top percentile quanta emissions to the median. * 10% of cases will have ${(log10NormalInstance.percentile(0.9)/(10**logCenter)).toLocaleString(undefined, {maximumSignificantDigits: 2})} times the median quanta. * 5% of cases will have ${(log10NormalInstance.percentile(0.95)/(10**logCenter)).toLocaleString(undefined, {maximumSignificantDigits: 2})} times the median quanta. This can explain cases of: * SuperSpreading even when masks are worn and precuations are taken. Including in hospitals and careful nursing homes. * The [low secondary attack rate in households](https://github.com/chonghorizons/covid-estimators/wiki/Secondary-Attack-Rates-at-home). Masks typically aren't worn at home and infectivity is highest before symptom onset. The logNormal distribution can explain the dispersion: everyone gets it in some households. Nobody else gets it in most households. _{^{NOTE: These numbers auto-update if you change the logNormal parameters below}} ` }34:T463,{ const chart1 = vl.markCircle() .title(`Infectivity as LogNormal Distribution, center=10^${logCenter.toPrecision(2)}, log StdDev=center=10^${logStdDev.toPrecision(2)}`) .data(data_percentiles) .encode( vl.x().fieldQ("p").title('percentile'), vl.y().fieldQ("x").scale({type: 'log'}).title('infectivity [GE/vol]').axis({format: ".2e"}), vl.tooltip(["p","x"]) ) const chart2 = vl.markLine({color: 'firebrick'}) .data(data_percentiles) .encode( vl.x().fieldQ("p").title('percentile'), vl.y().fieldQ("averageInfectivity").scale({type: 'log'}).title(''), vl.tooltip({format: ".2e"}).fieldN("averageInfectivity") ) const chart3 = vl.markLine({color: 'yellow'}) .data(data_percentiles) .encode( vl.x().fieldQ("p").title('percentile'), vl.y().fieldQ("median").scale({type: 'log'}).title(''), vl.tooltip({format: ".2e"}).fieldN("median") ) return vl.layer(chart1, chart2, chart3).render(); } 35:T1014,{ const chart1 = vl.markCircle().title("Infectivty Share vs Percentile") // Make a scatter chart .data(data_percentiles) .encode( vl.x().fieldQ("p").title('percentile').axis({format: ".0%"}), vl.y().fieldQ("infectivityshare").title('infectivityshare'), vl.tooltip(["infectivityshare","x"]) ) return md` ### Infectivity comparisons The previous figure shows the "infectivity" across ${numPts} representative people. The figures below show the infectivity share. The individual infectivity of all the people are added up to get total infectivity. The share is the proportion of the total infectivity. ${ await vl.markCircle().title("Infectivty Share vs Viral Concentration") // Make a scatter chart .data(data_percentiles) .encode( vl.x().fieldQ("x").scale({type: 'log'}).title('RNA concentration [GE/vol]').axis({format: ".2e"}), vl.y().fieldQ("infectivityshare").title('Infectivity Share'), vl.tooltip(["infectivityshare","x"]) ).render() } ${await vl.layer(chart1).render()} The individual infectivity is assumed to be proportional to the concentration of the virus in a PCR sample, as done in [G. Buonanno et al. (Aug 2020)](https://doi.org/10.1016/j.envint.2020.105794), [G. Buonanno et al. (Dec 2020)](https://doi.org/10.1016/j.envint.2020.106112), and in the Jimenez covid risk calculator [link1, google spreadsheet original](http://tinyurl.com/covid-estimator) and [link2, Observable Notebook](https://observablehq.com/@chonghorizons/covid-individual-event-risk-calculator-v0-2) In particular, Buonanno et al write: > In the present paper we have considered an average c_v value of 10^7 copies/mL. after finding a range of > researchers have recently found c_v values in the range 10^3-10^11 copies mL−1 Though 10^7 may represent the median viral concentration and they are careful to state their assumption, this analysis of the logNormal distribution suggests that the median value may underestimate infectivity. ### Adequate protection Specifically, one might choose a set of preventative measures (ventilation, masks, etc) that protects against the median concentration, currently set at ${(10**logCenter).toExponential(2)}. However, this may not protect against quanta in the edge cases. * 10% of cases will have ${(log10NormalInstance.percentile(0.9)/(10**logCenter)).toLocaleString(undefined, {maximumSignificantDigits: 2})} times the median quanta. _{^{concentrations of at least 10^${Math.log10(log10NormalInstance.percentile(0.9)).toFixed(2)} or ${log10NormalInstance.percentile(0.9).toLocaleString(undefined, {maximumSignificantDigits: 3})}}} * 5% of cases will have ${(log10NormalInstance.percentile(0.95)/(10**logCenter)).toLocaleString(undefined, {maximumSignificantDigits: 2})} times the median quanta. _{^{concentrations of at least 10^${Math.log10(log10NormalInstance.percentile(0.95)).toFixed(2)} or ${log10NormalInstance.percentile(0.95).toLocaleString(undefined, {maximumSignificantDigits: 3})}}} * 1% of cases will have ${(log10NormalInstance.percentile(0.99)/(10**logCenter)).toLocaleString(undefined, {maximumSignificantDigits: 2})} times the median quanta. _{^{concentrations of at least 10^${Math.log10(log10NormalInstance.percentile(0.99)).toFixed(2)} or ${log10NormalInstance.percentile(0.99).toLocaleString(undefined, {maximumSignificantDigits: 3}) }}} ### Infectivity Shares The top 10% of people carry the bulk of the infectivity/viral concentration. * The top 5% have an infectivity share of ${data_percentiles.reduce( (acc,row) => { if (row.p<0.95) return acc; return acc+ row.infectivityshare; },0).toLocaleString(undefined, { style: "percent", maximumSignificantDigits: 3 })} * The top 10% have an infectivity share of ${data_percentiles.reduce( (acc,row) => { if (row.p<0.90) return acc; return acc+ row.infectivityshare; },0).toLocaleString(undefined, { style: "percent", maximumSignificantDigits: 3 })} ` }36:T990,md`## Notes ##### LogNormal parameter values The the LogNormal, I use saliva values approximated from the figures in https://www.medrxiv.org/content/10.1101/2021.01.04.21249236v1. The author does not make the leap of translating saliva concentrations to infectivity, but others have as discussed above. ##### LogNormal parameterization The logNormal distribution as [commonly defined](https://en.wikipedia.org/wiki/Log-normal_distribution) uses mu (μ) and sigma (σ) in base e. I have used the log10 version I've seen more commonly in journal articles reporting genetic material concentrations. One has to scale by a factor ln(10). * The log10 center of the distribution of ${logCenter} converts to a mu of ${Math.log(10**logCenter).toFixed(2)}. * The log10 StdDev of the distribution of ${logStdDev} converts to a sigma of ${Math.log(10**logStdDev).toFixed(2)}. ##### True average of logNormal The logNormal's average is actually higher than what is numerically estimated from a 100 point (actually 99 point) distribution. LogNormal may be too skewed right. By using 1st to 99th percentile, I'm essentially truncating the distribution at the 99th percentile. The real mean of the logNormal is ${realMean.toExponential(2)} compared to the ${numPts} point estimate of ${data_percentiles[0].averageInfectivity.toExponential(2)}. ##### Use the actual empirical distribution My recommendation is to use the actual empirical distribution from nasopharyngeal or saliva from COVID tests (NAAT / PCR tests) which will get the tail distribution more correct. However, the logNormal distribution demonstrates the overwhelming share of infectivity that is in the top 5-10% of COVID positive people. This is also seen in https://www.medrxiv.org/content/10.1101/2021.01.04.21249236v1. Note that the COVID PCR test empirical distribution will mix different people at different *stages* of infection. Ideally, one would want to look at a time profile of the RNA concentration over the roughly 15 days of infectivity. And then, one would want to look at the distribution of maximum infectiveness and also infectiveness over time separately. ##### Download numerical array You can download the numerical percentile array for the logNormal to verify these calculations yourself. From a desktop browser (mobile probably doesn't work), click on the \`data_percentiles\` below. To the left, click on the 3 vertical dots. Then select download as CSV or JSON. `37:T736,class Log10NormalHC { constructor(center=1,log10sigma=1,n=100) { this.mu = Math.log(10**center); // Math.log is ln, or log with base e this.sigma = Math.log((10**log10sigma)); this.mean = Math.exp(this.mu+(this.sigma**2)*0.5); // this.n = n; // this.min = 1.0; } pdf(x) { const fact=1.0/(x*this.sigma*Math.sqrt(2.0*Math.PI)); return fact*Math.exp(-0.5*((Math.log(x)-this.mu)/this.sigma)**2.0); } cdf(x) { // https://en.wikipedia.org/wiki/Log-normal_distribution for formula return 0.5 + 0.5 * erf( (Math.log(x) - this.mu) / (2**0.5 * this.sigma) ) } percentile(p) { if (p<=0 || p>=1) return NaN; // https://en.wikipedia.org/wiki/Log-normal_distribution for formula return Math.exp(this.mu + this.sigma *Math.sqrt(2)* erfinv(2*p-1)); } // plot({ xmax = null, xmin=null}={}) { // var data = []; // if (xmin==null || xmin<=0) // xmin=Math.exp(this.mu-4*this.sigma) // if (xmax==null) // xmax=Math.exp(this.mu+4*this.sigma) // let logXmax=Math.log10(xmax); // let logXmin=Math.log10(xmin); // const logdx = (logXmax-logXmin)/this.n; // for (let i=0; i`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":527,"value":"md`### Infectivity / Viral Concentration\nThis plot shows the infectivity of the viral RNA concentration [GE/vol] of 100 people who are COVID positive. This assumes a logNormal density, [calibrated to saliva data](https://twitter.com/linseymarr/status/1348386351654252550?s=20). The parameters can be adjusted below.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":347,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":656,"value":"html`

LEGEND:
\nred is the average infectivity.\nThe average RNA concentration is ${data_percentiles[0].averageInfectivity.toExponential(2)} GE copies/vol (red line).\n
\nyellow is the median infectivity. This is 10^${logCenter} GE copies/vol, which equals ${(10**logCenter).toExponential(2)}
\nblue circles are the percentiles\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":9,"value":"viewof logCenter = sliderBoxButtons({title: \"log center of distribution\", description: \"The default of 4 means 10^4 copies/vol. 10^4 GE/ul was observed in saliva data. Note that 10^4 GE/ul = 10^7 GE/ml. See below for calibration data source.\", min: 0.1, sliderMax: 10, max: 100, value: 4.0, step: 0.1})\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":23,"value":"viewof logStdDev = sliderBoxButtons({title: \"log stdDev of distribution\", description: `The default of 1.7 is the \\\"stdDev of the log10\\'d concentration\\\" In other words, 68% of the distribution lies between 10^5.7 and 10^2.3 \"centered\" around a median at 10^4. This was observed in saliva data. See below.`, min: 0.1, sliderMax: 7, max: 100, value: 1.7, step: 0.1})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":599,"value":"md`##### Instructions\n\nUse the sliders or boxes to enter your parameters for the logNormal distribution.\nAll charts *and numbers in the text* will autoupdate.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":449,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":484,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":412,"value":"data_percentiles = {\n let a = log10NormalInstance\n let mostdata= rangeN(numPts-1).map(p => p/numPts).map((p)=> {\n let x=a.percentile(p)\n return {p: p, x: x}\n })\n const AccumulateX = (accumulator, row) => accumulator + row.x;\n let totalInfectivity=mostdata.reduce(AccumulateX,0);\n return mostdata.map( ({p, x}) => {\n return {p, x, infectivityshare: x/totalInfectivity, log10x: Math.log10(x), totalInfectivity: totalInfectivity, averageInfectivity: totalInfectivity/mostdata.length, median: Math.exp(a.mu)}; \n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":429,"value":"viewof numPts=number( {value:100, step:1, precision:0, title: \"Number of points for the simulation.\", max:10000, description: \"100 means percentiles 1st to 99th. Max value is 10000 to prevent long calcution time.\", submit: true});\n\n// You can remove the max to run larger values. But generating and array of 100000 points and plotting 100000 points is slow.","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":14,"value":"html`

Imports and code.

\n\n

You can ignore what is below here unless you want to edit the underlying code. The code below is for ObservableHQ to work.\n

\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":697,"value":"md`### ChangeLog\n\nJan 27, 2021: Initial publication.\n\n---\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":581,"value":"log10NormalInstance = new Log10NormalHC(logCenter,logStdDev)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":5,"value":"html`

 No longer using import { LogNormal } from \"@sscollis/distributions\"\n-- but thank you. I based my code off yours.

`\n\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":276,"value":"// tenpts= rangeN(9).map(x=>data4_2percentiles[x*10-1])","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":282,"value":"{ \n // check my math\n // let row=data4_2percentiles[9*10-1]\n // let ninetyShouldbe1pt28x2=row.log10x-4\n // let multipleof128=ninetyShouldbe1pt28x2/1.28\n // // works. 90th percentile is 2 log10 stadard deviations away.\n // return { row, ninetyShouldbe1pt28x2, multipleof128 }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":193,"value":"{\n \n // data4_2percentiles = {\n // let a = new Log10NormalHC(4,2)\n // let mostdata= rangeN(99).map(p => p/100).map((p)=> {\n // let x=a.percentile(p)\n // return {p: p, x: x}\n // })\n // const AccumulateX = (accumulator, row) => accumulator + row.x;\n // let totalInfectivity=mostdata.reduce(AccumulateX,0);\n // return mostdata.map( ({p, x}) => {\n // return {p, x, infectivityshare: x/totalInfectivity, log10x: Math.log10(x), totalInfectivity: totalInfectivity, averageInfectivity: totalInfectivity/mostdata.length, median: Math.exp(a.mu)}; \n // });\n // }\n\nreturn \"OLD: Testing data locked to the mu=4 and sigma=2 distribution\";\n\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":232,"value":"{\n return \"commented\" \n \n let a = new Log10NormalHC(1,1)\n return rangeN(9).concat([9.9]).map((p)=> {\n let x=a.percentile(p/10)\n return {percentile: p/10, x: x, log10x: Math.log10(x)}\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":248,"value":"{ return \"commented\" \n // test erfinv\n return rangeN(9).map((x)=> {\n let p = x/10\n return {percentile: p, x: Math.sqrt(2) * erfinv(2*p-1)}\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":212,"value":"rangeN = (n) => Array(n).fill(1).map((x, y) => x + y)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":467,"value":"realMean = { \n let a= new Log10NormalHC(logCenter,logStdDev)\n return a.mean \n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":62,"value":"$37","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":177,"value":"import { erf, erfinv, plot as plotBostock} from \"@mbostock/error-function\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":10,"value":"import {number, slider} from \"@severo/inputs-setter\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":16,"value":"import { sliderBoxButtons } from \"@chonghorizons/slider-with-number-input-box-and-quick-select-buttons\" ","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":474,"value":"import {number as numberWithSubmit} from \"@jashkenas/inputs\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":57,"value":"import { vl } from \"@vega/vega-lite-api\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":101,"value":"d3 = require(\"d3@6\")\n","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":105,"value":"vegalite = require('@observablehq/vega-lite')","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W38","stdlib":"1","titling":"auto"}} 32:"$31:body"