[2]:
import numpy as npimport polars as plimport numbaimport bebi103import bokehbokeh.io.output_notebook()
"}}; function display_loaded(error = null) { const el = document.getElementById("a2a7efa2-9d5b-4385-90cc-56f72bdcbd86"); if (el != null) { const html = (() => { if (typeof root.Bokeh === "undefined") { if (error == null) { return "BokehJS is loading ..."; } else { return "BokehJS failed to load."; } } else { const prefix = `BokehJS ${root.Bokeh.version}`; if (error == null) { return `${prefix} successfully loaded.`; } else { return `${prefix} encountered errors while loading and may not function as expected.`; } } })(); el.innerHTML = html; if (error != null) { const wrapper = document.createElement("div"); wrapper.style.overflow = "auto"; wrapper.style.height = "5em"; wrapper.style.resize = "vertical"; const content = document.createElement("div"); content.style.fontFamily = "monospace"; content.style.whiteSpace = "pre-wrap"; content.style.backgroundColor = "rgb(255, 221, 221)"; content.textContent = error.stack ?? error.toString(); wrapper.append(content); el.append(wrapper); } } else if (Date.now() < root._bokeh_timeout) { setTimeout(() => display_loaded(error), 100); } } function run_callbacks() { try { root._bokeh_onload_callbacks.forEach(function(callback) { if (callback != null) callback(); }); } finally { delete root._bokeh_onload_callbacks } console.debug("Bokeh: all callbacks have finished"); } function load_libs(css_urls, js_urls, callback) { if (css_urls == null) css_urls = []; if (js_urls == null) js_urls = []; root._bokeh_onload_callbacks.push(callback); if (root._bokeh_is_loading > 0) { console.debug("Bokeh: BokehJS is being loaded, scheduling callback at", now()); return null; } if (js_urls == null || js_urls.length === 0) { run_callbacks(); return null; } console.debug("Bokeh: BokehJS not loaded, scheduling load and callback at", now()); root._bokeh_is_loading = css_urls.length + js_urls.length; function on_load() { root._bokeh_is_loading--; if (root._bokeh_is_loading === 0) { console.debug("Bokeh: all BokehJS libraries/stylesheets loaded"); run_callbacks() } } function on_error(url) { console.error("failed to load " + url); } for (let i = 0; i < css_urls.length; i++) { const url = css_urls[i]; const element = document.createElement("link"); element.onload = on_load; element.onerror = on_error.bind(null, url); element.rel = "stylesheet"; element.type = "text/css"; element.href = url; console.debug("Bokeh: injecting link tag for BokehJS stylesheet: ", url); document.body.appendChild(element); } for (let i = 0; i < js_urls.length; i++) { const url = js_urls[i]; const element = document.createElement('script'); element.onload = on_load; element.onerror = on_error.bind(null, url); element.async = false; element.src = url; console.debug("Bokeh: injecting script tag for BokehJS library: ", url); document.head.appendChild(element); } }; function inject_raw_css(css) { const element = document.createElement("style"); element.appendChild(document.createTextNode(css)); document.body.appendChild(element); } const js_urls = ["https://cdn.bokeh.org/bokeh/release/bokeh-3.4.1.min.js", "https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.4.1.min.js", "https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.4.1.min.js", "https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.4.1.min.js", "https://cdn.bokeh.org/bokeh/release/bokeh-mathjax-3.4.1.min.js", "https://unpkg.com/@holoviz/panel@1.4.4/dist/panel.min.js"]; const css_urls = []; const inline_js = [ function(Bokeh) { Bokeh.set_log_level("info"); },function(Bokeh) { } ]; function run_inline_js() { if (root.Bokeh !== undefined || force === true) { try { for (let i = 0; i < inline_js.length; i++) { inline_js[i].call(root, root.Bokeh); } } catch (error) {display_loaded(error);throw error; }if (force === true) { display_loaded(); }} else if (Date.now() < root._bokeh_timeout) { setTimeout(run_inline_js, 100); } else if (!root._bokeh_failed_load) { console.log("Bokeh: BokehJS failed to load within specified timeout."); root._bokeh_failed_load = true; } else if (force !== true) { const cell = $(document.getElementById("a2a7efa2-9d5b-4385-90cc-56f72bdcbd86")).parents('.cell').data().cell; cell.output_area.append_execute_result(NB_LOAD_WARNING) } } if (root._bokeh_is_loading === 0) { console.debug("Bokeh: BokehJS loaded, going straight to plotting"); run_inline_js(); } else { load_libs(css_urls, js_urls, function() { console.debug("Bokeh: BokehJS plotting callback run at", now()); run_inline_js(); }); }}(window));
We continue our analysis of the drone sperm quality data set. Let’s load it and remind ourselves of the content.
[3]:
df = pl.read_csv( os.path.join(data_path, "bee_sperm.csv"), null_values=["NO SPERM", "No Sperm"], comment_prefix="#",)df.head()
[3]:
shape: (5, 18)
| Specimen | Treatment | Environment | TreatmentNCSS | Sample ID | Colony | Cage | Sample | Sperm Volume per 500 ul | Quantity | ViabilityRaw (%) | Quality | Age (d) | Infertil | AliveSperm | Quantity Millions | Alive Sperm Millions | Dead Sperm Millions |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| i64 | str | str | i64 | str | i64 | i64 | i64 | i64 | i64 | f64 | f64 | i64 | i64 | i64 | f64 | f64 | f64 |
| 227 | "Control" | "Cage" | 1 | "C2-1-1" | 2 | 1 | 1 | 2150000 | 2150000 | 96.726381 | 96.726381 | 14 | 0 | 2079617 | 2.15 | 2.079617 | 0.070383 |
| 228 | "Control" | "Cage" | 1 | "C2-1-2" | 2 | 1 | 2 | 2287500 | 2287500 | 96.349808 | 96.349808 | 14 | 0 | 2204001 | 2.2875 | 2.204001 | 0.083499 |
| 229 | "Control" | "Cage" | 1 | "C2-1-3" | 2 | 1 | 3 | 87500 | 87500 | 98.75 | 98.75 | 14 | 0 | 86406 | 0.0875 | 0.086406 | 0.001094 |
| 230 | "Control" | "Cage" | 1 | "C2-1-4" | 2 | 1 | 4 | 1875000 | 1875000 | 93.287421 | 93.287421 | 14 | 0 | 1749139 | 1.875 | 1.749139 | 0.125861 |
| 231 | "Control" | "Cage" | 1 | "C2-1-5" | 2 | 1 | 5 | 1587500 | 1587500 | 97.792506 | 97.792506 | 14 | 0 | 1552456 | 1.5875 | 1.552456 | 0.035044 |
Correlation
We might wish to investigate how two measured quantities are correlated. For example, if the number of dead sperm and the number of alive sperm are closely correlated, this would mean that a given drone produces some quantity of sperm and some fraction tend to be dead. Let’s take a look at this.
[4]:
# Set up plot on log scalep = bokeh.plotting.figure( frame_height=300, frame_width=300, x_axis_label="alive sperm (millions)", y_axis_label="dead sperm (millions)", x_axis_type="log", y_axis_type="log",)# Only use values greater than zero for log scaleinds = (pl.col("Alive Sperm Millions") > 0) & (pl.col("Dead Sperm Millions") > 0)# Populate glyphsfor color, ((treatment,), sub_df) in zip( bokeh.palettes.Category10_3, df.filter(inds).group_by("Treatment")): p.scatter( source=sub_df.to_dict(), x="Alive Sperm Millions", y="Dead Sperm Millions", color=color, legend_label=treatment, )p.legend.location = "bottom_right"bokeh.io.show(p)
There seems to be some correlation (on a log scale), but it is difficult to tell. We can compute the correlation with the bivariate correlation coefficient, also known as the Pearson correlation. It is the plug-in estimate of the correlation between variables (in this case alive and dead sperm). The correlation is the covariance divided by the geometric mean of the individual variances
The bivariate correlation coefficient is implemented with np.corrcoef(), but we will code our own and JIT it for speed.
[5]:
@numba.njitdef bivariate_r(x, y): """ Compute plug-in estimate for the bivariate correlation coefficient. """ return ( np.sum((x - np.mean(x)) * (y - np.mean(y))) / np.std(x) / np.std(y) / np.sqrt(len(x)) / np.sqrt(len(y)) )
We can use it to compute the bivariate correlation coefficient for the logarithm of alive and dead sperm.
[6]:
bivariate_r( np.log(df.filter(inds).get_column("Alive Sperm Millions").to_numpy()), np.log(df.filter(inds).get_column("Dead Sperm Millions").to_numpy()),)
[6]:
0.5219944217488051
Pairs bootstrap confidence intervals
How can we get a confidence interval on a correlation coefficient? We can again apply the bootstrap, but this time, the replicate is a pair of data, in this case a dead sperm count/alive sperm count pair. The process of drawing pairs of data points from an experiment and then computing bootstrap replicates from them is called pairs bootstrap. Let’s code it up for this example with the bivariate correlation.
Our strategy in coding up the pairs bootstrap is to draw bootstrap samples of the indices of measurement and use those indices to select the pairs.
[7]:
@numba.njitdef draw_bs_sample(data): """Draw a bootstrap sample from a 1D data set.""" return np.random.choice(data, size=len(data))@numba.njitdef draw_bs_pairs(x, y): """Draw a pairs bootstrap sample.""" inds = np.arange(len(x)) bs_inds = draw_bs_sample(inds) return x[bs_inds], y[bs_inds]
With our pairs sampling function in place, we can write a function to compute replicates.
[8]:
@numba.njit(parallel=True)def draw_bs_pairs_reps_bivariate(x, y, size=1): """ Draw bootstrap pairs replicates. """ out = np.empty(size) for i in numba.prange(size): out[i] = bivariate_r(*draw_bs_pairs(x, y)) return out
Finally, we can put it all together to compute confidence intervals on the correlation. To start, we extract all of the relevant measurements as Numpy arrays to allow for faster resampling (and that’s what our Numba’d functions require).
[9]:
# Extract NumPy arrays (only use values greater than zero for logs)inds = (pl.col("Alive Sperm Millions") > 0) & (pl.col("Dead Sperm Millions") > 0)dfs = df.filter(inds).partition_by('Treatment', as_dict=True)alive_ctrl = dfs[('Control',)].get_column("Alive Sperm Millions").to_numpy()alive_pest = dfs[('Pesticide',)].get_column("Alive Sperm Millions").to_numpy()dead_ctrl = dfs[('Control',)].get_column("Dead Sperm Millions").to_numpy()dead_pest = dfs[('Pesticide',)].get_column("Dead Sperm Millions").to_numpy()
Now we can compute the bootstrap replicates using our draw_bs_pairs_reps_bivariate() function.
[10]:
# Get repsbs_reps_ctrl = draw_bs_pairs_reps_bivariate( np.log(alive_ctrl), np.log(dead_ctrl), size=10_000)bs_reps_pest = draw_bs_pairs_reps_bivariate( np.log(alive_pest), np.log(dead_pest), size=10_000)
And from the replicates, we can compute and print the 95% confidence interval.
[11]:
# Get the confidence intervalsconf_int_ctrl = np.percentile(bs_reps_ctrl, [2.5, 97.5])conf_int_pest = np.percentile(bs_reps_pest, [2.5, 97.5])# Plot confidence intervalssummaries = [ dict( label="control", estimate=bivariate_r(np.log(alive_ctrl), np.log(dead_ctrl)), conf_int=conf_int_ctrl, ), dict( label="treated", estimate=bivariate_r(np.log(alive_pest), np.log(dead_pest)), conf_int=conf_int_pest, ),]bokeh.io.show( bebi103.viz.confints(summaries, x_axis_label="bivariate correlation of logs"))
We see a clear correlation in both samples, with a wide, but positive, confidence interval. Note that we did this analysis on a log scale, since the data span several orders of magnitude.
Computing environment
[12]:
%load_ext watermark%watermark -v -p numpy,polars,numba,bokeh,bebi103,jupyterlab
Python implementation: CPythonPython version : 3.12.4IPython version : 8.25.0numpy : 1.26.4polars : 1.1.0numba : 0.59.1bokeh : 3.4.1bebi103 : 0.1.21jupyterlab: 4.0.13
