40:null 41:["$","title",null,{"children":"Loess with optimal selection of span | Simon Jackman | Observable"}] 47:T6cd, ## Parameterizing smoothness - form a neighbourhood around a target point ${tex`x_0`}, ${tex`N(x_0)`} - the size of ${tex`N(x_0)`} is parameterised by ${tex`\lambda \in (0,1]`} the proportion of the data set encompassed by ${tex`N(x_0)`} and usually referred to as the **span** of the smoother. - inside that neighbourhood, fit a low-order polynomial regression (usually linear or quadratic), the regression fit to minimise a weighted sum of residuals criterion, weighting observation ${tex`i \in N(x_0)`} by a function of the distance of ${tex`x_i`} from the target point ${tex`x_0`}. - the functional form of the weighting function is the **kernel** of the smoother; for _loess_, a tricube weighting function is the standard default, i.e., ${tex`w(d) = (1 - |d|^3)^3`} where ${tex`d`} is the distance between ${tex`i \in N(x_0)`} and ${tex`x_0`}, normalised to range from 1 at the edge of ${tex`N(x_0)`} to 0 for ${tex`x_i = x_0`}. - the predicted value of ${tex`y`} at ${tex`x_0`} from the locally weighted regression supplies ${tex`\hat{y} | x_0 = \hat{f_\lambda}(x_0)`}. - given a sufficiently fine grid of target points, the locus of resulting output points is plotted as a line (or a series of connected line segments), a _scatterplot smooth_ of ${tex`y`} against ${tex`x`}. - As ${tex`\lambda \in (0,1] \rightarrow 1`}, _loess_ tends towards linear regression; as ${tex`\lambda \rightarrow 0`}, _loess_ merely linearly interpolates each unique ${tex`(x_i, y_i)`} point ordered by ${tex`x`}. Smoothness is therefore parameterised by ${tex`\lambda`}. Most implementations here on observablehq.com include an `Input` slider to permit the user to select or experiment with different values of ${tex`\lambda \in (0,1]`}. 48:T419,## Example On 14 October 2023, Australia held a referendum on whether to amend its Constitution to "recognise the First Peoples of Australia by establishing an Aboriginal and Torres Straight Islander Voice", an advisory body that "may make representations to the Parliament and the Executive Government of the Commonwealth on matters relating to Aboriginal and Torres Strait Islander peoples". Vote shares (Yes as a proportion of valid votes) are available at the polling place level, along with various polling place level covariates, including vote shares for political parties cast at the 2022 federal election and estimates of demographic characteristics of the group of voters utilising each polling place. We use _loess_ to summarise the relationship between ${tex`y`} (Yes %) and a set of covariates ${tex`x`} (accessible with the "Correlate" selection), at the the polling place level, by Australian state or territory. As is conventional an `Input` slider is presented to permit users to experiment with the _loess_ `span` parameters. 49:T42b,Plot.plot({ //grid: true, y: { label: "Yes %", percent: true, domain: [0, 100] }, x: { label: form.correlate.value, percent: form.correlate.name !== "median_age", }, r: { range: [.5, 10] }, marks: [ //Plot.frame(), Plot.dot( pp_res.filter(d => d.state === form.theState), { x: form.correlate.name, y: "yes_per", r: "formal", strokeWidth: 0.5, alpha: 0.2 } ), Plot.line( myLoess2( pp_res .filter(d => d.state === form.theState & !isNaN(d.yes_per) &!isNaN(d[form.correlate.name])) .slice() .sort(function(x,y){return d3.ascending(x[form.correlate.name],y[form.correlate.name]);}), { x: form.correlate.name, y: "yes_per", span: form.span, degree: form.degree } ), { x: "x", y: "y", strokeWidth: 1.5, stroke: "blue", } ) ] } ) 4a:T4ec,## _loess_ is slow in Javascript The _loess_ implementation I use here is a [`JS`](https://github.com/yongjun21/loess) port of Cleveland, Grosse and Shyu's (1992) implementation in `C`. The function I define below, `myLoess2` takes a data set and some options as input, calls `Loess` and returns an object suitable for passing to `Plot.line`. For the larger states with polling places numbering in the thousands (NSW, VIC and QLD), computation of a single _loess_ fit _to the entire data set_ is _slow_, taking well over two seconds on my hardware and browser (MacBook Pro M1 Max, latest Google Chrome) and updates to `span` are laggy. For instance, the NSW section of the data set is reasonably large, with over 3,000 observations, encompassing all polling places in the state of NSW. Producing predicted values for each data point takes 8 to 9 seconds on my hardware; repeated fits to the entire data set, exploring different values of `span`, is simply infesible. We attempt to overcome this by computing predictions at a much smaller number of _evaluation points_ (defaulting to 50 points). Searching over different `span` parameters is also done by evaluating each loess fit at this smaller number of evaluations points too, as elaborated below.4b:T51c,## Optimal smoothness via cross-validation How to choose the _loess_ `span` parameter? As ${tex`\lambda \rightarrow 0`} the fitted values produced by _loess_ ${tex`\hat{y}_i | x_i = \hat{f_\lambda}(x_i) \rightarrow y_i`} and the _loess_ residuals ${tex`e_i = y_i - \hat{y}_i \rightarrow 0`}. Setting ${tex`\lambda`} too small therefore _overfits_ the data, while setting ${tex`\lambda`} too large _oversmooths_ the data (the fit tending towards a simple linear regression), masking local features of the data. A popular approach to selecting smoothing parameters is to pose the problem as one of minimising out-of-sample predictive error. This idea encapsulates the idea that bandwidth selection is akin to negotiating a bias/variance tradeoff with respect to the fit of the model to the sample data. With ${tex`\lambda`} small, _loess_ generates low bias predictions due to the _local_ nature of the fitting, but with higher variance for the same reasons. Widening the bandwidth puts more data into the local neighbourhoods underlying the smoother, driving down variance, but at the cost of more bias as the fitting becomes less "local". Optimally negotiating this bias-variance tradeoff corresponds to choosing ${tex`\lambda`} so as to minimize the _cross-validation sum of squared residuals_4c:T5b7,where the subscript ${tex`(-i)`} denotes "leave out observation ${tex`i`} and ${tex`h_i`} is the ${tex`i`}-th element of the diagonal of ${tex`H(x)`}. Note that remarkable fact the "leave-out-${tex`i`}" residuals can be computed from the full-sample residuals ${tex`e_i`} divided by ${tex`1 - h_i`}. Further, among the interesting and useful properties of ${tex`H(x)`} is that its _trace_ ${tex`\sum_i h_i = k`}, the number of predictors in the regression (including any included intercept term). This means that ${tex`\bar{h} = \sum_i h_i/n = k/n`} and suggests that the terms in the denominator ${tex`1 - h_i`} be approximated by ${tex`1 - \bar{h} = 1-k/n`}, obviating direct computation of the hat matrix. Note too that ${tex`k=\text{trace}(H)`} is the degrees of freedom consumed by the regression. _Generalized cross-validation_ (GCV) results from replacing ${tex`1 - h_i`} with the approximation ${tex`1 - \text{trace}(H)/n = 1 - k/n`}, involving no additional computation beyond that ordinarily required for linear regression. ## Application to _loess_ For scatterplot smoothing, we don't wind up with a _hat_ matrix but rather a _smoothing_ matrix, ${tex`S_\lambda`} that - is a function of ${tex`x`} and the smoothing parameter ${tex`\lambda`} - ${tex`\text{trace}(S)`} is a good approximation for the degrees of freedom of the smoother. ### GCV for _loess_ The GCV sum of squared errors for _loess_ with span ${tex`\lambda`} is therefore 4d:T97a,But how to compute the trace of the smoother matrix ${tex`S_\lambda`}? For _loess_, a computational trick for obtaining the ${tex`i`}-th element of the diagonal of the smoother matrix given a candidate value of ${tex`\lambda`} is as follows: - sort ${tex`n`} from low to high - for ${tex`i`} = 1, ... ${tex`n`} (the length of ${tex`x`}) - construct a vector ${tex`z`} of length ${tex`n`}, with ${tex`z_i = 1`} and 0 otherwise. - with the candidate `span` parameter ${tex`\lambda`}, use _loess_ to predict ${tex`z`} as a function of ${tex`x`} - ${tex`\hat{z}_i = \hat{f}_{\lambda}(x_i)`} is the ${tex`i`}-th element of the diagonal of the smoother matrix ${tex`S_\lambda`}. - ${tex`\sum_i \hat{z}_i = \text{trace} S_\lambda`}. Incidentally, note that as is the case with the hat matrix in the linear regression context, the smoother matrix is solely a function of ${tex`x`} and ${tex`\lambda`}. ## An approximation of the trace of the smoother matrix Computing ${tex`\sum_i \hat{z}_i = \text{trace} S_\lambda`} this way requires _n_ runs of _loess_ for each candidate value of the span parameter ${tex`\lambda`}; the timing experiments reported above make clear that we will not be adopting this approach here. A compromise is that for applications involving more than ${tex`n`} = ${tex`M`} data points we perform the steps above at ${tex`M`} points in the domain of ${tex`x`}, evenly spaced across the empirical cumulative distribution function of ${tex`x`}. Experiments show as few as ${tex`M= 50`} evaluations provide reasonable performance in terms of speed and accuracy, over-estimating ${tex`\text{trace} S_\lambda`} by a trivial amount, around 1 or 2%. The approximation gets better with more evaluation points and even as many as 250 evaluations can produce good results in a tolerable amount of clock time. This supplies a reasonable approximation to the degrees of freedom used by the _loess_ smooth. Note that implementations of _loess_ in statistical programs like `R` use compiled `C` and `Fortran` that run extremely quickly relative to the `js` implementation used here. Further, the `Fortran` code underlying the implementation in `R` uses good approximations to ${tex`\text{trace} S_\lambda`} developed by Cleveland and Grosse (1991). Porting these approximations to `JS` implementations of _loess_ would be a useful step, or else having the ability to call `R` server-side for these computations.4e:T880,test_of_span = async () => { const testdata = pp_res .filter(d => d.state=="NSW") .map(d => [d.uni, d.yes_per]) .filter(d => !isNaN(d[0]) & !isNaN(d[1])) .slice() .sort(function(x,y){return d3.ascending(x[0],y[0]);}) const bw = .15; const npoints = [ 10, 30, 50, 80, 100, 250, 500, 1000, testdata.length]; function loess_trace(loess_data,bandwidth,M){ console.log("loess_trace has bandwidth " + bandwidth); const x_vec = loess_data.map(d => d[0]); //console.log("loess_trace x"); //console.log(x_vec); const n = x_vec.length const nindx = n > M ? M : n; var index = new Array(nindx); if(n>M){ for(let i = 0; i < nindx; i++){ index[i] = Math.floor(i/nindx * n); // indices of evaluation points } } else { for(let i = 0; i < nindx; i++){ // every point is an evaluation point index[i] = i; } } var tmpData = []; // copy data for(let i=0; i d[0]), y: tmpData.map(d => d[1]) }, { span: bandwidth, degree: 1 } ); const x0 = { x: [ x_vec[index[i]] ] }; //console.log(x0); const traceContribution = f.predict(x0).fitted[0]; if(isNaN(traceContribution)){ console.log("bad trace contribution at observation " + i); console.log(tmpData[i]); } else { trace += traceContribution; } tmpData[i][1] = 0.0; } return trace; } var startTime, endTime, elapsedTime; var z; var out = []; for(let i=0; i < npoints.length; i++){ startTime = new Date(); z = loess_trace(testdata,bw,npoints[i]); endTime = new Date(); elapsedTime = endTime - startTime; out.push( { elapsedTime: elapsedTime, npoints: npoints[i], z: z } ); } return out; }4f:T7fc,function loess_gcv_local( data, {x,y,...options} ){ const startTime = new Date(); const span = options.span ?? 0.5; const degree = options.degree ?? 1; const M = options.eval ?? 50; const model = new Loess( { x: data.map(d => d[x]), y: data.map(d => d[y]) }, { span: span, degree: degree } ); // set up evaluation grid const x_vec = data.map(d => d[x]); const n = x_vec.length; const nindx = n > M ? M : n; var index = new Array(nindx); if(n>M){ for(let i = 0; i < nindx; i++){ index[i] = Math.floor(i/nindx * n); // indices of evaluation points } } else { for(let i = 0; i < nindx; i++){ // every point is an evaluation point index[i] = i; } } var tmpData = []; // copy data for(let i=0; i d[0]), y: tmpData.map(d => d[1]) }, { span: span, degree: degree } ); const x0 = { x: [ x_vec[index[i]] ] }; const traceContribution = f.predict(x0).fitted[0]; const y_actual = model.y[index[i]]; const e = y_actual - model.predict(x0).fitted[0]; const sseContribution = Math.pow(e,2); if(isNaN(traceContribution)){ console.log("bad trace contribution at observation " + i); console.log(tmpData[i]); } else { trace += traceContribution; sse += sseContribution; } tmpData[i][1] = 0.0; } return { span: span, gcv: (sse * (n/nindx))/(n*(Math.pow(1-trace/n,2))), sse: sse * n/nindx, trace: trace, evalTime: new Date() - startTime, M: M, n: n } }50:[{"id":"7f6974d40afc8454bb5f0c88f951b3500307002e31e174f48df67e3e88918f0498e545ce9a9a746188f67b819c493ceafebe85cff9963c205e1db54e5f56c828","href":"https://static.observableusercontent.com/files/7f6974d40afc8454bb5f0c88f951b3500307002e31e174f48df67e3e88918f0498e545ce9a9a746188f67b819c493ceafebe85cff9963c205e1db54e5f56c828","name":"pp_res_raw.csv","type":"text/csv","lastModified":1698455106193.429,"size":2768283}] 45:{"id":"538e9ef107239a3f","slug":null,"title":"Loess with optimal selection of span","thumbnail":null,"access_level":"public","create_time":"$D2023-10-28T00:58:28.971Z","update_time":"$D2024-05-26T05:19:34.618Z","view_time":null,"first_public_time":"$D2025-10-21T11:49:13.132Z","has_unpublished_edits":null,"fork":false,"stars":1,"starred":false,"owner":{"type":"team","id":"e453ea60d29ec23e","login":"simonjackman","name":"Simon Jackman","avatar_url":"https://avatars-next.observableusercontent.com/c14bb6609a47c8c937d2b9cfa6d34263984b5724ee02d5600987e325903bc377","tier":"starter_2024","bio":"","following":false},"version":786,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"da2057c4cea8c9e4","login":"simonjackman","name":"Simon Jackman","avatar_url":"https://avatars-next.observableusercontent.com/c14bb6609a47c8c937d2b9cfa6d34263984b5724ee02d5600987e325903bc377","tier":null,"bio":"","following":false}],"members":[{"type":"individual","id":"da2057c4cea8c9e4","login":"simonjackman","name":"Simon Jackman","avatar_url":"https://avatars-next.observableusercontent.com/c14bb6609a47c8c937d2b9cfa6d34263984b5724ee02d5600987e325903bc377","tier":null,"role":"editor"},{"type":"team","id":"e453ea60d29ec23e","login":"simonjackman","name":"Simon Jackman","avatar_url":"https://avatars-next.observableusercontent.com/c14bb6609a47c8c937d2b9cfa6d34263984b5724ee02d5600987e325903bc377","tier":"starter_2024","role":"viewer"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":1,"collection":null,"type":"notebook","body":{"titling":"auto","title":"Loess with optimal selection of span","cells":[{"id":0,"value":"# Loess with optimal selection of span","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":473,"value":"_loess_ (local regression) is a popular scatterplot smoother, a tool for visualizing the dependence in expectation between an outcome ${tex`y`} and a correlate ${tex`x`} (and here we restrict attention to the univariate case, smoothing ${tex`y`} against a single correlate at a time). \n\nA _hyper-parameter_ ${tex`\\lambda`} governs the _bandwidth_ or _span_ of smoother, and in turn, the smoothness of the resulting line overlaid on the ${tex`(x,y)`} scatterplot, the _loess_ approximation to the conditional mean function ${tex`E(y|x; \\lambda)`}.\n\nBut how to choose ${tex`\\lambda`}?","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":476,"value":"ToC()","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":74,"value":"$47","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":273,"value":"$48","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":282,"value":"\nviewof form = Inputs.form(\n {\n\n correlate: Inputs.select(\n correlateStrings,\n {\n label: \"Correlate:\",\n format: x => x.value, \n value: correlateStrings.find(t => t.name === \"ALP\")\n }\n ), \n\n theState: Inputs.select(\n [ \"NSW\", \"VIC\", \"QLD\", \"SA\", \"WA\", \"TAS\", \"ACT\", \"NT\", \"Australia\" ],\n {\n label: \"State/Territory:\"\n }\n ),\n \n span: Inputs.range([0.01,1], {\n label: \"span parameter for loess\",\n value: 0.25,\n step: 0.01\n }),\n \n degree: Inputs.select([1,2],{label: \"degree of local polynomial\"})\n },\n {\n template: (inputs) => htl.html`

\n ${Object.values(inputs)}\n

`\n }\n )","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":301,"value":"$49","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":336,"value":"$4a","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":427,"value":"The following timing experiments summarise the result of fitting a series of `loess` models with three different values of `span` by varying the number of evaluation points from 25 to 1,000. \n\n_This next cell may take a while to compute; over 10 seconds for me._","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":382,"value":"// timing experiments\nexperiments = async () => {\n const M = [25, 50, 100, 250, 500, 1000]\n const spans = [0.1, 0.5, 0.9]\n var experiments = [];\n\n for(const thisEval of M){\n for(const thisSpan of spans){\n \n var startTime = new Date(); \n\n const zzz = myLoess2(\n pp_res\n .filter(d => d.state === \"NSW\" & !isNaN(d.yes_per) &!isNaN(d.uni))\n .slice()\n .sort(function(x,y){return d3.ascending(x[\"uni\"],y[\"uni\"]);}),\n {\n x: \"uni\",\n y: \"yes_per\",\n span: Number(thisSpan),\n degree: 1,\n eval: Number(thisEval)\n }\n ) \n var endTime = new Date();\n var elapsedTime = endTime - startTime;\n\n experiments.push(\n {\n span: Number(thisSpan), \n eval: Number(thisEval),\n elapsedTime: elapsedTime\n }\n ); \n }\n }\n \nreturn experiments;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":434,"value":"Execution time is largely a function of the number of evaluation points, with high values of `span` noticably increasing compute time when the number of evaluation points is large.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":411,"value":"experiments()\n .then(\n (experiments) => \n Plot.plot({\n y: {\n label: \"Compute time (millisceconds)\"\n },\n marks: [\n Plot.line(experiments,\n {\n x: \"eval\",\n y: \"elapsedTime\",\n z: \"span\",\n tip: true\n })\n ]\n })\n)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":487,"value":"We now consider the statistical theory of choosing an optimal span parameter for `loess`, keeping in mind that for larger data sets we will be needing some computational short cuts.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":270,"value":"$4b","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":98,"value":"\\sum_{i=1}^n (y_i - \\hat{f}_{-i}(x_i; \\lambda))^2","pinned":false,"hidden":false,"mode":"tex","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":101,"value":"where ${tex`f_{-i}`} denotes the _loess_ fit that results from dropping observation ${tex`i`}, making clear the way we approach this as an out-of-sample prediction problem.\n\nA naive approach is to use brute-force computation to compute the \"leave-out-one\" sum of squared errors. This is not feasible for large data sets and where — as is the case with `javascript` visualisation libraries — we are driving the computation client-side to the user's hardware through a browser.\n\nFor the linear regression case, we can exploit properties of the celebrated ${tex`n`}-by-${tex`n`} _hat_ matrix ${tex`H(x)`} that projects ${tex`y`} into ${tex`\\hat{y} = H(x)y`}, so-called as it \"puts the hat\" on ${tex`y`}. In particular, for linear regression we wind up with the result that \n\n","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":123,"value":"\\sum_{i=1}^n (y_i - x_i' \\hat{\\beta}_{(-i}))^2 = \\sum_{i=1}^n \\frac{e_i^2}{(1-h_i)^2}","pinned":false,"hidden":false,"mode":"tex","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":120,"value":"$4c","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":588,"value":"\\sum_{i=1}^n (y_i - \\hat{f}_{-i}(x_i))^2 = \\sum_{i=1}^n \\frac{e_i^2}{(1 - \\text{trace}S_\\lambda/n)^2}","pinned":false,"hidden":false,"mode":"tex","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":159,"value":"$4d","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":553,"value":"### Example\n\nThe next cell produces estimates of ${tex`\\text{trace} S_\\lambda`} for our running example of the _loess_ smooth of Voice referendum Yes% on the proportion of voters with a university degree, across NSW polling places. We set ${tex`\\lambda`} = .15 for the purposes of testing. Estimating ${tex`\\text{trace} S_\\lambda`} by re-fitting _loess_ to all 3,034 data points takes over 8 seconds and exactly reproduces the value returned in `R` on my laptop in 20 milliseconds; forming the estimate with 50 evaluation points ","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":557,"value":"test_of_span()\n .then(\n (out) => \n Plot.plot({\n \n y: {\n label: \"Estimate of trace of smoother matrix\"\n },\n \n x: {\n ticks: 5,\n transform: (d) => Math.pow(10, d),\n type: \"log\",\n tickFormat: \"~s\",\n label: \"Number of evaluation points (log scale)\",\n },\n \n marks: [\n Plot.line(out,\n {\n x: d => Math.log10(d.npoints),\n y: \"z\",\n z: \"span\",\n tip: true\n }),\n\n Plot.ruleY([12.75779])\n ]\n })\n)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":178,"value":"$4e","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":585,"value":"## Grid search for GCV-optimal span\n\nWith some infrastructure for computing the trace of a _loess_ smoother matrix under at hand, we now turn back to the task at hand, selecting ${tex`\\lambda`} so as to minimise generalised cross-validation error. That is, our problem is to find","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":602,"value":"\\argmin_{\\lambda \\in (0,1]} \\sum_{i=1}^n \\frac{e_i^2}{(1 - \\text{trace}S_\\lambda/n)^2}","pinned":false,"hidden":false,"mode":"tex","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":608,"value":"where ${tex`e_i = y_i - \\hat{f}(x_i; \\lambda)`} are the residuals from the _loess_ fit resulting from `span` set to ${tex`\\lambda`}.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":612,"value":"We face a similar computation hurdle as before: computing the residuals ${tex`e_i`} is computationally intensive for large ${tex`n`}, as long as 8 or 9 seconds of clock time on a data set of 3,000 observations. ","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":613,"value":"$4f","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":634,"value":"loess_gcv_local(testdata, {x: \"uni\", y: \"yes_per\", eval: 50, span: .151, degree: 2})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":617,"value":"testdata = pp_res\n .filter(d => d.state==\"NSW\")\n .filter(d => !isNaN(d[\"uni\"]) & !isNaN(d[\"yes_per\"]))\n .slice()\n .sort(function(x,y){return d3.ascending(x[\"uni\"],y[\"uni\"]);})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":619,"value":"test_grid_search = {\n const spans = d3.ticks(.05,.5,45);\n var out = [];\n for(const theSpan of spans){\n const z = loess_gcv_local(testdata, {x: \"uni\", y: \"yes_per\", eval: 100, span: theSpan, degree: 2});\n out.push(z)\n }\n out.sort(function(x,y){return d3.ascending(x[\"gcv\"],y[\"gcv\"]);})\nreturn(out)\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":699,"value":"Plot.plot({\n marks: [\n Plot.auto(test_grid_search, {x: \"span\", y: \"gcv\"})\n ]\n})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":751,"value":"loess_optimal_bw(testdata, {x: \"uni\", y: \"yes_per\"})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":739,"value":"function loess_optimal_bw(\n data,\n {x,y,...options}\n){\n\n const degree = options.degree ?? 1;\n //console.log(\"x \" + x);\n //console.log(\"y \" + y);\n const spans = d3.ticks(.05,.95,10);\n var out = [];\n for(const theSpan of spans){\n const z = loess_gcv_local(data, {x, y, eval: 100, span: theSpan, degree: degree});\n console.log(z);\n out.push(z)\n }\n out.sort(function(x,y){return d3.ascending(x[\"gcv\"],y[\"gcv\"]);}) \n return(out[0].span);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":29,"value":"viewof theState = Inputs.select(\n [ \"NSW\", \"VIC\", \"QLD\", \"SA\", \"WA\", \"TAS\", \"ACT\", \"NT\" ],\n {\n label: \"State/Territory:\"\n }\n )","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":32,"value":"viewof correlate = Inputs.select(\n correlateStrings,\n {\n label: \"Correlate:\",\n format: x => x.value, \n value: correlateStrings.find(t => t.name === \"ALP\")\n }\n) ","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":34,"value":"Plot.plot({\n //grid: true,\n \n width: 800,\n\n y: {\n label: \"Yes %\",\n percent: true,\n domain: [0, 100]\n },\n\n x: {\n label: correlate.value,\n percent: correlate.name !== \"median_age\",\n },\n\n r: {\n range: [.5, 10]\n },\n \n \n marks: [\n //Plot.frame(),\n \n Plot.dot(\n pp_res.filter(d => d.state === theState), \n {\n x: correlate.name,\n y: \"yes_per\",\n r: \"formal\",\n strokeWidth: 0.5,\n alpha: 0.2\n }\n ),\n \n Plot.line(\n myLoess2(\n pp_res\n .filter(d => d.state === theState && isFinite(d[correlate.name]) && isFinite(d.yes_per)),\n {\n x: correlate.name,\n y: \"yes_per\"\n }\n ),\n {\n x: \"x\",\n y: \"y\",\n strokeWidth: 1.5,\n stroke: \"blue\"\n }\n )\n \n ]\n \n }\n) ","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":483,"value":"## Wrapping `Loess` for `Plot.line`\n\nThe following function simply wraps `Loess` with a key option being the `M` argument, the number of evaluation points. An object or arrays with named elements is returned, which can be passed to `Plot.line`.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":359,"value":"function myLoess2(data,{x,y,...options}){\n \n const degree = options.degree ?? 1;\n const M = options.eval ?? 50;\n\n const theSpan = options.span ?? loess_optimal_bw(data, {x, y, eval: M, degree: degree});\n\n const model = new Loess(\n {\n x: data.map(d => d[x]),\n y: data.map(d => d[y])\n },\n {\n span: theSpan,\n degree: degree\n }\n );\n\n const x0 = model.grid([ M ]);\n \n const fit = model.predict(x0).fitted;\n\n var out = [];\n for(let i = 0; i < M; i++){\n out.push({x: x0.x[i], y: fit[i] });\n }\n return(out)\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":783,"value":"myLoess2(\n pp_res\n .filter(d => d.state === theState && isFinite(d[correlate.name]) && isFinite(d.yes_per)),\n {\n x: correlate.name,\n y: \"yes_per\"\n }\n )","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":37,"value":"pp_res = await FileAttachment(\"pp_res_raw.csv\").csv({typed: true})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":31,"value":"correlateStrings = [\n { name: \"ALP\", value: \"Labor 2022 1st preferences (%)\" } ,\n { name: \"Coalition\", value: \"Coalition 2022 1st preferences (%)\" },\n { name: \"GRN\", value: \"Greens 2022 1st preferences (%)\" },\n { name: \"UAP_ON\", value: \"UAP + One Nation 2022 1st preferences (%)\" },\n { name: \"Teal\", value: \"Teal 2022 1st preferences (%)\" },\n { name: \"uni\", value: \"Tertiary degree, 2021 Census (%)\" },\n { name: \"indig\", value: \"Indigenous, 2021 Census (%)\" },\n { name: \"english_only\", value: \"English only, 2021 Census (%)\" },\n { name: \"median_age\", value: \"Median age, 2021 Census\" }\n ]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":744,"value":"import {loess_best, loess_grid, loess_trace, loess_gcv} from \"@simonjackman/loess-autobandwidth\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":64,"value":"//import {addLoessRegression} from \"@simonjackman/loess-autobandwidth\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":3,"value":"// Plot = addLoessRegression(await require(\"@observablehq/plot\"))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":187,"value":"//d3 = require(\"d3@7\", \"d3-regression@1\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":196,"value":"Loess = (await import(\"https://cdn.skypack.dev/loess\")).default","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":198,"value":"myLoess = ({x, y, ...options }) => {\n const span = options.span ?? 0.5;\n const degree = options.degree ?? 1;\n const M = options.eval ?? 50;\n\n const [X, setX] = Plot.column(x);\n const [Y, setY] = Plot.column(y); \n\n return {\n ...Plot.transform(options, function(data, facets) {\n const X = setX(Plot.valueof(data, x));\n const Y = setY(Plot.valueof(data, y));\n\n const model = new Loess(\n {\n x: X,\n y: Y\n },\n {\n span: span,\n degree: degree\n }\n );\n\n const x0 = model.grid([ M ]);\n\n const fit = model.predict(x0);\n\n var out = [];\n for(let i = 0; i < M; i++){\n out.push({x: x0[i], y: fit[i] });\n }\n \n return { out, facets };\n })\n }\n};","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":295,"value":"correlateFormats = (\n {\n \"ALP\": d3.format(\".1%\"),\n \"Coalition\": d3.format(\".1%\"),\n \"GRN\": d3.format(\".1%\"),\n \"UAP_ON\": d3.format(\".1%\"),\n \"Teal\": d3.format(\".1%\"),\n \"uni\": d3.format(\".1%\"),\n \"indig\": d3.format(\".1%\"),\n \"english_only\": d3.format(\".1%\"),\n \"median_age\": d3.format(\".0\")\n }\n)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":469,"value":"import {ToC} from \"@rmw4269/contents\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W50","stdlib":"1"}} 46:"$45:body"