32:T9f9,# Visualizing Local Moran's I using dual-density plots The Local Moran's I statistic is a useful tool for detecting patterns in spatial data. It is used to detect both spatial clusters and spatial outliers. The Local Moran's I statistic requires a dataset comprised of spatial locations (e.g. US counties) each with an associated value for a particular measure (e.g. cancer mortality rate). The statistic is a property of a single location, so each location in the dataset will have its own Local Moran's I value. The dataset's values must be normalized so that the mean is equal to 0 (i.e. values above the mean are positive, and values below the mean are negative). Local Moran's I is based on the concept of **spatial lag**: a weighted mean of the values of a location's spatial neighbors. For the sake of simplicity, we will be assuming equal weights in this notebook, meaning each location's spatial lag will be equal to the mean of its location's neighboring values. Taken on its own, the Local Moran's I statistic is capable of detecting two basic spatial patterns: * **Positive spatial autocorrelation:** If the location’s value and spatial lag both deviate from the dataset’s mean in the same direction – either both higher or both lower. This is a potential spatial cluster, i.e. either a hot spot or cold spot. * **Negative spatial autocorrelation:** If the location’s value and spatial lag deviate from the mean in opposite directions – one higher and one lower. This is a potential spatial outlier. A result is deemed significant only if a location’s spatial lag deviates from the mean more than would be expected by random chance. For all other results, we say that the location exhibits no spatial autocorrelation. To learn more about Local Moran's I (and other similar statistics) we recommend [this excellent resource](https://geodacenter.github.io/workbook/6a_local_auto/lab6a.html) from the GeoDa group. ## Plotting Local Moran's I This notebook is primarily to discuss a new type of plot we have designed for visualizing Local Moran's I. We call these plots "Moran dual-density plots" because they consist of two density plots joined together, and marked with other information relevant to a single Local Moran's I calculation. Basically, Moran dual-density plots show the main elements that go into a Local Moran's I calculation and how these elements are related to each other. To explain how these plots work, we will first introduce an expanded and labelled version, before introducing a denser version.33:T67c,{ const featureMap = d3.index(counties.features, d => d.id) const resultMap = d3.index(results, d => d.id) const focusFeature = featureMap.get(result.id) const neighborMap = d3.group(weights, d => d[0]) const neighborWeights = neighborMap.get(focusFeature.id) const neighborFeatures = neighborWeights.map(d => featureMap.get(d[1])) const neighborsCombinedFeature = turf.combine(turf.featureCollection(neighborFeatures)) const dist = estimateDistribution(results.map(d => d.z)) const lagDist = estimateDistribution(results.map(d => d.lag)) return Plot.plot({ width: 300, color: { scheme: "RdYlBu", legend: true, pivot: d3.mean(results, d => d.z), reverse: true, label: "Mortality Rate (Z-normalized)" }, projection: {type: "albers-usa", domain: neighborsCombinedFeature }, marks: [ Plot.geo(counties.features.filter(d => d.id != focusFeature.id), { stroke: "lightgrey", fill: d => resultMap.get(d.id)?.z, opacity: .2 }), Plot.geo(neighborFeatures , { stroke: "lightgrey", strokeWidth: 1.5, fill: d => resultMap.get(d.id)?.z, }), Plot.geo(counties.features.filter(d => d.id == focusFeature.id) , { stroke: "yellow", strokeWidth: 3, fill: d => resultMap.get(d.id)?.z, }), // Plot.dot(neighborCentroids, { x: d => d[0], y: d => d[1] }), // Plot.text([focusFeature, ...neighborFeatures], // Plot.centroid({ // text: (_,i) => i, // fontSize: 15, // textAnchor: "middle", // stroke: "white", // fill: "black" // }) // ) ] }) }34:T174f,{ const div = document.createElement("div") div.style.display = "flex" div.style.alignItems = "middle" const featureMap = d3.index(counties.features, d => d.id) const resultMap = d3.index(results, d => d.id) const focusFeature = featureMap.get(result.id) const neighborMap = d3.group(weights, d => d[0]) const neighborWeights = neighborMap.get(focusFeature.id) const neighborFeatures = neighborWeights.map(d => featureMap.get(d[1])) const neighborsCombinedFeature = turf.combine(turf.featureCollection(neighborFeatures)) const dist = estimateDistribution(results.map(d => d.z)) const lagDist = estimateDistribution(results.map(d => d.lag)) const mapPlot = Plot.plot({ width: 500, color: { scheme: "RdYlBu", legend: true, pivot: d3.mean(results, d => d.z), reverse: true, label: "Mortality Rate (Z-normalized)" }, projection: {type: "albers-usa", domain: neighborsCombinedFeature }, marks: [ Plot.geo(counties.features.filter(d => d.id != focusFeature.id), { stroke: "lightgrey", fill: d => resultMap.get(d.id)?.z, opacity: .2 }), Plot.geo(neighborFeatures , { stroke: "lightgrey", strokeWidth: 1.5, fill: d => resultMap.get(d.id)?.z, }), Plot.geo(counties.features.filter(d => d.id == focusFeature.id) , { stroke: "yellow", strokeWidth: 3, fill: d => resultMap.get(d.id)?.z, }), // Plot.dot(neighborCentroids, { x: d => d[0], y: d => d[1] }), Plot.text([focusFeature, ...neighborFeatures], Plot.centroid({ text: (_,i) => i, fontSize: 15, textAnchor: "middle", stroke: "white", fill: "black" }) ) ] }) const resultIndex = d3.index(results, d => d.id) const neighborResults = neighborWeights.map(d => resultIndex.get(d[1])) const zExtent = d3.extent(results, d => d.z) const zDistExtent = d3.extent(dist, d => d[0]) const moranExtent = zDistExtent.map(d => d * result.z) const scaleMap = d3.scaleLinear( moranExtent, zDistExtent) const explainPlot = Plot.plot({ grid: true, height: 300, marginTop: 40, marginBottom: 40, marginLeft: 100, x: { axis: "top", domain: zDistExtent, label: "Z-normalized value" }, style: { fontSize: 14, color: "grey" }, y: { domain: [ "Focal value", "Neighbors", "Lag", "Statistic" ] }, marks: [ Plot.ruleX([0], {stroke: "#ebebeb", strokeWidth: 2}), Plot.dot(neighborResults, {x: "z", y: d => "Neighbors", fill: "black"}), Plot.text(neighborResults, {text: (_,i) => i+1, x: "z", y: d => "Neighbors", dy: -12, fontSize: 14, fill: "black"}), Plot.link(neighborResults, { x1: "z", y1: d => "Neighbors", x2: d => result.lag, y2: d => "Lag", stroke: "black", strokeDasharray: "1,3", strokeOpacity: .4, curve: "bump-y" }), Plot.link([result], { x1: "lag", y1: d => "Lag", x2: "lag", y2: d => "Statistic", stroke: "blue", strokeDasharray: "1,3", strokeOpacity: 1, strokeWidth: 2 }), Plot.dot([result], {x: "z", y: d => "Focal value", fill: "black", r: 3}), Plot.arrow([result], {x2: "z", x1: d => 0, y: d => "Focal value", stroke: "black", insetEnd: 5, strokeDasharray: "2,2", strokeWidth: 1.2, headLength: 8 }), Plot.text([result], {text: (_,i) => i, x: "z", y: d => "Focal value", dy: -14, fontSize: 14, fill: "black"}), Plot.dot([result], { x: "lag", fill: "blue", y: d => "Lag",}), Plot.text([result], { text: d => "Spatial lag", x: "lag", fill: "blue", y: d => "Lag", dx: 10, textAnchor: "start"}), Plot.axisX(scaleMap.ticks(), { anchor: "bottom", label: "Local Moran's I", x: scaleMap, tickFormat: scaleMap.tickFormat(), labelArrow: result.z > 0 ? "right" : "left"}), Plot.dot([result], Plot.mapX((D) => D.map(scaleMap), { x: d => d.statistic, y: d => "Statistic", fill: "blue", dy: 3, })), Plot.text([result], Plot.mapX((D) => D.map(scaleMap), { text:d => "Local Moran's I", x: d => d.statistic, y: d => "Statistic", fill: "blue", dx: 10, textAnchor: "start" })), // Plot.frame(), ] }) const explainDiv = document.createElement("div") explainDiv.style.display = "flex" explainDiv.style.flexDirection = "column" const zDistPlot = Plot.plot({ height: 60, marginLeft: 90, x: { domain: zDistExtent, axis:null }, y: { axis: null }, marks: [ Plot.areaY(dist, {x: d=> d[0], y: d => d[1], curve: "basis", fill: "lightgrey"}), Plot.lineY(lagDist, {x: d=> d[0], y: d => d[1], curve: "basis", stroke: "blue", strokeDasharray: "2,2", opacity: .3}), // Plot.ruleX([0]) ] }) const statDistPlot = Plot.plot({ height: 60, marginLeft: 80, x: { domain: zDistExtent.map(d => d * result.z), axis:null }, y: { axis: null, reverse: true }, marks: [ Plot.areaY(result.permutationDistribution, {x: d=> d[0], y: d => d[1], curve: "basis", fill: "lightgrey"}), Plot.ruleX([result.lowerCutoff], { stroke: "purple", strokeDasharray: "2,2" }), Plot.ruleX([result.upperCutoff], { stroke: "green", strokeDasharray: "2,2" }), Plot.text([result], {text: d => "Cutoff (Negative)", x: "lowerCutoff", dx: result.z < 0 ? 5 : -5, textAnchor: result.z < 0 ? "start" : "end", fontSize: 13, fill: "purple"}), Plot.text([result], {text: d => "Cutoff (Positive)", x: "upperCutoff", dx: result.z >= 0 ? 5 : -5, textAnchor: result.z >= 0 ? "start" : "end", fontSize: 13, fill: "green"}), Plot.dot([result], { x: "statistic", fill: "blue", y: 0, dy: 3}) // Plot.ruleX([0]) ] }) explainDiv.appendChild(zDistPlot) explainDiv.appendChild(explainPlot) explainDiv.appendChild(statDistPlot) div.appendChild(mapPlot) div.appendChild(explainDiv) return div }35:T72a,estimateDistribution = { const gaussian = (await import('https://cdn.jsdelivr.net/npm/gaussian@1.3.0/+esm')).default function kde(X, K, h) { return function (x) { let sum = 0 for (let xi of X) { sum += K((x - xi) / h) } return (1/(X.length*h)) * sum } } function estimateDistribution(X, n=50, extent=null) { const normal = gaussian(0,1) const h = 0.9 * d3.deviation(X) * X.length ** (-1/5) const kernel = kde(X, d => normal.pdf(d), h) if (!extent) { extent = d3.extent(X) } const threshold = 0.001 const midpoint = (extent[1] + extent[0])/2 let points = [] let step = (extent[1] - extent[0])/n for (let i = 0; i < n; i++) { const x = midpoint - step*i const value = kernel(x) points.push([x,value]) if (value < threshold) { break } } points = points.reverse() for (let i = 1; i < n*2; i++) { const x = midpoint + step*i const value = kernel(x) points.push([x,value]) if (value < threshold) { break } } const max = d3.max(points, d => d[1]) return points.map(([x,y]) => [x, y]) } function interpolatePoints(points, xVals) { const resultPoints = [...points] for (const x of xVals) { let i = 1 while (i < points.length && points[i][0] < x) { i++ } if (i === 0 || i === points.length) { resultPoints.push([x, NaN]) continue } const x0 = points[i - 1][0] const y0 = points[i - 1][1] const x1 = points[i][0] const y1 = points[i][1] const y = y0 + (x - x0) * (y1 - y0) / (x1 - x0) resultPoints.push([x, y]) } return resultPoints.sort((a,b) => a[0] - b[0]) } return estimateDistribution }36:[{"id":"ced0e14962d45d008602d67e70ca82e30a14849927bb93cb4b06dc61450cd7c69f604036cc56ff3a6bf56783807c4bc8ce5a42eea3ea148f46e0f2cdb69a449d","href":"https://static.observableusercontent.com/files/ced0e14962d45d008602d67e70ca82e30a14849927bb93cb4b06dc61450cd7c69f604036cc56ff3a6bf56783807c4bc8ce5a42eea3ea148f46e0f2cdb69a449d","name":"stateFips.csv","type":"text/csv","lastModified":1720213150386.859,"size":814}] 30:{"title":"Visualizing Local Moran's I using dual-density plots","cells":[{"id":0,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":388,"value":"viewof result = {\n const featureIndex = d3.index(counties.features, d => d.id)\n return Inputs.select(results, {label: \"Select county\", format: d => `${featureIndex.get(d.id)?.properties.name}, ${stateFipsMap.get(d.id.slice(0,2)).short} (${resultTextShort(d)})`, value: results.find(d => d.id == \"56011\")})\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":770,"value":"Plot.legend({color: { scheme: \"RdYlBu\", pivot: d3.mean(results, d => d.z), reverse: true, label: \"Mortality Rate (Z-normalized)\" }})\n","pinned":true,"hidden":false,"mode":"ojs"},{"id":778,"value":"Plot.legend({color: {type: \"linear\", scheme: \"RdYlBu\", pivot: d3.mean(results, d => d.z), reverse: true, label: \"Mortality Rate (Z-normalized)\"}})","pinned":true,"hidden":false,"mode":"ojs"},{"id":761,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":364,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":584,"value":"This location exhibits ${resultText(result)}. \n\nTry selecting different locations using the selection input above and seeing how they are visualized (the example data here is US cancer mortality from 2010-2020).","pinned":false,"hidden":false,"mode":"md"},{"id":628,"value":"A Moran dual-density plot shows the following:\n\n* **Density plot of values:** An estimated density plot of all values in the dataset. In Moran's I analysis, these are typically z-score normalized as we have done here (for each value: subtract the mean, divide by the standard deviation).\n* **Density plot of spatial lags:** An estimated density plot of all spatial lags in the dataset (optional).\n* **The focal location value and the neighbor values.** Again, z-score normalized.\n* **The spatial lag.** Dashed grey lines explicitly show that this value is calculated from the neighbor values.\n* **The Local Moran's I statistic value**. A vertical blue line conveys how this is directly proportional to the spatial lag.\n* **Density plot of permuted statistic values:** To determine significance, the statistic's value is typically compared to an empirical distribution generated using bootstrapping. We plot this distribution in a density plot, alongside the cut-offs for significance.\n\n","pinned":false,"hidden":false,"mode":"md"},{"id":687,"value":"Something worth mentioning for anybody unfamiliar to with the Local Moran's I statistic: the magnitude of the focal location's z-score normalized value has no effect on the pseudo p-value calculation. However, it's SIGN determines the handedness of the distribution (if the value is negative then the distribution will be mirrored compared to if the value is positive).\n\nWe hope the Moran dual-density plot makes the core ideas of Local Moran's I clear. ","pinned":false,"hidden":false,"mode":"md"},{"id":700,"value":"### The final plot","pinned":false,"hidden":false,"mode":"md"},{"id":660,"value":"And now we condense and rearrange these elements a little to give the final proposed plot:","pinned":false,"hidden":false,"mode":"md"},{"id":724,"value":"viewof displayLag = Inputs.select([\"none\", \"stroke\", \"fill\"], {label: \"Display Lag Distribution\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":750,"value":"viewof areaMode = Inputs.select([\"hot_cold\", \"positive_negative\", \"cluster_labels\"], {label: \"Area Labels\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":733,"value":"viewof labelAreas = Inputs.toggle({label: \"Label Areas\", value: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":438,"value":"{\n const dualDensityPlot = new moranplot.MoranDualDensity(results, {\n width: 640, height: 200, \n fontSize: 14,\n lagDistribution: displayLag != \"none\" ? \"auto\" : null,\n lagDistributionDisplay: displayLag,\n labelAreas,\n areaMode,\n colors: {\n positiveAutocorrelation: \"#8fe39a\",\n negativeAutocorrelation: \"#9b81a6\",\n }\n })\n dualDensityPlot.focus(result.id)\n return dualDensityPlot.plot()\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":717,"value":"To use these plots in your own work, check out the [MoranPlot library](https://observablehq.com/@siliconjazz/moranplot-overview).\n\nSee the plot working into a [spatial dashboard here](https://episphere.github.io/moranplot/). For more information on the other plots in that dashboard, see the the [MoranPlot notebook](https://observablehq.com/@siliconjazz/moranplot-overview).","pinned":false,"hidden":false,"mode":"md"},{"id":569,"value":"## Working","pinned":false,"hidden":false,"mode":"md"},{"id":9,"value":"import { counties } from \"@observablehq/plot-mapping\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":92,"value":"import { results as oldResults } from \"@siliconjazz/moranplot-overview\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":529,"value":"results = {\n const results = spatialAutocorrelation.spatialAutocorrelation(oldResults, \"value\", { \n method: \"local_moran_i\", weightTuples: weights, permutationOutput: { granularity: 20 }, outputCutoffs: true, outputNeighborWeights: true\n })\n const neighborMap = d3.group(weights, d => d[0])\n results.forEach(result => {\n result.neighbors = result.neighborWeights\n })\n // results.forEach(result => {\n // const neighbors = neighborMap.get(result.id)\n // if (neighbors) {\n // result.neighbors = neighbors.map(d => d.slice(1))\n // }\n // })\n return results\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":396,"value":"stateFipsMap ={\n const data= await FileAttachment(\"stateFips.csv\").csv()\n return d3.index(data, d => d.fips)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":32,"value":"weights = spatialAutocorrelation.calculateEqualContiguousWeightTuples(counties)","pinned":false,"hidden":false,"mode":"ojs"},{"id":432,"value":"## Helper ","pinned":false,"hidden":false,"mode":"md"},{"id":434,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":590,"value":"function resultText(result) {\n if (result.statistic > 0 && result.p <= 0.05) {\n return htl.html`significant positive local spatial autocorrelation`\n } else if (result.statistic < 0 && result.p <= 0.05) {\n return htl.html`significant negative local spatial autocorrelation`\n } else {\n return htl.html`no significant local spatial autocorrelation`\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":602,"value":"function resultTextShort(result) {\n if (result.statistic > 0 && result.p <= 0.05) {\n return \"Positive SA, \" + result.label\n } else if (result.statistic < 0 && result.p <= 0.05) {\n return \"Negative SA\"\n } else {\n return \"No SA\"\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":26,"value":"## Imports","pinned":false,"hidden":false,"mode":"md"},{"id":465,"value":"moranplot = import(\"https://episphere.github.io/moranplot/src/export.js\")\n //import(\"http://127.0.0.1:5500/src/export.js?d=\"+new Date())","pinned":false,"hidden":false,"mode":"ojs"},{"id":45,"value":"turf = require('@turf/turf@7.0.0/turf.min.js')","pinned":false,"hidden":false,"mode":"ojs"},{"id":29,"value":"spatialAutocorrelation = import(\"https://episphere.github.io/spatial-autocorrelation/src/autocorrelation.js\")\n //import(\"http://127.0.0.1:5500/src/autocorrelation.js\")","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W36","stdlib":"1","titling":"auto"}