update

.pdm-python

@@ -0,0 +1 @@
+/Users/assuncaolfi/Projects/site/.venv/bin/python

_extensions/shafayetShafee/bsicons/_extension.yml

@@ -5,3 +5,4 @@ quarto-required: ">=1.2.0"
 contributes:
   shortcodes:
     - bsicons.lua
+

_extensions/shafayetShafee/bsicons/bsicons.lua

@@ -19,29 +19,29 @@ return {
     local color = str(kwargs["color"])
     local label = str(kwargs["label"])
     local class = str(kwargs["class"])
-
+    
     if not isEmpty(size) then
       size = "font-size: " .. size .. ";"
     else
       size = ''
     end
-
+    
     if not isEmpty(color) then
       color = "color: " .. color  .. ";"
     else
       color = ''
     end
-
+    
     local style = "style=\"" .. size .. color .. "\""
-
+    
     if not isEmpty(label) then
       label = " aria-label=\"" .. label  .. "\""
     end
-
+    
     if isEmpty(class) then
       class = ''
     end
-
+    
     local role = "role=\"img\""
     local aria_hidden = "aria-hidden=\"true\""
 

_freeze/blog/aging-curve/index/execute-results/html.json

@@ -0,0 +1,16 @@
+{
+  "hash": "98d7afe0f115e118baa265ed5b35e599",
+  "result": {
+    "engine": "jupyter",
+    "markdown": "---\ntitle: Additive aging curve (draft)\ndescription: At what age does working memory peak?\ndate: today\n---\n\n::: {.callout-warning}\nThis post is a work in progress.\n:::\n\n\nRecently, I was involved in designing an experiment where each participant\nreceived a treatment at a random time $t$, between 5 and 30 minutes. After the \ntreatment, each participant produced a binary response. Soon, we realized time \nhad more than one effect over the response rate: as $t$ increased, the rate of\npositive responses 1) increased; then 2) plateaued; and finally 3) decreased.\n\nThis kind of non-monotonic relationship is common in cognitive and sports research, \nparticularly in the relationship between age and performance, where it's \ncalled an aging curve. For an interesting review of aging curves, see [@Vaci2019],\nwhere the authors discuss modeling strategies and study the effect of aging over\nthe performance of NBA players.\n\nAndrew Gelman wrote about this topic a couple of times in his blog: see his posts\nfrom [2018](https://statmodeling.stat.columbia.edu/2018/09/07/bothered-non-monotonicity-heres-one-quick-trick-make-happy/)\nand [2023](https://statmodeling.stat.columbia.edu/2023/01/01/how-to-model-a-non-monotonic-relation/), \nwhere he suggests modeling these relationships using an additive function like\n\n$$g(t) = g_1(t) + g_2(t),$$\n\nwhere  \n$g_1(t)$ is a monotonically increasing function with a right asymptote; and  \n$g_2(t)$ is a monotonically decreasing function with a left asymptote.\n\nIn this post, I'll analyze an experimental dataset by fitting and comparing\ntwo different models: a non-parametric bootstrap and a decomposable curve \nlike $g(t)$.\n\n## The Digit Span test\n\nThe motivation for Gelman's post from 2018 was a study relating age to peak\ncognitive functioning [@Hartshorne2015]. According to the study, some of their\nexperiments were conducted through a large scale online experimentation platform:\n\n> Participants in Experiment 2 (N = 10,394; age range = 10–69 years old)\n> [...] were visitors to TestMyBrain.org, who took part in experiments in\n> order to contribute to scientific research and in exchange for performance-\n> related feedback. [...] We continued data collection for each experiment for\n> approximately 1 year, sufficient to obtain around 10,000 participants, which\n> allowed fine-grained age-of-peak-performance analysis.\n\nThe data produced by this experiment is available online [@Germine_Hartshorne_2016]. \nThis dataset contains results for multiple tests, but I'll focus on the Digit Span \ntest during this analysis. According to [Cambridge Cognition](https://cambridgecognition.com/digit-span-dgs/):\n\n> Digit Span (DGS) is a measure of verbal short term and working memory that can be used in two formats, Forward Digit Span and Reverse Digit Span. This is a verbal task, with stimuli presented auditorily, and responses spoken by the participant and scored automatically by the software. Participants are presented with a random series of digits, and are asked to repeat them in either the order presented (forward span) or in reverse order (backwards span). While superficially very similar tasks, forward and backwards span rely on somewhat separable cognitive capacities: the simpler forward span task requires verbal working memory and attention, while the backwards span task additionally tests cognitive control and executive function.\n\nParticipants are scored according to their longest correctly repeated list of digits.\n\n::: {#digit-span .cell execution_count=2}\n``` {.python .cell-code}\nimport polars as pl\n\ndigit_span = (\n    pl.read_csv(\"data/experiment-2.csv\")\n    .filter(pl.col(\"age\").is_between(10, 69))\n    .with_columns(\n        y=(pl.col(\"DigitSpan\") - pl.col(\"DigitSpan\").mean()) / pl.col(\"DigitSpan\").std()\n    )\n)\n```\n:::\n\n\nThe relationship between age and Digit Span performance for each participant is plotted below:\n\n::: {#cell-digit-span-plot .cell execution_count=3}\n\n::: {.cell-output .cell-output-display execution_count=3}\n![](index_files/figure-html/digit-span-plot-output-1.png){#digit-span-plot width=599 height=445}\n:::\n:::\n\n\nVisually, it's still unclear if this relationship follows an aging curve, but\nwe'll get back to this matter in the next section.\n\n## Bootstrap estimates\n\nIn the original paper, the authors describe a bootstrap resampling procedure\nto estimate the distribution of ages of peak performance:\n\n> Estimates and standard errors for age of peak performance were calculated using\n> a bootstrap resampling procedure identical to the one used in Experiment 1\n> but applied to raw performance data. To dampen noise, we smoothed means for each\n> age using a moving 3-year window prior to identifying age of peak performance\n> in each sample. Other methods of dampening noise provide similar results.\n\nLet's decompose this method (as I understand it) into steps:\n\n1. Sample, with replacement, $n$ observations from the dataset;\n2. Calculate the mean performance for each age within the sample;\n3. Repeat steps 1 and 2 $m$ times;\n4. Sort each sample by age and smooth age means using a 3-year rolling average;\n5. Find the age of peak performance for each sample.\n\n::: {#bootstrap .cell execution_count=4}\n``` {.python .cell-code}\ndef sample_bootstrap(data: pl.DataFrame):\n    samples = (\n        data.sample(n * m, with_replacement=True, seed=seed)\n        .with_columns(sample=pl.arange(1, n * m + 1) % m)\n        .group_by(\"sample\", \"age\")\n        .agg(mean=pl.col(\"y\").mean())\n        .sort(\"sample\", \"age\")\n        .with_columns(smoothed_mean=pl.col(\"mean\").rolling_mean(3).over(\"sample\"))\n    )\n    peak = samples.group_by(\"sample\").agg(\n        age=pl.col(\"age\").get(pl.col(\"smoothed_mean\").arg_max())\n    )\n    return samples, peak\n\n\nn = digit_span.height\nm = 10000\nseed = 37\nsamples, peak = sample_bootstrap(digit_span)\n```\n:::\n\n\nThis algorithm yields the following bootstrap distribution of ages of peak performance:\n\n::: {#cell-bootstrap-distribution .cell execution_count=5}\n\n::: {.cell-output .cell-output-display execution_count=5}\n![](index_files/figure-html/bootstrap-distribution-output-1.png){#bootstrap-distribution width=599 height=445}\n:::\n:::\n\n\nThis distribution suggests two important things:\n\n1. The most probable age of peak performance is 33;\n2. Peak performance could happen anywhere between the early 20s and late 30s, except during the late 20s.\n\nSuggestion 2 is probably not true. In fact, this distribution seems like a mixture of two distributions, but I'll get back to this point in the next section. For now, I'll use our bootstrap estimates to replicate figure 3a from the original paper. Using the samples obtained in step 4, for each age mean, I calculated its median and 90% interquantile range, yielding a nonparametric curve:\n\n::: {#cell-bootstrap-curve .cell execution_count=6}\n\n::: {.cell-output .cell-output-display execution_count=6}\n![](index_files/figure-html/bootstrap-curve-output-1.png){#bootstrap-curve width=599 height=445}\n:::\n:::\n\n\nSince this curve is empirical, there's not much more than meets the eye here. However, note that it follows the rising, plateauing and falling behavior of an aging curve. There's a steep increase during ages 10 to 20, followed by a plateau between 20 and 30, and a slow decline beginning at 40.\n\n### The language effect\n\n::: {#cell-performance-language .cell execution_count=7}\n\n::: {.cell-output .cell-output-display execution_count=7}\n![](index_files/figure-html/performance-language-output-1.png){#performance-language width=727 height=445}\n:::\n:::\n\n\n::: {#cell-english-span .cell execution_count=8}\n\n::: {.cell-output .cell-output-display execution_count=8}\n![](index_files/figure-html/english-span-output-1.png){#english-span width=599 height=445}\n:::\n:::\n\n\n## Additive functions\n\n$$\n\\begin{align}\ng(t) &= g_1(t) + g_2(t) \\\\\ny &\\sim \\mathrm{Normal}(g(t), \\sigma) \\\\\n\\end{align}\n$$\n\n### Double exponential\n\n$$\n\\begin{align}\ng_1(t) &= \\alpha + \\beta_1 \\exp(-\\lambda_1 t) \\\\\ng_2(t) &= \\beta_2 \\exp(\\lambda_2 t) \\\\\n\\end{align}\n$$\n\n::: {#cell-double-exponential .cell execution_count=9}\n``` {.python .cell-code}\nimport numpy as np\nimport pymc as pm\n\n\ndef g(x):\n    return g_1(x) + g_2(x)\n\n\ndef g_1(x):\n    return α + β[0] * pm.math.exp(-λ[0] * x)\n\n\ndef g_2(x):\n    return β[1] * pm.math.exp(λ[1] * x)\n\n\nage = english_span.get_column(\"age\")\ny = english_span.get_column(\"y\")\nage_range = np.arange(age.min(), age.max() + 1)\nwith pm.Model() as double_exponential:\n    t = pm.Data(\"t\", age)\n    α = pm.Normal(\"α\", 0, 1)\n    β = pm.Normal(\"β\", 0, 1, size=2)\n    λ = pm.HalfNormal(\"λ\", 0.004, size=2)\n    μ = pm.Deterministic(\"μ\", g(t))\n    σ = pm.HalfNormal(\"σ\", 1)\n    pm.Normal(\"y\", mu=μ, sigma=σ, observed=y)\n    curve = pm.Deterministic(\"curve\", g(age_range))\n    samples = pm.sample(progressbar=False, target_accept=0.95, random_seed=seed)\n```\n\n::: {.cell-output .cell-output-stderr}\n```\nWARNING:pytensor.tensor.blas:Using NumPy C-API based implementation for BLAS functions.\n```\n:::\n\n::: {#double-exponential .cell-output .cell-output-display}\n```{=html}\n<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">\n</pre>\n```\n:::\n:::\n\n\n::: {#cell-double-exponential-curve .cell execution_count=10}\n\n::: {.cell-output .cell-output-display execution_count=10}\n![](index_files/figure-html/double-exponential-curve-output-1.png){#double-exponential-curve width=599 height=445}\n:::\n:::\n\n\n::: {#cell-double-exponential-peak .cell execution_count=11}\n\n::: {.cell-output .cell-output-display execution_count=11}\n![](index_files/figure-html/double-exponential-peak-output-1.png){#double-exponential-peak width=599 height=445}\n:::\n:::\n\n\n### Double logistic\n\n[@Lipovetsky2010]\n\n$$\n\\begin{align}\ng_1(t) &= \\alpha_1 + \\frac{\\alpha_2 - \\alpha_1}{1 + \\exp(\\beta_1 - \\lambda_1 t)} \\\\\ng_2(t) &= \\frac{\\alpha_3 - \\alpha_2}{1 + \\exp(\\beta_2 + \\lambda_2 t)}\n\\end{align}\n$$\n\n::: {#cell-double-logistic .cell execution_count=12}\n``` {.python .cell-code}\ndef g_1(t):\n    return α[0] + (α[1] - α[0]) / (1 + pm.math.exp(β[0] - λ[0] * t))\n\n\ndef g_2(t):\n    return (α[2] - α[1]) / (1 + pm.math.exp(β[1] + λ[1] * t))\n\n\nwith pm.Model() as double_logistic:\n    t = pm.Data(\"t\", age)\n    α = pm.Normal(\"α\", 0, 1, size=3)\n    β = pm.Normal(\"β\", 0, 1, size=2)\n    λ = pm.HalfNormal(\"λ\", 1, size=2)\n    μ = pm.Deterministic(\"μ\", g(t))\n    σ = pm.HalfNormal(\"σ\", 1)\n    pm.Normal(\"y\", mu=μ, sigma=σ, observed=y)\n    curve = pm.Deterministic(\"curve\", g(age_range))\n    samples = pm.sample(progressbar=False, target_accept=0.95, random_seed=seed)\n```\n\n::: {#double-logistic .cell-output .cell-output-display}\n```{=html}\n<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\">\n</pre>\n```\n:::\n\n::: {.cell-output .cell-output-stderr}\n```\nERROR:pymc.stats.convergence:There were 1 divergences after tuning. Increase `target_accept` or reparameterize.\nERROR:pymc.stats.convergence:The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details\n```\n:::\n:::\n\n\n::: {#bc531661 .cell execution_count=13}\n\n::: {.cell-output .cell-output-display execution_count=13}\n![](index_files/figure-html/cell-14-output-1.png){width=599 height=445}\n:::\n:::\n\n\n::: {#21b54688 .cell execution_count=14}\n\n::: {.cell-output .cell-output-display execution_count=14}\n![](index_files/figure-html/cell-15-output-1.png){width=599 height=445}\n:::\n:::\n\n\n",
+    "supporting": [
+      "index_files"
+    ],
+    "filters": [],
+    "includes": {
+      "include-in-header": [
+        "<script src=\"https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js\" integrity=\"sha512-c3Nl8+7g4LMSTdrm621y7kf9v3SDPnhxLNhcjFJbKECVnmZHTdo+IRO05sNLTH/D3vA6u1X32ehoLC7WFVdheg==\" crossorigin=\"anonymous\"></script>\n<script src=\"https://cdnjs.cloudflare.com/ajax/libs/jquery/3.5.1/jquery.min.js\" integrity=\"sha512-bLT0Qm9VnAYZDflyKcBaQ2gg0hSYNQrJ8RilYldYQ1FxQYoCLtUjuuRuZo+fjqhx/qtq/1itJ0C2ejDxltZVFg==\" crossorigin=\"anonymous\" data-relocate-top=\"true\"></script>\n<script type=\"application/javascript\">define('jquery', [],function() {return window.jQuery;})</script>\n"
+      ]
+    }
+  }
+}