31:T4d0,html`

LRT Application:
Signal Detection from FDA Adverse Event Reporting System (FAERS)

This is a demo of querying data from FAERS database and identifying signals from adverse evenets by using Likelihood Ratio Test (LRT) method.

Part 2: Signal Detection and Visualization

Several statistical methods that are available in the literature to analyze postmarket safety databases, such as the U.S. Federal Drug Administration’s (FDA) adverse event reporting system (FAERS), for identifying drug-event combinations with disproportionately high frequencies, are subject to high false discovery rates. Likelihood ratio test (LRT) based method, via an extensive simulation study, retains good power and sensitivity for identifying signals, as well as controls both the Type I error and false discovery rates.

`32:Ta4b,### **LRT Method Brief Review**
(Please visit basic LRT method introduction if the method is new for you.)
From LRT Method Application Part 1, we successfully build the raw data for LRT analysis. The dataset contains 4 variables: | Variable Name | Definition | |:-----------:|:-------------:| | ntrt (${tex`n_{ij}`}) | number of events for the ${tex`ith`} AE and ${tex`jth`} drug | | nidot | total number of events for ${tex`ith`} AE (marginal total for ${tex`ith`} row). | |ndotj | total number of events for ${tex`jth`} drug (the marginal total for ${tex`jth`} column) | | ndotdot | grand total |
We define the null hypothesis, ${tex`𝐻_0: 𝑝_𝑖=𝑞_𝑖=𝑝_0`} for all 𝑖; and the alternative hypothesis (one-sided), ${tex`𝐻_a: 𝑝_i > 𝑞_𝑖`} for at least one 𝑖 among the 𝐼 AEs. Under the null hypothesis, the maximum likelihood estimates for ${tex`𝑝_0`} is ${tex`\hat{p_0}= 𝑛_{.𝑗} × 𝑛_{..}`}, and the expected number of cases for 𝑖𝑡ℎ AE and 𝑗∗𝑡ℎ drug is ${tex` E_{ij*} = n_{i.} \times \frac{n_{.j*}}{n_{..}} = \frac{n_{i.} \times n_{.j*}}{n_{..}} `}

First, we calculate Log(LR) from observed data.

${tex` \large LR_{ij*}=\frac{L_a(\hat{p_i}, \hat{q_i})}{L_0(\hat{p_0})}= \frac{ (\frac{n_{ij*}}{n_{i.}})^{n_{ij*}} \times (\frac{n_{.j*}-n_{ij*}}{n_{..}-n_{i.}})^{(n_{.j*}-n_{ij*})} } { (\frac{n_{.j*}}{n_{..}})^{n_{.j*}} }, \small i=1,\cdots,I. `} The distribution of 𝑀𝐿𝑅=𝑚𝑎𝑥(𝐿𝑜𝑔𝐿𝑅) under ${tex`𝐻_0`}, is not analytically tractable and is obtained using Monte Carlo simulation. First, the number of cases for each AE, for a given 𝑗∗𝑡ℎ drug, are simulated under ${tex`𝐻_0`}. The joint distribution of ${tex`(𝑛_{1𝑗∗},...,𝑛_{𝐼𝑗∗})`}, conditioning on ${tex`𝑛_{.𝑗∗} `} and ${tex`(𝑛_{1.},...,𝑛_{𝐼.})`} is

${tex`\left( n_{1j}, ..., n_{Ij} \right ) | n_{.j} \sim multinomial \left( n_{.j}, \left( \frac{n_{1.}}{n_{..}}, ..., \frac{n_{I.}}{n_{..}} \right) \right) `}
After running the simulation ${numberOfSimulation} times, Log(LR) and MaxLRs can be calculated for each column of the simulated data.

Threshold at 𝛼 = 0.05 is 95 percentile of the number of simulated 𝑀𝑎𝑥𝐿𝑅𝑠.
Reject ${tex`𝐻_0`} if observed 𝑀𝑎𝑥𝐿𝑅 > threshold. 33:T89f,funcLRT = py`import numpy as np import pandas as pd def flogLRv(x, fnidot, fndotj, fndotdot): logLRvr = x * (np.log(x) - np.log(fnidot)) logLRvr = logLRvr + (fndotj - x) * (np.log(fndotj - x) - np.log(fndotdot - fnidot)) logLRvr = logLRvr - fndotj * (np.log(fndotj) - np.log(fndotdot)); logLRvr.fillna(0) pj = x / fnidot qj = (fndotj - x) / (fndotdot - fnidot) logLRvr[ pj < qj ] = 0 return logLRvr def fLRT(qnum, number, alpha = 0.05, alphaSpending = 1): mylist = qnum mylist["ntrt"] = pd.to_numeric(mylist["ntrt"]) mylist["nidot"] = pd.to_numeric(mylist["nidot"]) mylist["ndotj"] = pd.to_numeric(mylist["ndotj"]) mylist["ndotdot"] = pd.to_numeric(mylist["ndotdot"]) mylist["RR"] = (mylist["ntrt"])*((mylist["ndotdot"]-mylist["nidot"])/(mylist["nidot"]*(mylist["ndotj"]-mylist["ntrt"]))) mylist["stdRR"] = 1/(mylist["ntrt"]) - 1/mylist["nidot"] + 1/(mylist["ndotj"]-mylist["ntrt"]) - 1/(mylist["ndotdot"]-mylist["nidot"]) mylist["stdRR"] = np.sqrt(mylist["stdRR"]) mylist["lowRR"] = np.log(mylist["RR"])- 1.96*(mylist["stdRR"]) mylist["lowRR"] = np.exp(mylist["lowRR"]) mylist["uppRR"] = np.log(mylist["RR"])+ 1.96*(mylist["stdRR"]) mylist["uppRR"] = np.exp(mylist["uppRR"]) mylist["logLR"] = flogLRv(mylist["ntrt"], mylist["nidot"], mylist["ndotj"][0], mylist["ndotdot"][0]) mylist = mylist.sort_values("logLR", ascending=False) cnstR = number np.random.seed(123) mylist["Pvector"] = mylist["nidot"] / mylist["ndotdot"] mylist["Pobs"] = mylist["ntrt"] / mylist["ndotj"] sim_cntdata = np.random.multinomial(mylist["ndotj"][0], mylist["Pvector"], size = cnstR) logLRsim = np.array([flogLRv(i, mylist["nidot"], mylist["ndotj"][0], mylist["ndotdot"][0]) for i in sim_cntdata]) maxlogLRsim = np.array([np.max(i) for i in logLRsim]) mylist["T_alpha"] = np.quantile(maxlogLRsim, 1 - alpha) mylist["pvalue"] = np.array([ np.mean(maxlogLRsim > i) for i in mylist["logLR"]]) if alphaSpending >= 1: mylist["is_signal"] = (mylist["logLR"] > mylist["T_alpha"]) else: mylist["alpha"] = alphaSpending mylist["is_signal"] = (mylist["pvalue"] < alphaSpending) mylist = mylist.sort_values("logLR", ascending=False) return mylist fLRT`34:Td41,heatmap = { var chartWidth = width - margin.left - margin.right, chartHeight = height - margin.top - margin.bottom; var svg = d3.create('svg') .attr("width", chartWidth + margin.left + margin.right) .attr("height", chartHeight + margin.top + margin.bottom); var vis = svg.append("g") .attr("transform", "translate(" + margin.left + "," + margin.top + ")"); var toolTip = d3.select("body").append("div").attr("class", "toolTip"); // Labels of row and columns var xAxisLabel = [...new Set(JSON.parse(LRTResultYear).map(n => n.year))]; var yAxisLabel = [...new Set(JSON.parse(LRTResultYear).map(n => n.AE))].sort(); // Build X scales and axis: var x = d3.scaleBand() .range([ 0, chartWidth ]) .domain(xAxisLabel) .padding(0.01); vis.append("g") .attr("transform", "translate(0," + chartHeight + ")") .call(d3.axisBottom(x)) // Build X scales and axis: var y = d3.scaleBand() .range([ chartHeight, 0 ]) .domain(yAxisLabel) .padding(0.01); vis.append("g") .call(d3.axisLeft(y)); // Build color scale var myColor = d3.scaleLinear() .range(["#fc836a", "#3BAFDA", "lightblue"]) .domain([0, 0.06, 1]); // .range(["lightblue", "#FC6E51"]) // .domain(d3.extent(JSON.parse(LRTResultYear), d => Math.floor(d.logLR/1000)*1000 )); var rect = vis.selectAll() .data(JSON.parse(LRTResultYear), function(d) {return d.year+':'+d.AE;}) .enter() .append("rect"); rect.attr("x", function(d) { return x(d.year) }) .attr("y", function(d) { return y(d.AE) }) .attr("width", x.bandwidth() ) .attr("height", y.bandwidth() ) .style("fill", "lightblue") // .style('opacity', 0.8) .on('mousemove', nodeMouseOver) .on('mouseout', nodeMouseOut ) .transition() .duration(600) .style("fill", function(d) { return myColor(d.pvalue)} ); vis.selectAll(".signalLabel") .data(JSON.parse(LRTResultYear), function(d) {return d.year+':'+d.AE;}) .enter().append("text") .attr("class", "signalLabel") .attr("x", function(d) { return x(d.year) + x.bandwidth()/2 ; }) .attr("y", function(d) { return y(d.AE) + y.bandwidth()/2 + 5 ; }) .text(d => d.is_signal == true ? '❕' : '' ) .style('font-size', '0.65rem') .style('opacity', 0.6); function nodeMouseOver(event, d){ toolTip.style("left", event.pageX + (18) + "px") .style("top", event.pageY + (18) + "px") .style("display", "block") .html(`${d.AE} Year: ${d.year} Is Signal: ${d.is_signal} logLR: ${parseFloat(d.logLR).toFixed(2)} P-value: ${d.pvalue} alpha: ${d.alpha}`); // Optional cursor change on target d3.select(event.target).style("cursor", "pointer"); // Optional highlight effects on target d3.select(event.target) .transition() .attr('stroke', '#aeaeae') .attr('stroke-width', 2); } function nodeMouseOut(event, d){ toolTip.style("display", "none"); d3.select(event.target).style("cursor", "default"); d3.select(event.target) .transition() .attr('stroke', '#fff') .attr('stroke-width', 0); } return svg.node() }35:T474,html ` pillCSS `36:T408,// https://getbootstrap.com/docs/4.1/getting-started/introduction/#css html`css `37:[] 2f:{"title":"LRT Method Application: Part 2","cells":[{"id":0,"value":"# **LRT Method Application: Part 2**","pinned":false,"hidden":false,"mode":"md"},{"id":930,"value":"pills = {\n const pillData = ['Pyodide', 'Python in JavaScript', 'Monte Carlo simulation', 'd3-heatmap'];\n return html `\n

\n ${pillData.map(d => '

'+d+'

')}\n

\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":20,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":382,"value":"viewof numberOfSimulation = slider({\n min: 100, \n max: 10000, \n step: 100,\n value: 100, \n title: \"Number of Simulation\", \n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":767,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":830,"value":"### **Overall LRT Analysis** ","pinned":false,"hidden":false,"mode":"md"},{"id":52,"value":"**LRT function**
\nThis function has been rewritten in Python since we can use the numpy.random.multinomial to easily conduct the Monte Carlo simulation.
\nPyodide library is used to run the Python codes in Observable.","pinned":false,"hidden":false,"mode":"md"},{"id":166,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":826,"value":"Run this function with two parameters:\n* Name of raw data\n* Number of simulation","pinned":false,"hidden":false,"mode":"md"},{"id":319,"value":"LRTResult = py`import pandas as pd\n${funcLRT}(pd.read_json(${JSON.stringify(LRTData)}),${numberOfSimulation}).to_json(orient=\"records\")`","pinned":true,"hidden":false,"mode":"ojs"},{"id":908,"value":"LRTData","pinned":true,"hidden":false,"mode":"ojs"},{"id":820,"value":"In the table below, \"logLR\" is the statistic from the observed data and \"T_alpha\" is the threshold value.","pinned":false,"hidden":false,"mode":"md"},{"id":919,"value":"Inputs.table(JSON.parse(LRTResult))","pinned":false,"hidden":false,"mode":"ojs"},{"id":396,"value":"LRTYearData","pinned":true,"hidden":false,"mode":"ojs"},{"id":834,"value":"### **LRT Trend Analysis by Year**","pinned":false,"hidden":false,"mode":"md"},{"id":843,"value":"**Iterate the LRT function by year** \n\nThe LRT method was also applied to longitudinal data from ${d3.min(yearList)} to ${d3.max(yearList)} where a monotone **alpha-spending** function was used to ensure the control of type I error at each look and overall for trend analysis.\n\n ${tex`\\large \\alpha_{(k)} = \\frac{0.05}{2^k}, \\small k=0, 1, 2, \\cdots `}
\n\n\nThat is, alpha over 2 to power k. Here k=0 (means 1st time period), 1 (2nd time peroid), and so on.","pinned":false,"hidden":false,"mode":"md"},{"id":387,"value":"LRTResultYear = py`import pandas as pd\nimport numpy as np\nLRT_result_Year = pd.DataFrame()\n\nfor k in range(len(${yearList})):\n LRTYearData = pd.read_json(${JSON.stringify(LRTYearData)})\n countAeListYeark = LRTYearData[LRTYearData[\"Year\"] <= ${yearList}[k]].groupby(\"AE\").agg({\"ntrt\": \"sum\"})\n countAeListTotalYeark = LRTYearData[LRTYearData[\"Year\"] <= ${yearList}[k]].groupby(\"AE\").agg({\"nidot\": \"sum\"})\n combine_data_Yeark = pd.concat([countAeListYeark, countAeListTotalYeark], axis=1, join=\"inner\")\n combine_data_Yeark[\"ndotj\"] = np.sum(combine_data_Yeark[\"ntrt\"])\n combine_data_Yeark[\"ndotdot\"] = np.sum(combine_data_Yeark[\"nidot\"])\n alphaSpending = 0.05 / (2**k)\n LRT_result_Yeark = ${funcLRT}(combine_data_Yeark, ${numberOfSimulation}, alphaSpending = alphaSpending)\n LRT_result_Yeark[\"year\"] = str(${yearList}[k])\n LRT_result_Year = LRT_result_Year.append(LRT_result_Yeark)\nLRT_result_Year['AE'] = LRT_result_Year.index\n\nLRT_result_Year.to_json(orient=\"records\")`","pinned":true,"hidden":false,"mode":"ojs"},{"id":753,"value":"aq.from(JSON.parse(LRTResultYear)).view({limit:JSON.parse(LRTResultYear).length})","pinned":false,"hidden":false,"mode":"ojs"},{"id":697,"value":"{\n return d3.selectAll('table tr').each(function() { \n var container = [];\n d3.select(this).selectAll('td').each(function() {\n container.push(d3.select(this).html());\n });\n console.log(container);\n container.indexOf('true') > -1 ? d3.select(this).selectAll('td').style('color', 'blue') :d3.select(this).selectAll('td').style('color', 'black')\n })\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":861,"value":"**Trend Analysis Visualization**\n\nThe heatmap consists of all adverse events and different years. The color code indicates the p-value of one adverse reaction in a certain year. The darker the red, the smaller the p-value; The lighter the blue, the greater the p-value. An exclamation mark (❕) represents that the adverse event is identified as a signal in a certain year (hover over on each cell to see the details). \n","pinned":false,"hidden":false,"mode":"md"},{"id":434,"value":"$34","pinned":true,"hidden":false,"mode":"ojs"},{"id":877,"value":"

\n\n### **Reference**
\n\n

Lan Huang, Jyoti Zalkikar, and Ram Tiwari. Likelihood ratio test based method for signal detection in drug classes using FDA’s AERS database. Journal of Biopharmaceutical Statistics, Special Issue on Drug Safety, Vol23(1), 2013

Jianjin Xu, Lan Huang, Zhihao Yao, et al. Statistical methods for clinical study site selection for inspection, Therapeutic Innovation & Regulatory Science (TIRS), 2018

Zhiheng Xu, Jianjin Xu, Zhihao Yao, Lan Huang, Mary Jung, Ram Tiwari (2020) Evaluating Medical Device Safety Signals Using a Likelihood Ratio Test Method, JBS, 28;1-10, DOI: 10.1080/10543406.2020.1783284

","pinned":false,"hidden":false,"mode":"md"},{"id":5,"value":"

\n### **Appendix**","pinned":false,"hidden":false,"mode":"md"},{"id":479,"value":"height = 1900","pinned":false,"hidden":false,"mode":"ojs"},{"id":481,"value":"width = 800 ","pinned":false,"hidden":false,"mode":"ojs"},{"id":914,"value":"import {drugName, dateFrom, dateTo, maxNumberOfEvent, LRTData, yearList, LRTYearData} from '@howardyao/lrt-method-application-part-1'","pinned":false,"hidden":false,"mode":"ojs"},{"id":485,"value":"margin = ({top: 30, right: 30, bottom: 50, left: 240})","pinned":false,"hidden":false,"mode":"ojs"},{"id":378,"value":"import {date, slider, text} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":36,"value":"import {form} from \"@mbostock/form-input\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":15,"value":"import { aq, op } from '@uwdata/arquero'","pinned":false,"hidden":false,"mode":"ojs"},{"id":163,"value":"import { py, pyodide } from '@gnestor/pyodide'","pinned":false,"hidden":false,"mode":"ojs"},{"id":927,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":12,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"d3 = require(\"d3@6\")","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W37","stdlib":"1","titling":"auto"}

LRT Application: Signal Detection from FDA Adverse Event Reporting System (FAERS)

Part 2: Signal Detection and Visualization

LRT Application:
Signal Detection from FDA Adverse Event Reporting System (FAERS)