diff --git a/examples/spatial/nyc_bym.ipynb b/examples/spatial/nyc_bym.ipynb
index 28377c4c4..baddd3a72 100644
--- a/examples/spatial/nyc_bym.ipynb
+++ b/examples/spatial/nyc_bym.ipynb
@@ -838,15 +838,36 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 6,
    "id": "8f8459d1-648d-48d6-ab24-1ec030dea611",
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "ename": "ModuleNotFoundError",
+     "evalue": "No module named 'pymc'",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
+      "Cell \u001b[1;32mIn[6], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mpymc\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpm\u001b[39;00m \u001b[38;5;66;03m# type: ignore\u001b[39;00m\n\u001b[0;32m      2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mnp\u001b[39;00m\n\u001b[0;32m      4\u001b[0m \u001b[38;5;66;03m# Assuming BYM_model is the model you have defined\u001b[39;00m\n\u001b[0;32m      5\u001b[0m \n\u001b[0;32m      6\u001b[0m \u001b[38;5;66;03m# Wrapping mixture and mu in pm.Deterministic in the model definition\u001b[39;00m\n",
+      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'pymc'"
+     ]
+    }
+   ],
    "source": [
-    "phi_pred = idata.posterior.phi.mean((\"chain\", \"draw\")).values\n",
-    "beta0_pred = idata.posterior.beta0.mean((\"chain\", \"draw\")).values\n",
-    "sigma_pred = idata.posterior.sigma.mean((\"chain\", \"draw\")).values\n",
-    "y_predict = np.exp(log_E + beta0_pred + sigma_pred * (1 / scaling_factor) * phi_pred)"
+    "import pymc as pm\n",
+    "import numpy as np\n",
+    "# Wrapping mixture and mu in pm.Deterministic in the model definition\n",
+    "with BYM_model:\n",
+    "    mixture = pm.Deterministic('mixture', some_mixture_expression)\n",
+    "    mu = pm.Deterministic('mu', some_mu_expression)\n",
+    "\n",
+    "# Use pm.do to condition on rho=1 and sample posterior predictive\n",
+    "with pm.do(BYM_model, {'rho': 1.0}):\n",
+    "    y_predict_rho_1 = pm.sample_posterior_predictive(idata, var_names=['mixture', 'mu'], predictions=True, extend_inferencedata=False)\n",
+    "\n",
+    "# Compute the mean of the predictions\n",
+    "y_predict = y_predict_rho_1.predictions.mu.mean(dim=['chain', 'draw'])\n"
    ]
   },
   {
@@ -1265,7 +1286,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.12.0"
   },
   "myst": {
    "substitutions": {
diff --git a/examples/spatial/nyc_bym.myst.md b/examples/spatial/nyc_bym.myst.md
index 5947f79f7..1a26e90f2 100644
--- a/examples/spatial/nyc_bym.myst.md
+++ b/examples/spatial/nyc_bym.myst.md
@@ -383,10 +383,20 @@ The payoff of all this work is that we can now visualize what it means to decomp
 We'll extract the means of several parameters to generate predictions. In the first case, we'll visualize only the predictions that come from the spatial component of the model. In other words, we are assuming $\rho = 1$ and we ignore $\theta$ and social fragmentation.
 
 ```{code-cell} ipython3
-phi_pred = idata.posterior.phi.mean(("chain", "draw")).values
-beta0_pred = idata.posterior.beta0.mean(("chain", "draw")).values
-sigma_pred = idata.posterior.sigma.mean(("chain", "draw")).values
-y_predict = np.exp(log_E + beta0_pred + sigma_pred * (1 / scaling_factor) * phi_pred)
+import pymc as pm
+import numpy as np
+# Wrapping mixture and mu in pm.Deterministic in the model definition
+with BYM_model:
+    mixture = pm.Deterministic('mixture', some_mixture_expression)
+    mu = pm.Deterministic('mu', some_mu_expression)
+
+# Use pm.do to condition on rho=1 and sample posterior predictive
+with pm.do(BYM_model, {'rho': 1.0}):
+    y_predict_rho_1 = pm.sample_posterior_predictive(idata, var_names=['mixture', 'mu'], predictions=True, extend_inferencedata=False)
+
+# Compute the mean of the predictions
+y_predict = y_predict_rho_1.predictions.mu.mean(dim=['chain', 'draw'])
+
 ```
 
 Then we'll overlay our predictions onto the same {ref}`adjacency map we built earlier <adjacency-map>`.