{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Correlation models and hypothesis testing\n", "This demo shows how to use the PCMPy toolbox to test hypotheses about the true (noisefree) correlation between activity patterns as described in Diedrichsen et al. (2026).\n", "We will show how to obtain individual and group-level maximum-likelihood estimates of the true correlation between activity patterns, and how to test hypotheses about the size of these correlations using a bootstrap approach. In the first part of this jupyter notebook, we will focus on the single-pattern setting, in which we want to estimate the correlation between **two** activity patterns. In the second part, we will consider a slightly more complex multi-pattern setting, in which we want to estimate the true correlation between two **sets** of activity patterns measured under two different conditions. \n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# First import necessary libraries\n", "#(make sure PcmPy is in your python path)\n", "import PcmPy as pcm\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from numpy import exp, sqrt\n", "import scipy.stats as sps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Single pattern setting: Estimating the correlation between two patterns\n", " The correlation between two activity patterns in fMRI can be quite low, even if the correspondence between these two patterns is perfect. This is because both patterns are measured with measurement noise. To solve this problem, we can obtain the maximum-likelihood estimator of the correlations, using repeated measurements of the same two activity patterns across multiple runs (repetitions) to estimate the signal and noise variance on each pattern. Note that this estimator is not unbiased, and can get unstable for very low signal to noise datasets - so we need to be caseful in the interpretation and inference (see Diedrichsen et al., 2026 for details).\n", "\n", "In the simplest case of 2 activity patterns, the model will have 3 model parameters: \n", "- $\\sigma^2_x$: The strength (variance across voxels) of the activity patterns associated with condition **X**. \n", "- $\\sigma^2_y$: The strength (variance across voxels) of the activity patterns associated with condition **Y**. \n", "- $r$: The true correlation between **A** and **B** \n", "\n", "And one additional noise parameter: \n", "- $\\sigma^2_{\\epsilon}$: The variance of the measurement noise across repeated measures of the same pattern. \n", "\n", "To ensure that the variance parameters are positive, and the correlation between $-1,1$, we fit the following nonlinear transforms:\n", "\n", "$\\theta_1 = log(\\sigma^2_x)$ \n", "\n", "$\\theta_2 = log(\\sigma^2_y)$\n", "\n", "$\\theta_3 = atanh(r)$ \n", "\n", "$\\theta_{\\epsilon} = log(\\sigma^2_{\\epsilon})$ \n", "\n", "These (hyper-)parameters are collectively referred to by $\\theta$ (thetas). The algorithm maximizes the Type-II maximum likelihood of the data given the parameters: $p(Y|\\theta)$ " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulating data\n", "In this example, we will use simulated data with a known correlation between the true (noiseless) activity patterns (`Mtrue`). In this example, we set that our two activity patterns are positively correlated with a true correlation of 0.7.\n", "\n", ">Note that in the single-pattern setting, we crate a `pcm.CorrelationModel` with `num_items` set to 1, as we have only 1 activity pattern per condition. We also set `cond_effect` to `False`, as we do not want to account for the overall condition effect.\n", "\n", "Because the correlation of the true mode (`Mtrue`) is fixedm the model has only 2 hyper-parameters reflecting the signal strength (or true pattern variance) for each activity pattern (`item`). \n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The true model has 2 hyper parameters\n" ] } ], "source": [ "Mtrue_sim = pcm.CorrelationModel('corr', num_items=1, corr=0.7, cond_effect=False)\n", "print(f\"The true model has {Mtrue_sim.n_param} hyper parameters\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now use the simulation module to create 20 datasets (e.g., one per simulated participant) with a relatively low signal-to-noise level (0.2:1). We will use a design with 2 conditions and 8 partitions/runs per dataset.\n", ">Note that the thetas are specified as log(variance)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Create the design. In this case it's 2 conditions, across 8 runs (partitions)\n", "n_cond=2\n", "n_part = 8\n", "cond_vec, part_vec = pcm.sim.make_design(n_cond, n_part)\n", "\n", "# Starting from the true model above, generate 30 datasets/participants\n", "# with relatively low signal-to-noise ratio (0.06:1)\n", "# Note that signal variance is drawn from a gamma distribution with shape=2 and scale=0.03 (mean 0.06)\n", "# The noise is by default set to 1.\n", "# For replicability, we are using a fixed random seed (100)\n", "n_subj = 30\n", "rng = np.random.default_rng(seed=102)\n", "signal = np.random.gamma(shape=2,scale=0.03,size=(n_subj,))\n", "D = pcm.sim.make_dataset(model=Mtrue_sim, \\\n", " theta=[0,0], \\\n", " cond_vec=cond_vec, \\\n", " part_vec=part_vec, \\\n", " n_sim=n_subj, \\\n", " signal=signal,\\\n", " rng=rng)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First let's look at the correlation that we get when we calculate the *naive* correlation—i.e. the correlation between the two estimated activity patterns. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated mean correlation: 0.2334\n" ] } ], "source": [ "def get_corr(X,cond_vec):\n", " \"\"\"\n", " Get normal correlation\n", " \"\"\"\n", " p1 = X[cond_vec==0,:].mean(axis=0)\n", " p2 = X[cond_vec==1,:].mean(axis=0)\n", " return np.corrcoef(p1,p2)[0,1]\n", "\n", "r = np.empty((20,))\n", "for i in range(20):\n", " data = D[i].measurements\n", " r[i] = get_corr(data, cond_vec)\n", "print(f'Estimated mean correlation: {r.mean():.4f}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can see, due to measurement noise, the estimated mean correlation is much lower than the true value of 0.7.\n", "\n", "This is not a problem if we just want to test the hypothesis that the true correlation is larger than zero. Then we can just calculate the individual correlations per subject and test them against zero using a t-test. \n", "\n", "However, if we want to test whether the true correlation is larger or smaller than a specific value (for example `true_corr=1`, indicating that the activity patterns are the same), or if we want to test whether the correlations are higher in one brain area than another, then this becomes an issue. \n", "\n", "Different brain regions measured with fMRI often differ dramatically in their signal-to-noise ratio. Thus, we need to take into account our level of measurement noise. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fitting models\n", "To get the maximum-likelihood estimate of the correlation of X and Y we generate a flexible model (`Mflex`) that has the correlation as a parameter that is being optimized. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The flexible model has 3 hyper parameters\n" ] } ], "source": [ "# Now make the flexible model\n", "Mflex_sim = pcm.CorrelationModel(\"flex\", num_items=1, corr=None, cond_effect=False)\n", "print(f\"The flexible model has {Mflex_sim.n_param} hyper parameters\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now fit the models to all datasets in one go. The resulting dataframe `T` has the log-likelihoods, noise variances, and iterations for the different dataset (rows). The second return argument `theta` contains the parameters for the model fit. \n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
variablelikelihoodnoiseiterations
modelflexflexflex
0-224.9288750.8675805.0
1-244.5788430.98819712.0
2-272.7725971.1031094.0
3-263.3164741.0854179.0
4-243.0561260.9572074.0
\n", "
" ], "text/plain": [ "variable likelihood noise iterations\n", "model flex flex flex\n", "0 -224.928875 0.867580 5.0\n", "1 -244.578843 0.988197 12.0\n", "2 -272.772597 1.103109 4.0\n", "3 -263.316474 1.085417 9.0\n", "4 -243.056126 0.957207 4.0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# In this case, we the mean of the block is not modelled\n", "# This is valid if the measurement noise for A and B within a block are independent.\n", "# We don't need to include a scale parameter given that we don't have a fixed model\n", "T, theta = pcm.fit_model_individ(D, Mflex_sim, fixed_effect=None, fit_scale=False, verbose=False)\n", "T.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also fit all the datasets together with one single model, allowing only for separate scale and noise parameter for each subject. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
variablelikelihoodnoisescaleiterations
modelflexflexflexflex
0-226.8877390.8762953.29794113
1-245.2936530.9934601.02944613
2-273.2121261.1137031.25339313
3-263.4730571.0874450.50321413
4-243.5976510.9666181.93799213
\n", "
" ], "text/plain": [ "variable likelihood noise scale iterations\n", "model flex flex flex flex\n", "0 -226.887739 0.876295 3.297941 13\n", "1 -245.293653 0.993460 1.029446 13\n", "2 -273.212126 1.113703 1.253393 13\n", "3 -263.473057 1.087445 0.503214 13\n", "4 -243.597651 0.966618 1.937992 13" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# We can either include or not include a individual scale factor.\n", "# Given that that the variance parameters (and correlation parameter) are the same across the entire group, but\n", "# since we cannot assume that the signal-to-noise ratio is the same across participants, we include a scale factor.\n", "T_gr, theta_gr = pcm.fit_model_group(D, Mflex_sim, fixed_effect=None, fit_scale=True, verbose=False)\n", "T_gr.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the likelhood of each of the datasets is a bit lower than for the individual fit. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpreting the individual fits\n", "We can now get the signal variance (theta[0], theta[1]), the correlation parameters (theta[2]), and the noise parameter (theta[3]) from the parameters. \n", "Note that all the variance parameter are internally represented as log(var) and the correlation parameters are fisherz(r), so we need to make sure we transform it. \n", "\n", "The overall estimated fSNR is $fSNR = \\sigma_1 \\sigma_2 N_{part}/\\sigma^2_{\\epsilon}$. \n", "\n", "Where $N_{part}$ is the number of partitions/runs per dataset.\n", "\n", "To make it easier, PCM has `get_correlation` and `get_fSNR` functions that do the transformation for us.\n", "\n", "As a diagnostic plot, we recommend plotting the individual correlation parameter estimates against the estimated fSNR. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Get the individual parameter estimates from the individual fits\n", "theta_ind = theta[0] # Get the parameters from the first (and only model)\n", "r_indiv = Mflex_sim.get_correlation(theta[0]) # Inverse fisher-z transform\n", "fSNR_indiv = Mflex_sim.get_fSNR(theta[0], n_part, separate=False)\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Add the group estimate as a horizontal line\n", "theta_g,_= pcm.group_to_individ_param(theta_gr[0],Mflex_sim,20)\n", "r_group = Mflex_sim.get_correlation(theta_g)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Correlation parameter')" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAGwCAYAAAC5ACFFAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWOBJREFUeJzt3Xl4U2X6N/Bvki5pCYSlS1qs7IoFBKG2FhxQqbSgLA6MICiLAsomUhg2hVpREHQEBYQB2fyBsviCA9UpSwEVWcoAVSqLA1RZ2rRAabrRLXnePzo9cGibJiVpkvb7ua5ckCf3Ob3PkpM75zx5jkIIIUBEREREFVI6OgEiIiIiZ8ZiiYiIiMgMFktEREREZrBYIiIiIjKDxRIRERGRGSyWiIiIiMxgsURERERkhpujE6gNTCYTUlNTUb9+fSgUCkenQ0RERBYQQiAnJweBgYFQKis/f8RiyQZSU1MRFBTk6DSIiIioGq5cuYIHHnig0tdZLNlA/fr1AZSu7AYNGjg4GyIiIrJEdnY2goKCpM/xyrBYsoGyS28NGjRgsURERORiqupCww7eRERERGawWCIiIiIyg8USERERkRksloiIiIjMYLFEREREZAaLJSIiIiIzWCwRERERmcFiiYiIiMgMFktEREREZnAEbxdgNAkkpmRCn12AGzmFuJVfCKVCifBWTfBEyyZQKS2/eW/ZvDJyCuBXX43QFo0tmr6609lqno7++/aev71zcTZ3L6+PxhMQwI28wjqx7FWpal+o7ftKbV++MnVlOWsLlyqWfvzxR3z00Uc4ceIE0tLSsGPHDgwYMMDsNAcPHkR0dDR+++03BAUF4Z133sHIkSNlMcuXL8dHH30EvV6Pjh07YunSpQgNDbXfglig7I2094we3yalIjOvqFzMsgMX0NDbHR/+tQOi2gdUOc/45DTE7jqDNEOB1BagVSOmb7DZ6as7na1ycfTft/f87Z2Ls6loee9Wm5e9KlXtC7V9X6nty1emrixnbaIQQghHJ2Gpf//73/j555/RpUsX/PWvf62yWEpJSUH79u3xxhtvYPTo0UhISMBbb72F7777DpGRkQCALVu2YPjw4Vi5ciXCwsKwZMkSbNu2DefPn4efn59FeWVnZ0Or1cJgMNjk3nBVfZhUZOXLnasseMZtPIl7N3bZ95gVlUxf3enMsWaejv771eHo5XNmlS3v3Wrrslelqn1hbPcWWPVjSq3dV+rKe6GuLKersPTz26WKpbspFIoqi6UZM2bgu+++Q3JystQ2ZMgQZGVlIT4+HgAQFhaGxx9/HMuWLQMAmEwmBAUFYdKkSZg5c6ZFuZSt7NTUVOh0OumGfEVFRSguLoabmxs8PT2l+Ly8PACAl5cXlMrSbmPFxcUoKipCwrkMTN52RnojmYpKCyaFuwcUitJYYSyBMJZAoVRC4eYBoPRbyZ5JT0CpANRqNVQqFQCgpKQE+bcL8OziH5Gef2dTm4oLAAEo3NyhVKqg06pxcGp3lBQXQalUwsvLC0aTwJML9+PaDQMgBBRu7lAoS+crTEagpBi6hl44/E5v6fTx7du3YTKZ4OnpCTe30hOXRqMRBQUFUCgU8FR74cmF+5FmKIAoKYIwmaBQuUGh+t9JTpMRfvWUSJj6NLy9vc3HChP8vBXYF/0UGtTXSMtWWFiIkpISuLu7w8PDQ9qut2/fhtEkELnsmFSIipJiCJMRCpUKCpU7FAD8G3hiz6QnoFIq4O3tXW573j1fIQTy8/MBAN7e3jAJ4MmF+5GamQNhNEKhVEHh5i7lJooK4K/1xM9v94ZSqbQ41t1NJdtPVCoV1Gq1FJufnw8hhGzbWxNbUlKCwsJCaduXKdue1sRasu1NxYUV7lOipBhQACp3NXRaNQ7NeAbFRYUwGo3w8PCAu7u7bHsCQL169aQcCgoKLI6taD+5e3taE1vRflLZ+76iWIVShZ5Lfpb2y4re9wpTCYwl8vd9WawCQECTBvh5VgRUSoVNtn119pOKtr0lsUXFJXhyfjz02YVQut/Jt+x9r1S5IaCxBodmPAMIk7RPeXt722zbV3c/seYY4an2wl8WHSg9phmLy73vFQB8vQT2RT8FTT1v6fOhOp8lrniMuHt7WhNrzba/N9bSYqlWd/A+cuQIIiIiZG2RkZE4cuQIgNId8MSJE7IYpVKJiIgIKaYihYWFyM7Olj0AIDAwEDdu3JDiPvroI2g0GkycOFE2vZ+fHzQaDS5fviy1LV++HBqNBq++Nlr2jePayldxZfEgFN+4IrXlnt6HK4sH4frORVJbmqEAD7V9BBqNBidPnpTat2zZAm2D+vh17duyHPQbonFl8SAUXv0N4n/TL167GRqNRlofiSmZSDMUIP2rmbiyeBBup9yZb8Gfv+Ly4kE4tXwiElMypfbevXtDo9Fgx44dUtvRo0eh0WjQsWNHaZ4AcH3HfFxZPAh5Zw5KsUXX/8Txd/uiZevWstgbcf/AlcWDkPtLvBRbfCsNx9/ti6ZNm8qW7fXXX4dGo8Gnn356Z/2kpUGj0aBx40ayM3aZ+7/AlcWDYDiyFQAgAKRm3IS2QX1oNBqUlJRIsW+//TY0Gg3efvvOuiwpKYFGo4FGo4HBYJByNhzZiiuLByFz/xey3C5/OhjH3+2Lfx87K8Vm/2cnriwehJu7l8lir3w+Asff7YsdP5yQ2latWgWNRoNXXnlFFtumTRtoNBqcPn1aatu0aRM0Gg0GDhwoi+3YsSM0Gg2OHj0qte3YsQMajQa9e/eWxYaHh0Oj0WD//v1S2+7du6HRaNC9e3dZbEREBDQaDeLi4qS2n376CRqNBo8//rh82/+/eaXb/uyPUmxR+kVcWTwIqV+Ml/bJxJRMDBkyBBqNBuvXr5diz549C41Gg+bNm8tyeO2116DRaLB8+fI76/zyZWg0mnJniidOnAiNRoOPPvpIartx44a0Pe82Y8YMaDQaxMbGSm35+flSbNmHIQDExsZCo9FgxowZsnmUxVZ0jBgycoxsv7y6bBiuLB4EY/Z1qc1w4jtcWTwIN/79qWy+11a+isuLB+Hypf9K78X169eXznfIEFlscHBwhccIjUaDfv36yWIff/xxaDQa/PTTT1JbXFyc7BhRpnv37tBoNNi9e7fUtn//fmg0GoSHh8tiKzpGrNuxB8ff7Yu0dZNksWXHiNwzB6X94fTp09BoNGjTpo0s9pVXXoFGo8GqVauktosXL0Kj0Vh1jGjYsKEsNjo6GhqNBvPnz5faDAaDtD2tOUbs//UPaTtXdIwQAI6/NwDaBvWRlpYmtX/66afQaDR4/fXXZbk1bdoUGo0GFy9elNpc+Rhxt379+kGj0WDLli1S28mTJ6HRaBAcHCyLvd9jhCVqdbGk1+vh7+8va/P390d2djZu376NGzduwGg0Vhij1+srne+CBQug1WqlR1BQkM1yvl1srPa0RtP9nSQ03Jb3i8rIsewyoKVx1sSaTMLyWCc6N2ppzjfyCiyOvZlbeD8pOQ1r9pP7mcYVFdzH+/5urrq+buWX75NZEVddvjLXc107/7qsVl+Ge+ihhzBq1CjMmjVLavv+++/x3HPPIT8/H7du3ULTpk1x+PBh2bef6dOn44cffsCxY8cqnG9hYSEKC+98gGVnZyMoKOi+L8Nt/88fmP7/ksudYgeqvgwHAOtefhShLRqXOx3649lUjFz/Hyjd7+Rw92W4sssg/zcqBJ0f0EinTo9cvImXVh+t4pKJAlvG90B4qyYAqj51+kvabby0uvTbSkWX1srmu+HVUKi9vM3HChNEcRHWj3ocz3R4ULZ9KjvFfuzSTby66c43q3svw5XOV2Ddy48irGUTqy/DHb2UiZdWH63wFPvd2/Prcd2hVCotju3WxlfaT1z1FPvd276qy3Bll2K+HvMEHmtar9ZfhvvPnwaM+PJUuW1vyfv+7tjNY7sivFUTl7sMd+j3DAxd+ZNs2wPl3/dfj3kCoc0buuxluF/SbmPoF8f+tz0rf9+vH/U4nmr3AC/DOdFlOJf6NZy1dDod0tPTZW3p6elo0KABvLy8oFKpoFKpKozR6XSVztfT01O2w5apV6+e9KYBAA8PD+kNc2/cvdzd3fGgX2PZQRAAlB7qcrGyfjv/E6BVo3vwA+V+eurm5oYe7YLQ1Oe/0BsKpEt8dx+QFAB0WjW6tvGTTR/aojECtGroDSjfGVGpgtKjtK9TaIvGUvvdb4wyKpVKWubQFl7/m2cB4OYBxT2xSqUKOt966B78gLRclcYqlND5NkSPdvIzexVtH6VSiXr16qFHO28EaO+sC4WbOxS4c6BSAAho6IUe7YLKrcuKtqdCoZBtT9k6U7njXiqP0r44T7TyuWv5LIsFSveTsoPB3e4+cFQn1s3NTTog3a2i7WlNbGXb/u7CvYxCqYLCo/RgW7ZPVvZz6rLtea+7D/hVxVa0n9y7PasTa837viz2ybbeCNCelfbLit73Kjc3CJVbufei0kMtW1+Abbb9/e4nd2/7qmLDW/uiqW9D2TEKABT/e9/fuz/YY9vbYj+p6hgR1tL7zjFN5S59QZNiATT93zFNqaz+Z4mrHiOqG2vNtq8o1hK1+jJceHg4EhISZG179+6VziJ5eHigS5cushiTyYSEhIRy19lrQtkHbXVG2ojpG1zpGB0qpQIxfUuv8d4bUfa8oumrO5051szT0X+/Ohy9fM7M3PLerTYue1Wq2hcUAMb8pUWlrwOuvb7qynuhrixnbeRSxVJubi6SkpKQlJQEoHRogKSkJKmj9KxZszB8+HAp/o033sClS5cwffp0nDt3Dp9//jm2bt2KKVOmSDHR0dFYvXo1NmzYgLNnz2LcuHHIy8vDqFGjanTZAMs/TO7WyNu9ymEDACCqfQBWvNwZOq28qtZp1WZ/qlrd6WyVi6P/fnU4evmcWWXLe7fauuxVqWpfmNUnuFbvK3XlvVBXlrO2cak+SwcPHsTTTz9drn3EiBFYv349Ro4ciT/++AMHDx6UTTNlyhScOXMGDzzwAObMmVNuUMply5ZJg1J26tQJn332GcLCwizOqybGWWpczx39OwYisKE3R/DmCN4ujyN4V44jeNfu5StTV5bT2dX6cZacia2LJYBvJCIiIntjB28Xp1IqpF+YERERkeO4VJ8lIiIioprGYomIiIjIDBZLRERERGawWCIiIiIyg8USERERkRksloiIiIjMYLFEREREZAaLJSIiIiIzWCwRERERmcFiiYiIiMgMFktEREREZrBYIiIiIjKDxRIRERGRGSyWiIiIiMxgsURERERkBoslIiIiIjNYLBERERGZwWKJiIiIyAwWS0RERERmsFgiIiIiMoPFEhEREZEZLJaIiIiIzGCxRERERGQGiyUiIiIiM1gsEREREZnBYomIiIjIDBZLRERERGawWCIiIiIyg8USERERkRksloiIiIjMcLliafny5WjevDnUajXCwsKQmJhYaexTTz0FhUJR7vHcc89JMSNHjiz3elRUVE0sChEREbkAN0cnYI0tW7YgOjoaK1euRFhYGJYsWYLIyEicP38efn5+5eK3b9+OoqIi6fnNmzfRsWNH/O1vf5PFRUVFYd26ddJzT09P+y0EERERuRSXOrP0ySefYMyYMRg1ahSCg4OxcuVKeHt7Y+3atRXGN27cGDqdTnrs3bsX3t7e5YolT09PWVyjRo1qYnGIiIjIBbhMsVRUVIQTJ04gIiJCalMqlYiIiMCRI0csmseaNWswZMgQ1KtXT9Z+8OBB+Pn54eGHH8a4ceNw8+ZNs/MpLCxEdna27EFERES1k8sUSzdu3IDRaIS/v7+s3d/fH3q9vsrpExMTkZycjNGjR8vao6Ki8OWXXyIhIQELFy7EDz/8gN69e8NoNFY6rwULFkCr1UqPoKCg6i0UEVEdYDQJHLl4E/9KuoYjF2/CaBKOTonIKi7VZ+l+rFmzBh06dEBoaKisfciQIdL/O3TogEcffRStWrXCwYMH0bNnzwrnNWvWLERHR0vPs7OzWTAREVUgPjkNsbvOIM1QILUFaNWI6RuMqPYBDsyMyHIuc2bJx8cHKpUK6enpsvb09HTodDqz0+bl5WHz5s147bXXqvw7LVu2hI+PDy5cuFBpjKenJxo0aCB7EBGRXHxyGsZtPCkrlABAbyjAuI0nEZ+c5qDMiKzjMsWSh4cHunTpgoSEBKnNZDIhISEB4eHhZqfdtm0bCgsL8fLLL1f5d65evYqbN28iIIDfeIjIubjS5SyjSSB21xlUlGFZW+yuM069DERlXOoyXHR0NEaMGIGQkBCEhoZiyZIlyMvLw6hRowAAw4cPR9OmTbFgwQLZdGvWrMGAAQPQpEkTWXtubi5iY2MxcOBA6HQ6XLx4EdOnT0fr1q0RGRlpfYJ5eYBKVb5dpQLUanlcZZRKwMurerH5+YCo5MCjUMCo9kJiSiYycgqgczMhpFkjqJSKCmPh7X3n+e3bgMlUeR53d5i3JragADDTN8yqWG/v0rwBoLAQKCmxTayXV+l6BoCiIqC42DaxavWdfcWa2OLi0vjKeHoCbm7Wx5aUlK6Lynh4AO7u1scajaXbrjLu7qXx1saaTKX7mi1i3dxK1wVQ+v7Jz7dNrDXvewti9/6mx/x/n0VqdhEK3UtzCNCq8V7P5ni2XSVn1608Rsje99bEVvC+/8+lm8i6ngUvALc97iybZ3EhlP+bb9b1AvzntysIa3nXsZnHiFI8RpSqiWOEJYSLWbp0qXjwwQeFh4eHCA0NFUePHpVe69GjhxgxYoQs/ty5cwKA2LNnT7l55efni169eglfX1/h7u4umjVrJsaMGSP0er1VORkMBgFAGEoPLeUfffrIJ/D2rjgOEKJHD3msj0/lsSEh8thmzSqNzW71kHhi/j7RbEacaDYjTpxv8mDl823WTD7fkJDKY3185LE9elQe6+0tj+3Tp/LYe3fNQYPMx+bm3okdMcJ8bEbGndjx483HpqTciZ02zXxscvKd2JgY87GJiXdiFy0yH3vgwJ3YZcvMx8bF3Yldt8587Natd2K3bjUfu27dndi4OPOxy5bdiT1wwHzsokV3YhMTzcfGxNyJTU42Hztt2p3YlBTzsePH34nNyDAfe/fxJTfXfOygQULGXKwVx4gjQe2l93HzGXHihleDyudrxTFCBAfLY4ODK4+14hhxw6uBlG+zGXHiSFD7yufLY8SdB48RpQ87HyOkz2+DQZjjUmeWAGDixImYOHFiha8dPHiwXNvDDz8MIUSF8V5eXti9e7ct03NaaVkF5foNEJFrq/jIRkS2phCVVRJksezsbGi1WhhSUyvu7O3Ay3BGk0DEJweRll2IAvc7OaiLC6AQgAKAv9YT+6KfunNJjpfh7uAp9lKufIrdxS/DHbt0EyPXHZeemxQK6TIcAHgVla6z9aMel1/OAhx6Ga7s2JNuKER+BZfhKjz2ADxGlOExopSdjxHS57fBYPbHWi53Zsmp1asnf/Oai7Nmnpa6++D1P4kXbyLltgK4q1ACICucUm4DiRmFCG/V5N7JS919sK2KNbFqddUx1Yn19LzzgWbLWA8Py69x2yvW3f3OQcaWsW5udw6KtoxVqSzfh62JVSrtE6tQ2CcWqHas3pgl6/Nzr7LX9EYL1l8FxwibxFbwvlcBmDGwC8ZtPAkF7pwFK3T3RFlpNGNgZ6jqawCUFlelfSqz4FdfjdAWjaHiMcL6WB4jSlnzvrcAi6VaLiPHsktvlsYRUc3yq29ZEWBpXE2Kah+AFS93LjfOku6ecZY4FhM5OxZLtZwrH2iJCAht0RgBWjX0hoIK+ygpUFp8hLZoXNOpWSSqfQCeDdZJv8SVzhr979Jb2VhM9y5b2VhMK17uzIKJHM5lxlmi6ik70FYwQACA0gNtgBMfaInqOpVSgZi+wQBQ7n1c9jymb3DFw4A4CZVSgfBWTdC/U1OEt2oi5cqxmMhVsFiq5WrDgZaoriu7nKXTys8A67Rqlz7zkpiSafZXugJAmqEAiSmZNZcUUQV4Ga4OsLTfADmnOx1fy1/CoLqjqstZroh9KslVsFiqI2rjgbYuYMdXulvZ5azagn0qyVWwWKpDatuBtrZjx1eq7Vy98zrVHeyzROSE2PGV6gL2qSRXwWKJyAmx4yvVFbW18zrVLrwMR+SE2PGV6hL2qSRnx2KJyAmx4yvVNexTSc6Ml+GInBAHEyUich4sloicEDu+EhE5DxZLRE6KHV+JiJwD+ywROTF2fCUicjwWS0ROjh1fiYgci5fhiIiIiMxgsURERERkBoslIiIiIjNYLBERERGZwWKJiIiIyAwWS0RERERmcOgAIgcymgTHUCIicnIslogcJD45DbG7ziDNUCC1BWjViOkbzNG5iYicCC/DETlAfHIaxm08KSuUAEBvKMC4jScRn5zmoMyIiOheLJaIapjRJBC76wxEBa+VtcXuOgOjqaIIIiKqaSyWiGpYYkpmuTNKdxMA0gwFSEzJrLmkiIioUiyWiGpYRk7lhVJ14oiIyL5YLBHVML/6apvGERGRfblcsbR8+XI0b94carUaYWFhSExMrDR2/fr1UCgUsodaLf8AEkJg7ty5CAgIgJeXFyIiIvDf//7X3otBdVhoi8YI0KpR2QABCpT+Ki60ReOaTIuIiCrhUsXSli1bEB0djZiYGJw8eRIdO3ZEZGQkMjIyKp2mQYMGSEtLkx5//vmn7PVFixbhs88+w8qVK3Hs2DHUq1cPkZGRKCjgJRCyD5VSgZi+wQBQrmAqex7TN5jjLREROQmXKpY++eQTjBkzBqNGjUJwcDBWrlwJb29vrF27ttJpFAoFdDqd9PD395deE0JgyZIleOedd9C/f388+uij+PLLL5Gamopvv/22BpaI6qqo9gFY8XJn6LTyM506rRorXu7McZaIiJyIywxKWVRUhBMnTmDWrFlSm1KpREREBI4cOVLpdLm5uWjWrBlMJhM6d+6M+fPno127dgCAlJQU6PV6RERESPFarRZhYWE4cuQIhgwZUuE8CwsLUVhYKD3Pzs6+38WjOiiqfQCeDdZxBG8iIifnMmeWbty4AaPRKDszBAD+/v7Q6/UVTvPwww9j7dq1+Ne//oWNGzfCZDKha9euuHr1KgBI01kzTwBYsGABtFqt9AgKCrqfRaM6TKVUILxVE/Tv1BThrZqwUCIickIuUyxVR3h4OIYPH45OnTqhR48e2L59O3x9ffHPf/7zvuY7a9YsGAwG6XHlyhUbZUxERETOxmWKJR8fH6hUKqSnp8va09PTodPpLJqHu7s7HnvsMVy4cAEApOmsnaenpycaNGggexAREVHt5DLFkoeHB7p06YKEhASpzWQyISEhAeHh4RbNw2g04vTp0wgIKO0826JFC+h0Otk8s7OzcezYMYvnSURERLWby3TwBoDo6GiMGDECISEhCA0NxZIlS5CXl4dRo0YBAIYPH46mTZtiwYIFAID33nsPTzzxBFq3bo2srCx89NFH+PPPPzF69GgApb+Ue+utt/D++++jTZs2aNGiBebMmYPAwEAMGDDAUYtJRERETsSliqXBgwfj+vXrmDt3LvR6PTp16oT4+Hipg/bly5ehVN45WXbr1i2MGTMGer0ejRo1QpcuXXD48GEEBwdLMdOnT0deXh7Gjh2LrKwsPPnkk4iPjy83eCURERHVTQohBG9tfp+ys7Oh1WphMBjYf4mIiMhFWPr57TJ9loiIiIgcwapiqbi4GK1atcLZs2ftlQ/VYUaTwJGLN/GvpGs4cvEmjCae9CQiIsezqs+Su7s775lGdhGfnIbYXWeQZrizfwVo1YjpG8xbfxARkUNZfRluwoQJWLhwIUpKSuyRD9VB8clpGLfxpKxQAgC9oQDjNp5EfHKagzIjIiKqxq/hjh8/joSEBOzZswcdOnRAvXr1ZK9v377dZslR7Wc0CcTuOoOKLrgJAAoAsbvO4NlgHW8FQkREDmF1sdSwYUMMHDjQHrlQHZSYklnujNLdBIA0QwESUzIR3qpJzSVGRET0P1YXS+vWrbNHHlRHZeRY1gfO0jgiIiJbq9bQASUlJdi3bx/++c9/IicnBwCQmpqK3NxcmyZHtZ9ffcsG/7Q0joiIyNasPrP0559/IioqCpcvX0ZhYSGeffZZ1K9fHwsXLkRhYSFWrlxpjzyplgpt0RgBWjX0hoIK+y0pAOi0aoS2aFzTqREREQGoxpmlyZMnIyQkBLdu3YKXl5fU/sILL8huSEtkCZVSgZi+pbefubf7dtnzmL7B7NxNREQOY3Wx9NNPP+Gdd96Bh4eHrL158+a4du2azRKjuiOqfQBWvNwZOq38UptOq8aKlztznCUiInIoqy/DmUwmGI3Gcu1Xr15F/fr1bZIU1T1R7QPwbLAOiSmZyMgpgF/90ktvPKNERESOZvWZpV69emHJkiXSc4VCgdzcXMTExKBPnz62zI3qGJVSgfBWTdC/U1OEt2rCQomIiJyCQghh1Q24rl69isjISAgh8N///hchISH473//Cx8fH/z444/w8/OzV65Oy9K7FhMREZHzsPTz2+piCSgdOmDLli345ZdfkJubi86dO2PYsGGyDt91CYslIiIi12O3YunHH39E165d4eYm7+5UUlKCw4cPo3v37tXL2IWxWCIiInI9ln5+W91n6emnn0ZmZma5doPBgKefftra2RERERE5NauLJSEEFIryHW9v3rxZ7qa6RERERK7O4qED/vrXvwIo/fXbyJEj4enpKb1mNBrx66+/omvXrrbPkIhqFaNJcIgIInIpFhdLWq0WQOmZpfr168s6c3t4eOCJJ57AmDFjbJ8hEdUa8clpiN11BmmGOzdGDtCqEdM3mIOPElE5zvLlyuJiad26dQBKR+qeNm0aL7kRkVXik9MwbuPJcvcA1BsKMG7jSY7WTkQyzvTlyuo+SzExMfD09MS+ffvwz3/+Ezk5OQCA1NRU5Obm2jxBInJ9RpNA7K4zFd4suawtdtcZGE1Wj2RCRLVQ2Zeruwsl4M6Xq/jktBrNx+pi6c8//0SHDh3Qv39/TJgwAdevXwcALFy4ENOmTbN5gkTk+hJTMssd9O4mAKQZCpCYUv6XtkRUtzjjlyuri6XJkycjJCQEt27dkvVbeuGFF5CQkGDT5IiodsjIqbxQqk4cEdVezvjlyuob6f700084fPgwPDw8ZO3NmzfHtWvXbJYYEdUefvXVNo0jotrLGb9cWX1myWQywWg0lmu/evUq6tevb5OkiKh2CW3RGAFaNSr7DYsCpR03Q1s0rsm0iMgJOeOXK6uLpV69emHJkiXSc4VCgdzcXMTExKBPnz62zI2IagmVUoGYvsEAUK5gKnse0zeY4y0RkVN+ubK6WPrHP/6Bn3/+GcHBwSgoKMDQoUOlS3ALFy60R45EVAtEtQ/Aipc7Q6eVfxvUadUcNoCIJM745crqG+kCpTfN3bx5M3799Vfk5uaic+fOGDZsmKzDd13CG+kSWc5ZBpkjIudWE+MsWfr5Xa1iieRYLBEREdmevb9cWfr5bfWv4YDSASgPHTqEjIwMmEwm2WtvvvlmdWZJREREJKNSKhDeqomj07C+z9L69evRokULvPbaa/j444+xePFi6XF3x297Wb58OZo3bw61Wo2wsDAkJiZWGrt69Wr85S9/QaNGjdCoUSNERESUix85ciQUCoXsERUVZe/FICIiIhdhdbE0Z84czJ07FwaDAX/88QdSUlKkx6VLl+yRo2TLli2Ijo5GTEwMTp48iY4dOyIyMhIZGRkVxh88eBAvvfQSDhw4gCNHjiAoKAi9evUqNx5UVFQU0tLSpMfXX39t1+UgIiIi12F1n6UmTZogMTERrVq1sldOlQoLC8Pjjz+OZcuWASgd8ykoKAiTJk3CzJkzq5zeaDSiUaNGWLZsGYYPHw6g9MxSVlYWvv3222rnxT5LRERErsfSz2+rzyy99tpr2LZt230lVx1FRUU4ceIEIiIipDalUomIiAgcOXLEonnk5+ejuLgYjRvLx2Y4ePAg/Pz88PDDD2PcuHG4efOm2fkUFhYiOztb9iAiIqLayeoO3gsWLMDzzz+P+Ph4dOjQAe7u7rLXP/nkE5sld7cbN27AaDTC399f1u7v749z585ZNI8ZM2YgMDBQVnBFRUXhr3/9K1q0aIGLFy9i9uzZ6N27N44cOQKVSlXhfBYsWIDY2NjqLwwRERG5jGoVS7t378bDDz8MoHQE7zJ3/9/ZfPjhh9i8eTMOHjwItfrOoHhDhgyR/t+hQwc8+uijaNWqFQ4ePIiePXtWOK9Zs2YhOjpaep6dnY2goCD7JU9EREQOY3Wx9I9//ANr167FyJEj7ZBO5Xx8fKBSqZCeni5rT09Ph06nMzvtxx9/jA8//BD79u3Do48+aja2ZcuW8PHxwYULFyotljw9PeHp6WndAhAREZFLsrrPkqenJ7p162aPXMzy8PBAly5dkJCQILWZTCYkJCQgPDy80ukWLVqEefPmIT4+HiEhIVX+natXr+LmzZsICOCtF4iIiKgaxdLkyZOxdOlSe+RSpejoaKxevRobNmzA2bNnMW7cOOTl5WHUqFEAgOHDh2PWrFlS/MKFCzFnzhysXbsWzZs3h16vh16vR25uLgAgNzcXf//733H06FH88ccfSEhIQP/+/dG6dWtERkY6ZBmJiIjIuVh9GS4xMRH79+9HXFwc2rVrV66D9/bt222W3L0GDx6M69evY+7cudDr9ejUqRPi4+OlTt+XL1+GUnmn/luxYgWKioowaNAg2XxiYmLw7rvvQqVS4ddff8WGDRuQlZWFwMBA9OrVC/PmzeNlNiIiIgJQjXGWys7iVGbdunX3lZAr4jhLRERErsdu94ari8UQERER1V1W91kiIiIiqkusPrMEAN988w22bt2Ky5cvo6ioSPbayZMnbZIYkbMzmgQSUzKRkVMAv/pqhLZoDJXSeccaIyKi6rH6zNJnn32GUaNGwd/fH6dOnUJoaCiaNGmCS5cuoXfv3vbIkcjpxCen4cmF+/HS6qOYvDkJL60+iicX7kd8cpqjUyMiIhuzulj6/PPPsWrVKixduhQeHh6YPn069u7dizfffBMGg8EeORI5lfjkNIzbeBJphgJZu95QgHEbT7JgIiKqZawuli5fvoyuXbsCALy8vJCTkwMAeOWVV/D111/bNjsiJ2M0CcTuOoOKfkJa1ha76wyMJqt+ZEpERE7M6mJJp9MhMzMTAPDggw/i6NGjAICUlBRYOQoBkctJTMksd0bpbgJAmqEAiSmZNZcUERHZldXF0jPPPIOdO3cCKB1zacqUKXj22WcxePBgvPDCCzZPkMiZZORUXihVJ46IiJyf1b+GW7VqFUwmEwBgwoQJaNKkCQ4fPox+/frh9ddft3mCRM7Er77apnFEROT8rCqWSkpKMH/+fLz66qt44IEHAABDhgzBkCFD7JIckbMJbdEYAVo19IaCCvstKQDotKXDCBARUe1g1WU4Nzc3LFq0CCUlJfbKh8ipqZQKxPQNBlBaGN2t7HlM32COt0REVItY3WepZ8+e+OGHH+yRC5FLiGofgBUvd4ZOK7/UptOqseLlzohqH+CgzIiIyB6s7rPUu3dvzJw5E6dPn0aXLl1Qr1492ev9+vWzWXJEziqqfQCeDdZxBG8iojpAIaz8vb9SWfnJKIVCAaPReN9JuRpL71pMREREzsPSz2+rzyyV/RKOiIiIqC6wus8SERERUV1i9ZklAMjLy8MPP/yAy5cvo6ioSPbam2++aZPEiIiIiJyB1cXSqVOn0KdPH+Tn5yMvLw+NGzfGjRs34O3tDT8/PxZLREREVKtYfRluypQp6Nu3L27dugUvLy8cPXoUf/75J7p06YKPP/7YHjkSEREROYzVxVJSUhKmTp0KpVIJlUqFwsJCBAUFYdGiRZg9e7Y9ciQiIiJyGKuLJXd3d2n4AD8/P1y+fBkAoNVqceXKFdtmR0RERORgVvdZeuyxx3D8+HG0adMGPXr0wNy5c3Hjxg383//9H9q3b2+PHImIiIgcxuozS/Pnz0dAQOntHD744AM0atQI48aNw/Xr17Fq1SqbJ0hERETkSFaP4E3lcQRvIiIi12O3EbzLZGRk4Pz58wCAtm3bwtfXt7qzIqq1jCbB+8cREbk4q4ulnJwcjB8/Hps3b5buA6dSqTB48GAsX74cWq3W5kkSuaL45DTE7jqDNEOB1BagVSOmbzCi2gc4MDMiIrKG1X2WRo8ejWPHjiEuLg5ZWVnIyspCXFwc/vOf/+D111+3R45ELic+OQ3jNp6UFUoAoDcUYNzGk4hPTnNQZkREZC2r+yzVq1cPu3fvxpNPPilr/+mnnxAVFYW8vDybJugK2GeJ7mY0CTy5cH+5QqmMAoBOq8ahGc/wkhwRkQNZ+vlt9ZmlJk2aVHipTavVolGjRtbOjqjWSUzJrLRQAgABIM1QgMSUzJpLioiIqs3qYumdd95BdHQ09Hq91KbX6/H3v/8dc+bMsWlyRK4oI6fyQqk6cURE5FhWd/BesWIFLly4gAcffBAPPvggAODy5cvw9PTE9evX8c9//lOKPXnypO0yJXIRfvXVNo0jIiLHsrpYGjBggB3SsNzy5cvx0UcfQa/Xo2PHjli6dClCQ0Mrjd+2bRvmzJmDP/74A23atMHChQvRp08f6XUhBGJiYrB69WpkZWWhW7duWLFiBdq0aVMTi0O1UGiLxgjQqqE3FKCiDoFlfZZCWzSu6dSIiKgaXGpQyi1btmD48OFYuXIlwsLCsGTJEmzbtg3nz5+Hn59fufjDhw+je/fuWLBgAZ5//nl89dVXWLhwIU6ePCndmmXhwoVYsGABNmzYgBYtWmDOnDk4ffo0zpw5A7Xasm/+7OBN9yr7NRwAWcFU1p17xcudOXwAEZGDWfr57VLFUlhYGB5//HEsW7YMAGAymRAUFIRJkyZh5syZ5eIHDx6MvLw8xMXFSW1PPPEEOnXqhJUrV0IIgcDAQEydOhXTpk0DABgMBvj7+2P9+vUYMmSIRXmVrezU1FTodDooFKUfiUVFRSguLoabmxs8PT2l+LJfDHp5eUk3JS4uLkZRURFUKpWsSLMmNj8/H0IIqNVqqFQqAEBJSQkKCwuhVCrh5eVVrdjbt2/DZDLB09MTbm6lJyONRiMKCgqsilUoFPD29pZiCwoKYDQa4eHhAXd3d6tjTSYTbt++DaD0V5plCgsLUVJSAnd3d3h4eFgdK4RAfn4+AMDb27vc9rQ09vukK1iw+3ek55mkv+fnJTC79yPoF9JC2p622E8q2p622E/Ktuf97if3bs/73U8q2573u5/cvT3vdz+pbHtaE8tjRO0+Rli67XmMsN8xwuKTHcJFFBYWCpVKJXbs2CFrHz58uOjXr1+F0wQFBYnFixfL2ubOnSseffRRIYQQFy9eFADEqVOnZDHdu3cXb775ZqW5FBQUCIPBID2uXLkiUHoCQWRkZEhx77//vgAgRo8eLZve29tbABApKSlS2+LFiwUAMXToUFmsj4+PACCSk5OltlWrVgkAon///rLYZs2aCQAiMTFRatu4caMAICIiImSxwcHBAoA4cOCA1LZjxw4BQHTt2lUWGxISIgCIuLg4qW3Pnj0CgOjYsaMstkePHgKA2Lp1q9R26NAhAUC0bt1aFtunTx8BQKxbt05qO3XqlAAgAgMDZbGDBg0SAMSyZcuktt9//10AEFqtVhY7YsQIAUAsWrRIart69aoAINzc3GSx48ePFwBETEyM1Hbr1i1pexYVFUnt06ZNEwDEtGnTpLaioiIp9tatW1J7TEyMACDGjRsnDl+4Ib49dVUcvnBDuLm5CQDi6tWrUuyiRYsEADFixAhZblqtVgAQv//+u9S2bNkyAUAMGjRIFhsYGFhuX163bp0AIPr06SOLbd26tQAgDh06JLVt3bpVABA9evSQxXbs2FEAEHv27JHa4uLiBAAREhIii+3atasAIHuPHjhwQAAQwcHBstiIiAgBQGzcuFFqS0xMFABEs2bNZLH9+/cXAMSqVauktuTkZAFA+Pj4yGKHDh0qAMje9ykpKQKA8Pb2lsWOHj1aABDvv/++1JaRkSFtz7tNnjxZABCzZ8+W2nJzc6XY3NxcqX327NkCgJg8ebJsHjxGlOIxolTZMWL8+PGyv8djRKmaOkYYDAYBQBgMBmFOtW93UtNu3LgBo9EIf39/Wbu/vz/OnTtX4TR6vb7C+LJf8pX9ay6mIgsWLEBsbKzVy0B1j0KhQHirJo5Og4iI7oPLXIZLTU1F06ZNcfjwYYSHh0vt06dPxw8//IBjx46Vm8bDwwMbNmzASy+9JLV9/vnniI2NRXp6Og4fPoxu3bohNTUVAQF3+o+8+OKLUCgU2LJlS4W5FBYWorCwUHqenZ2NoKAgXoarIpan2HmKvTqxvAzHYwSPETxG2GrbV/cynMucWfLx8YFKpUJ6erqsPT09HTqdrsJpdDqd2fiyf9PT02XFUnp6Ojp16lRpLp6enrIdtky9evWkNwJQWqyVvWHujbuXu7u7tJGrG3v3zlPGzc1N2tGqG3v3zl5GpVJVmJs1sRV1oLcmVqlUVhhb0faxJlahUFQYW9H2tCYWqHh72mI/qWh72mI/qWh7OsN+Utn2vN/9pLLteb/7CWC/bc9jROWxPEbcwWNE5bGWsHpQSqPRiDVr1mDo0KGIiIjAM888I3vYi4eHB7p06YKEhASpzWQyISEhQXam6W7h4eGyeADYu3evFN+iRQvodDpZTHZ2No4dO1bpPImIiKhusfrM0uTJk7F+/Xo899xzaN++vexMir1FR0djxIgRCAkJQWhoKJYsWYK8vDyMGjUKADB8+HA0bdoUCxYskHLt0aMH/vGPf+C5557D5s2b8Z///AerVq0CUFrxv/XWW3j//ffRpk0baeiAwMBAh48nRURERM7B6mJp8+bN2Lp1q2xgx5oyePBgXL9+HXPnzoVer0enTp0QHx8vddC+fPmydN0WALp27YqvvvoK77zzDmbPno02bdrg22+/lcZYAkr7POXl5WHs2LHIysrCk08+ifj4+GqfqiMiIqLaxeoO3oGBgTh48CAeeughe+XkcjgoJRERkeux9PPb6j5LU6dOxaeffgoX+REdERER0X2x+jLcoUOHcODAAfz73/9Gu3btyvWk3759u82SIyIiInI0q4ulhg0b4oUXXrBHLkREREROx+piad26dfbIg4iIiMgpVXtQyuvXr+P8+fMAgIcffhi+vr42S4qIiIjIWVjdwTsvLw+vvvoqAgIC0L17d3Tv3h2BgYF47bXXpCHdiYiIiGoLq4ul6Oho/PDDD9i1axeysrKQlZWFf/3rX/jhhx8wdepUe+RIRERE5DBWj7Pk4+ODb775Bk899ZSs/cCBA3jxxRdx/fp1W+bnEjjOEhERkeux2zhL+fn50ojZd/Pz8+NlOCIiIqp1rC6WwsPDERMTg4KCAqnt9u3biI2N5c1niYiIqNax+tdwn376KSIjI/HAAw+gY8eOAIBffvkFarUau3fvtnmCRERERI5kdZ8loPRS3KZNm3Du3DkAwCOPPIJhw4bBy8vL5gm6AvZZIiIicj2Wfn5Xa5wlb29vjBkzptrJEREREbkKi4qlnTt3onfv3nB3d8fOnTvNxvbr188miRERERE5A4suwymVSuj1evj5+UGprLxPuEKhgNFotGmCroCX4YiIiFyPTS/DmUymCv9PREREVNtZPXTAl19+icLCwnLtRUVF+PLLL22SFBEREZGzsPrXcCqVCmlpafDz85O137x5E35+frwMx8twRERELsFuI3gLIaBQKMq1X716FVqt1trZERERETk1i4cOeOyxx6BQKKBQKNCzZ0+4ud2Z1Gg0IiUlBVFRUXZJkoiIiMhRLC6WBgwYAABISkpCZGQkNBqN9JqHhweaN2+OgQMH2jxBIiIiIkeyuFiKiYkBADRv3hyDBw+GWq22W1JEREREzsLqEbxHjBhhjzyIiIiInJLVxZLRaMTixYuxdetWXL58GUVFRbLXMzMzbZYcERERkaNZ/Wu42NhYfPLJJxg8eDAMBgOio6Px17/+FUqlEu+++64dUiQiIiJyHKuLpU2bNmH16tWYOnUq3Nzc8NJLL+GLL77A3LlzcfToUXvkSEREROQwVhdLer0eHTp0AABoNBoYDAYAwPPPP4/vvvvOttkREREROZjVxdIDDzyAtLQ0AECrVq2wZ88eAMDx48fh6elp2+yIiIiIHMzqYumFF15AQkICAGDSpEmYM2cO2rRpg+HDh+PVV1+1eYJEREREjmT1veHudeTIERw5cgRt2rRB3759bZWXS+G94YiIiFyPpZ/fVg8dcK/w8HCEh4ff72yIiIiInJJFxdLOnTstnmG/fv2qnYw5mZmZmDRpEnbt2gWlUomBAwfi008/ld125d74mJgY7NmzB5cvX4avry8GDBiAefPmyW74W9FNgb/++msMGTLELstBRERErsWiYqnsvnBVUSgUMBqN95NPpYYNG4a0tDTs3bsXxcXFGDVqFMaOHYuvvvqqwvjU1FSkpqbi448/RnBwMP7880+88cYbSE1NxTfffCOLXbdunewmwA0bNrTLMhAREZHrue8+SzXh7NmzCA4OxvHjxxESEgIAiI+PR58+fXD16lUEBgZaNJ9t27bh5ZdfRl5eHtzcSutEhUKBHTt2WFwQVoR9loiIiFyPpZ/fVv8a7m4FBQX3M7nFjhw5goYNG0qFEgBERERAqVTi2LFjFs+nbGWUFUplJkyYAB8fH4SGhmLt2rWoqn4sLCxEdna27EFERES1k9XFktFoxLx589C0aVNoNBpcunQJADBnzhysWbPG5gkCpQNh+vn5ydrc3NzQuHFj6PV6i+Zx48YNzJs3D2PHjpW1v/fee9i6dSv27t2LgQMHYvz48Vi6dKnZeS1YsABarVZ6BAUFWbdARERE5DKsLpY++OADrF+/HosWLYKHh4fU3r59e3zxxRdWzWvmzJlQKBRmH+fOnbM2xXKys7Px3HPPITg4uNz96+bMmYNu3brhsccew4wZMzB9+nR89NFHZuc3a9YsGAwG6XHlypX7zpGIiIick9VDB3z55ZdYtWoVevbsiTfeeENq79ixo9WFzdSpUzFy5EizMS1btoROp0NGRoasvaSkBJmZmdDpdGanz8nJQVRUFOrXr48dO3bA3d3dbHxYWBjmzZuHwsLCSkck9/T05GjlREREdYTVxdK1a9fQunXrcu0mkwnFxcVWzcvX1xe+vr5VxoWHhyMrKwsnTpxAly5dAAD79++HyWRCWFhYpdNlZ2cjMjISnp6e2LlzJ9RqdZV/KykpCY0aNWIxRERERACqcRkuODgYP/30U7n2b775Bo899phNkrrXI488gqioKIwZMwaJiYn4+eefMXHiRAwZMkT6Jdy1a9fQtm1bJCYmAigtlHr16oW8vDysWbMG2dnZ0Ov10Ov10vAGu3btwhdffIHk5GRcuHABK1aswPz58zFp0iS7LAcRERG5HqvPLM2dOxcjRozAtWvXYDKZsH37dpw/fx5ffvkl4uLi7JEjAGDTpk2YOHEievbsKQ1K+dlnn0mvFxcX4/z588jPzwcAnDx5Uvql3L1nwlJSUtC8eXO4u7tj+fLlmDJlCoQQaN26NT755BOMGTPGbstBRERErqVa4yz99NNPeO+99/DLL78gNzcXnTt3xty5c9GrVy975Oj0OM4SERGR67HLveFKSkowf/58vPrqq9i7d+99J0lERETk7Kzqs+Tm5oZFixahpKTEXvkQERERORWrO3j37NkTP/zwgz1yISIiInI6Vnfw7t27N2bOnInTp0+jS5cuqFevnuz1fv362Sw5IiIiIkezuoO3Uln5ySiFQiH9LL8uYQdvIiIi12OXDt5A6eCTRERERHWFVX2WiouL4ebmhuTkZHvlQ0RERORUrCqW3N3d8eCDD9bJS21ERERUN1n9a7i3334bs2fPRmZmpj3yISIiInIqVvdZWrZsGS5cuIDAwEA0a9as3K/hTp48abPkiIiIiBzN6mJpwIABdkiDiIiIyDlV695wJMehA4iIiFyP3YYOKHPixAmcPXsWANCuXTs89thj1Z0VERERkdOyuljKyMjAkCFDcPDgQTRs2BAAkJWVhaeffhqbN2+Gr6+vrXMkIiIichirfw03adIk5OTk4LfffkNmZiYyMzORnJyM7OxsvPnmm/bIkYiIiMhhrO6zpNVqsW/fPjz++OOy9sTERPTq1QtZWVm2zM8lsM8SERGR67H089vqM0smkwnu7u7l2t3d3XkrFCIiIqp1rC6WnnnmGUyePBmpqalS27Vr1zBlyhT07NnTpskREREROZrVxdKyZcuQnZ2N5s2bo1WrVmjVqhVatGiB7OxsLF261B45EhERETmM1b+GCwoKwsmTJ7Fv3z6cO3cOAPDII48gIiLC5skRERERORoHpbQBdvAmIiJyPTbv4L1//34EBwcjOzu73GsGgwHt2rXDTz/9VL1siYiIiJyUxcXSkiVLMGbMmAorL61Wi9dffx2ffPKJTZMjIiIicjSLi6VffvkFUVFRlb7eq1cvnDhxwiZJERERETkLi4ul9PT0CsdXKuPm5obr16/bJCkiIiIiZ2FxsdS0aVMkJydX+vqvv/6KgIAAmyRFRERE5CwsLpb69OmDOXPmoKCgoNxrt2/fRkxMDJ5//nmbJkdERETkaBYPHZCeno7OnTtDpVJh4sSJePjhhwEA586dw/Lly2E0GnHy5En4+/vbNWFnxKEDiIiIXI+ln98WD0rp7++Pw4cPY9y4cZg1axbKaiyFQoHIyEgsX768ThZKREREVLtZNYJ3s2bN8P333+PWrVu4cOEChBBo06YNGjVqZK/8iIiIiBzK6tudAECjRo3w+OOP2zoXIiIiIqdj9Y10HSUzMxPDhg1DgwYN0LBhQ7z22mvIzc01O81TTz0FhUIhe7zxxhuymMuXL+O5556Dt7c3/Pz88Pe//x0lJSX2XBQiIiJyIdU6s+QIw4YNQ1paGvbu3Yvi4mKMGjUKY8eOxVdffWV2ujFjxuC9996Tnnt7e0v/NxqNeO6556DT6XD48GGkpaVh+PDhcHd3x/z58+22LEREROQ6XOJGumfPnkVwcDCOHz+OkJAQAEB8fDz69OmDq1evIjAwsMLpnnrqKXTq1AlLliyp8PV///vfeP7555Gamip1Tl+5ciVmzJiB69evw8PDw6L8+Gs4IiIi12PzG+k60pEjR9CwYUOpUAKAiIgIKJVKHDt2zOy0mzZtgo+PD9q3b49Zs2YhPz9fNt8OHTrIfsUXGRmJ7Oxs/Pbbb5XOs7CwENnZ2bIHERER1U4ucRlOr9fDz89P1ubm5obGjRtDr9dXOt3QoUPRrFkzBAYG4tdff8WMGTNw/vx5bN++XZrvvcMdlD03N98FCxYgNja2uotDRERELsShxdLMmTOxcOFCszFnz56t9vzHjh0r/b9Dhw4ICAhAz549cfHiRbRq1ara8501axaio6Ol59nZ2QgKCqr2/IiIiMh5ObRYmjp1KkaOHGk2pmXLltDpdMjIyJC1l5SUIDMzEzqdzuK/FxYWBgC4cOECWrVqBZ1Oh8TERFlMeno6AJidr6enJzw9PS3+u0REROS6HFos+fr6wtfXt8q48PBwZGVl4cSJE+jSpQsAYP/+/TCZTFIBZImkpCQAkG74Gx4ejg8++AAZGRnSZb69e/eiQYMGCA4OtnJpiIiIqDZyiQ7ejzzyCKKiojBmzBgkJibi559/xsSJEzFkyBDpl3DXrl1D27ZtpTNFFy9exLx583DixAn88ccf2LlzJ4YPH47u3bvj0UcfBQD06tULwcHBeOWVV/DLL79g9+7deOeddzBhwgSeOSIiIiIALlIsAaW/amvbti169uyJPn364Mknn8SqVauk14uLi3H+/Hnp124eHh7Yt28fevXqhbZt22Lq1KkYOHAgdu3aJU2jUqkQFxcHlUqF8PBwvPzyyxg+fLhsXCYiIiKq21xinCVnx3GWiIiIXE+tGmeJiIiIyFFYLBERERGZwWKJiIiIyAwWS0RERERmsFgiIiIiMoPFEhEREZEZLJaIiIiIzGCxRERERGQGiyUiIiIiM1gsEREREZnBYomIiIjIDBZLRERERGawWCIiIiIyg8USERERkRksloiIiIjMYLFEREREZAaLJSIiIiIzWCwRERERmcFiiYiIiMgMFktEREREZrBYIiIiIjKDxRIRERGRGSyWiIiIiMxgsURERERkBoslIiIiIjNYLBERERGZwWKJiIiIyAwWS0RERERmsFgiIiIiMoPFEhEREZEZLJaIiIiIzHCZYikzMxPDhg1DgwYN0LBhQ7z22mvIzc2tNP6PP/6AQqGo8LFt2zYprqLXN2/eXBOLRERERC7AzdEJWGrYsGFIS0vD3r17UVxcjFGjRmHs2LH46quvKowPCgpCWlqarG3VqlX46KOP0Lt3b1n7unXrEBUVJT1v2LChzfMnIiIi1+QSxdLZs2cRHx+P48ePIyQkBACwdOlS9OnTBx9//DECAwPLTaNSqaDT6WRtO3bswIsvvgiNRiNrb9iwYblYIiIiIsBFLsMdOXIEDRs2lAolAIiIiIBSqcSxY8csmseJEyeQlJSE1157rdxrEyZMgI+PD0JDQ7F27VoIIczOq7CwENnZ2bIHERER1U4ucWZJr9fDz89P1ubm5obGjRtDr9dbNI81a9bgkUceQdeuXWXt7733Hp555hl4e3tjz549GD9+PHJzc/Hmm29WOq8FCxYgNjbW+gUhIiIil+PQM0szZ86stBN22ePcuXP3/Xdu376Nr776qsKzSnPmzEG3bt3w2GOPYcaMGZg+fTo++ugjs/ObNWsWDAaD9Lhy5cp950hERETOyaFnlqZOnYqRI0eajWnZsiV0Oh0yMjJk7SUlJcjMzLSor9E333yD/Px8DB8+vMrYsLAwzJs3D4WFhfD09KwwxtPTs9LXiIiIqHZxaLHk6+sLX1/fKuPCw8ORlZWFEydOoEuXLgCA/fv3w2QyISwsrMrp16xZg379+ln0t5KSktCoUSMWQ0RERATARfosPfLII4iKisKYMWOwcuVKFBcXY+LEiRgyZIj0S7hr166hZ8+e+PLLLxEaGipNe+HCBfz444/4/vvvy813165dSE9PxxNPPAG1Wo29e/di/vz5mDZtWo0tGxERETk3lyiWAGDTpk2YOHEievbsCaVSiYEDB+Kzzz6TXi8uLsb58+eRn58vm27t2rV44IEH0KtXr3LzdHd3x/LlyzFlyhQIIdC6dWt88sknGDNmjN2Xh4iIiFyDQlT1O3mqUnZ2NrRaLQwGAxo0aODodIiIiMgCln5+u8Q4S0RERESOwmKJiIiIyAwWS0RERERmsFgiIiIiMoPFEhEREZEZLJaIiIiIzGCxRERERGQGiyUiIiIiM1gsEREREZnBYomIiIjIDBZLRERERGawWCIiIiIyg8USERERkRksloiIiIjMYLFEREREZAaLJSIiIiIzWCwRERERmcFiiYiIiMgMFktEREREZrBYIiIiIjKDxRIRERGRGSyWiIiIiMxgsURERERkBoslIiIiIjNYLBERERGZwWKJiIiIyAwWS0RERERmsFgiIiIiMoPFEhEREZEZLJaIiIiIzHBzdAJkHaNJIDElExk5BfCrr0Zoi8ZQKRVWxxAREZFlXKZY+uCDD/Ddd98hKSkJHh4eyMrKqnIaIQRiYmKwevVqZGVloVu3blixYgXatGkjxWRmZmLSpEnYtWsXlEolBg4ciE8//RQajcaOS1M98clpiN11BmmGAqktQKtGTN9gRLUPsDiGiIiILOcyl+GKiorwt7/9DePGjbN4mkWLFuGzzz7DypUrcezYMdSrVw+RkZEoKLhTSAwbNgy//fYb9u7di7i4OPz4448YO3asPRbhvsQnp2HcxpOyIggA9IYCjNt4EvHJaRbFEBERkXUUQgjh6CSssX79erz11ltVnlkSQiAwMBBTp07FtGnTAAAGgwH+/v5Yv349hgwZgrNnzyI4OBjHjx9HSEgIACA+Ph59+vTB1atXERgYaFFO2dnZ0Gq1SE1NhU6ng0JResmrqKgIxcXFcHNzg6enpxSfl5cHAPDy8oJSWVqvFhcXo6ioCCqVCmq1WhZrNAn0+uwI9DlFpctmLIEwlkChVELh5gEFAJ1WjZLC20jPLoTCzR0Kpao01mSEKCmGQqFAoI8Wh2Y8A5VSgfz8fAghoFaroVKVxpaUlKCwsBBKpRJeXl5SDrdv34bJZIKnpyfc3EpPRhqNRhQUFFgVq1Ao4O3tLcUWFBTAaDTCw8MD7u7uVseaTCbcvn0bAFCvXj0ptrCwECUlJXB3d4eHh4fVsUII5OfnAwC8vb3LbU9rYi3Z9rbYTyrantbEVrXt73c/uXd73u9+Utn2vN/95O7teb/7SWXbs7r7ibljhKWx1mx7HiMqj+UxovYcI8o+vw0GAxo0aIBKCRezbt06odVqq4y7ePGiACBOnTola+/evbt48803hRBCrFmzRjRs2FD2enFxsVCpVGL79u2VzrugoEAYDAbpceXKFQFAABAZGRlS3Pvvvy8AiNGjR8um9/b2FgBESkqK1LZ48WIBQAwdOlQW6+PjIwCIgFeXi2Yz4kSzGXGiceREAUB4tXlCams2I06oGvgJAEI3/BOprcnzUwUAoW7WSTSbEScOX7ghhBAiODhYABAHDhyQ/taOHTsEANG1a1dZDiEhIQKAiIuLk9r27NkjAIiOHTvKYnv06CEAiK1bt0pthw4dEgBE69atZbF9+vQRAMS6deuktlOnTgkAIjAwUBY7aNAgAUAsW7ZMavv9998FgHL7w4gRIwQAsWjRIqnt6tWrAoBwc3OTxY4fP14AEDExMVLbrVu3pO1ZVFQktU+bNk0AENOmTZPaioqKpNhbt25J7TExMQKAGD9+vOzvubm5CQDi6tWrUtuiRYsEADFixAhZrFarFQDE77//LrUtW7ZMABCDBg2SxQYGBpbb39etWycAiD59+shiW7duLQCIQ4cOSW1bt24VAESPHj1ksR07dhQAxJ49e6S2uLg4AUCEhITIYrt27SoAiB07dkhtBw4cEABEcHCwLDYiIkIAEBs3bpTaEhMTBQDRrFkzWWz//v0FALFq1SqpLTk5WQAQPj4+stihQ4cKAGLx4sVSW0pKigAgvL29ZbGjR48WAMT7778vtWVkZEjb826TJ08WAMTs2bOlttzcXCk2NzdXap89e7YAICZPniybh72PEcnJyVLbqlWrBADRv39/WWyzZs0EAJGYmCi1bdy4UQAQERERslgeI0rxGFGqth4jDAaDACAMBoMwx2Uuw1lLr9cDAPz9/WXt/v7+0mt6vR5+fn6y193c3NC4cWMppiILFiyAVquVHkFBQTbO3n4ycgqqDiIiIiKJQy/DzZw5EwsXLjQbc/bsWbRt21Z6bulluMOHD6Nbt25ITU1FQMCdjs0vvvgiFAoFtmzZgvnz52PDhg04f/68bFo/Pz/ExsZW2j+qsLAQhYWF0vPs7GwEBQXZ7TLcsUs3MWrjL1AoSmPvvQxXxlRcAAhUeBkOCgWU7p74eswTCG/VhKfYzcQKnmIvtz1d9RT7/Wx7XobjMaKyWB4jas8xwtLLcA4tlq5fv46bN2+ajWnZsqW00wGWF0uXLl1Cq1atcOrUKXTq1Elq79GjBzp16oRPP/0Ua9euxdSpU3Hr1i3p9ZKSEqjVamzbtg0vvPCCRcth8TXPajKaBJ5cuB96QwEq2lhlfZaEEEjPLjQbU9ZniYiIqK6z9PPboUMH+Pr6wtfX1y7zbtGiBXQ6HRISEqRiKTs7G8eOHZPOGIWHhyMrKwsnTpxAly5dAAD79++HyWRCWFiYXfKqDpVSgZi+wRi38SQUgKwYKit7YvoGA0CVMSyUiIiIrOMyfZYuX76MpKQkXL58GUajEUlJSUhKSkJubq4U07ZtW+zYsQMAoFAo8NZbb+H999/Hzp07cfr0aQwfPhyBgYEYMGAAAOCRRx5BVFQUxowZg8TERPz888+YOHEihgwZYvEv4WpKVPsArHi5M3Rataxdp1VjxcudEdU+wKIYIiIiso7LDEo5d+5cbNiwQXr+2GOPAQAOHDiAp556CgBw/vx5GAwGKWb69OnIy8vD2LFjkZWVhSeffBLx8fGya7ObNm3CxIkT0bNnT2lQys8++6xmFspKUe0D8Gywzuzo3JbEEBERkeVcbpwlZ2TvPktERERke5Z+frvMZTgiIiIiR2CxRERERGQGiyUiIiIiM1gsEREREZnBYomIiIjIDBZLRERERGawWCIiIiIyg8USERERkRksloiIiIjMcJnbnTizskHQs7OzHZwJERERWarsc7uqm5mwWLKBnJwcAEBQUJCDMyEiIiJr5eTkQKvVVvo67w1nAyaTCampqahfvz4UiopvWJudnY2goCBcuXKF94+rANePeVw/VeM6Mo/rp2pcR+bVxvUjhEBOTg4CAwOhVFbeM4lnlmxAqVTigQcesCi2QYMGtWYnsweuH/O4fqrGdWQe10/VuI7Mq23rx9wZpTLs4E1ERERkBoslIiIiIjNYLNUQT09PxMTEwNPT09GpOCWuH/O4fqrGdWQe10/VuI7Mq8vrhx28iYiIiMzgmSUiIiIiM1gsEREREZnBYomIiIjIDBZLRERERGawWLKh5cuXo3nz5lCr1QgLC0NiYqLZ+G3btqFt27ZQq9Xo0KEDvv/++xrK1DGsWT+rV6/GX/7yFzRq1AiNGjVCRERElevT1Vm7/5TZvHkzFAoFBgwYYN8EnYC16ygrKwsTJkxAQEAAPD098dBDD9Xq95m162fJkiV4+OGH4eXlhaCgIEyZMgUFBQU1lG3N+vHHH9G3b18EBgZCoVDg22+/rXKagwcPonPnzvD09ETr1q2xfv16u+fpKNaun+3bt+PZZ5+Fr68vGjRogPDwcOzevbtmknUEQTaxefNm4eHhIdauXSt+++03MWbMGNGwYUORnp5eYfzPP/8sVCqVWLRokThz5ox45513hLu7uzh9+nQNZ14zrF0/Q4cOFcuXLxenTp0SZ8+eFSNHjhRarVZcvXq1hjOvGdaunzIpKSmiadOm4i9/+Yvo379/zSTrINauo8LCQhESEiL69OkjDh06JFJSUsTBgwdFUlJSDWdeM6xdP5s2bRKenp5i06ZNIiUlRezevVsEBASIKVOm1HDmNeP7778Xb7/9tti+fbsAIHbs2GE2/tKlS8Lb21tER0eLM2fOiKVLlwqVSiXi4+NrJuEaZu36mTx5sli4cKFITEwUv//+u5g1a5Zwd3cXJ0+erJmEaxiLJRsJDQ0VEyZMkJ4bjUYRGBgoFixYUGH8iy++KJ577jlZW1hYmHj99dftmqejWLt+7lVSUiLq168vNmzYYK8UHao666ekpER07dpVfPHFF2LEiBG1vliydh2tWLFCtGzZUhQVFdVUig5l7fqZMGGCeOaZZ2Rt0dHRolu3bnbN0xlYUgxMnz5dtGvXTtY2ePBgERkZacfMnIMl66ciwcHBIjY21vYJOQFehrOBoqIinDhxAhEREVKbUqlEREQEjhw5UuE0R44ckcUDQGRkZKXxrqw66+de+fn5KC4uRuPGje2VpsNUd/2899578PPzw2uvvVYTaTpUddbRzp07ER4ejgkTJsDf3x/t27fH/PnzYTQaayrtGlOd9dO1a1ecOHFCulR36dIlfP/99+jTp0+N5Ozs6tIx2hZMJhNycnJq5TEa4I10beLGjRswGo3w9/eXtfv7++PcuXMVTqPX6yuM1+v1dsvTUaqzfu41Y8YMBAYGljt41QbVWT+HDh3CmjVrkJSUVAMZOl511tGlS5ewf/9+DBs2DN9//z0uXLiA8ePHo7i4GDExMTWRdo2pzvoZOnQobty4gSeffBJCCJSUlOCNN97A7NmzayJlp1fZMTo7Oxu3b9+Gl5eXgzJzTh9//DFyc3Px4osvOjoVu+CZJXJ6H374ITZv3owdO3ZArVY7Oh2Hy8nJwSuvvILVq1fDx8fH0ek4LZPJBD8/P6xatQpdunTB4MGD8fbbb2PlypWOTs0pHDx4EPPnz8fnn3+OkydPYvv27fjuu+8wb948R6dGLuarr75CbGwstm7dCj8/P0enYxc8s2QDPj4+UKlUSE9Pl7Wnp6dDp9NVOI1Op7Mq3pVVZ/2U+fjjj/Hhhx9i3759ePTRR+2ZpsNYu34uXryIP/74A3379pXaTCYTAMDNzQ3nz59Hq1at7Jt0DavOPhQQEAB3d3eoVCqp7ZFHHoFer0dRURE8PDzsmnNNqs76mTNnDl555RWMHj0aANChQwfk5eVh7NixePvtt6FU1u3v0pUdoxs0aMCzSnfZvHkzRo8ejW3bttXKM/9l6va7wUY8PDzQpUsXJCQkSG0mkwkJCQkIDw+vcJrw8HBZPADs3bu30nhXVp31AwCLFi3CvHnzEB8fj5CQkJpI1SGsXT9t27bF6dOnkZSUJD369euHp59+GklJSQgKCqrJ9GtEdfahbt264cKFC1IhCQC///47AgICalWhBFRv/eTn55criMoKS8FbhtapY3R1ff311xg1ahS+/vprPPfcc45Ox74c3cO8tti8ebPw9PQU69evF2fOnBFjx44VDRs2FHq9XgghxCuvvCJmzpwpxf/888/Czc1NfPzxx+Ls2bMiJiam1g8dYM36+fDDD4WHh4f45ptvRFpamvTIyclx1CLYlbXr51514ddw1q6jy5cvi/r164uJEyeK8+fPi7i4OOHn5yfef/99Ry2CXVm7fmJiYkT9+vXF119/LS5duiT27NkjWrVqJV588UVHLYJd5eTkiFOnTolTp04JAOKTTz4Rp06dEn/++acQQoiZM2eKV155RYovGzrg73//uzh79qxYvnx5rR46wNr1s2nTJuHm5iaWL18uO0ZnZWU5ahHsisWSDS1dulQ8+OCDwsPDQ4SGhoqjR49Kr/Xo0UOMGDFCFr9161bx0EMPCQ8PD9GuXTvx3Xff1XDGNcua9dOsWTMBoNwjJiam5hOvIdbuP3erC8WSENavo8OHD4uwsDDh6ekpWrZsKT744ANRUlJSw1nXHGvWT3FxsXj33XdFq1athFqtFkFBQWL8+PHi1q1bNZ94DThw4ECFx5SydTJixAjRo0ePctN06tRJeHh4iJYtW4p169bVeN41xdr106NHD7PxtY1CCJ5vJSIiIqoM+ywRERERmcFiiYiIiMgMFktEREREZrBYIiIiIjKDxRIRERGRGSyWiIiIiMxgsURERERkBoslIiIiIjNYLBFRrSeEwNixY9G4cWMoFAokJSU5OiUiciEsloio1ouPj8f69esRFxeHtLQ0NG3aFOPGjcODDz4IT09P6HQ6REZG4ueff5amad68ORQKBY4ePSqb11tvvYWnnnpKev7uu+9CoVBAoVBApVIhKCgIY8eORWZmZk0tHhHZmZujEyAisreLFy8iICAAXbt2BQB0794dRUVF2LBhA1q2bIn09HQkJCTg5s2bsunUajVmzJiBH374wez827Vrh3379sFoNOLs2bN49dVXYTAYsGXLFrstExHVHBZLRFSrjRw5Ehs2bAAAKBQKaLVaGAwGHDx4ED169AAANGvWDKGhoeWmHTt2LFauXInvv/8effr0qfRvuLm5QafTAQCaNm2Kv/3tb1i3bp0dloaIHIGX4YioVvv000/x3nvv4YEHHkBaWhr++9//QqPR4Ntvv0VhYaHZaVu0aIE33ngDs2bNgslksujv/fHHH9i9ezc8PDxskT4ROQEWS0RUq2m1WtSvXx8qlQo6nQ6+vr5Yv349NmzYgIYNG6Jbt26YPXs2fv311wqnf+edd5CSkoJNmzZV+jdOnz4NjUYDLy8vtGjRAr/99htmzJhhr0UiohrGYomI6pyBAwciNTUVO3fuRFRUFA4ePIjOnTtj/fr15WJ9fX0xbdo0zJ07F0VFRRXO7+GHH0ZSUhKOHz+OGTNmIDIyEpMmTbLzUhBRTWGxRER1klqtxrPPPos5c+bg8OHDGDlyJGJiYiqMjY6Oxu3bt/H5559X+LqHhwdat26N9u3b48MPP4RKpUJsbKw90yeiGsRiiYgIQHBwMPLy8ip8TaPRYM6cOfjggw+Qk5NT5bzeeecdfPzxx0hNTbV1mkTkACyWiKhOuXnzJp555hls3LgRv/76K1JSUrBt2zYsWrQI/fv3r3S6sWPHQqvV4quvvqryb4SHh+PRRx/F/PnzbZk6ETkIiyUiqlM0Gg3CwsKwePFidO/eHe3bt8ecOXMwZswYLFu2rNLp3N3dMW/ePBQUFFj0d6ZMmYIvvvgCV65csVXqROQgCiGEcHQSRERERM6KZ5aIiIiIzGCxRERERGQGiyUiIiIiM1gsEREREZnBYomIiIjIDBZLRERERGawWCIiIiIyg8USERERkRksloiIiIjMYLFEREREZAaLJSIiIiIz/j9CgvCHwddGKQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Make the diagnostic plot of the correlation parameter against the fSNR\n", "plt.scatter(fSNR_indiv, r_indiv)\n", "plt.axhline(r_group[0], color='red', linestyle='--')\n", "plt.axhline(1, color='k', linestyle=':')\n", "plt.axhline(-1, color='k', linestyle=':')\n", "plt.axhline(0, color='k', linestyle=':')\n", "plt.xlabel('fSNR')\n", "plt.ylabel('Correlation parameter')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Figure show the individual correlation estimates plotted against the fSNR, As we can see, the group estimate (red line) is quite close to the true (simulated) value. Individual estimates show a downward bias that and increased variance that becomes worse as the fSNR decreases. Correlation estimates also increasingly hit the boundaries of -1 and 1 as the fSNR decreases." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### One-sample inference via bootstrap \n", "To test whether the estimated group correlation is smaller than a specific value (for example `true_corr=1`, indicating that the activity patterns are the same), we can use a group bootstrap. In this approach, we resample the datasets (participants) with replacement, and for each resampled dataset, we calculate the group correlation estimate. This gives us a distribution of group correlation estimates. If the hypothesized value (x) is outside the 95% confidence interval (for a two-sided test) or 90% confidence interval (for a one-sided test) we can reject the hypothesis that r = x. " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bootstrapping 200 samples: \n", "Bootstrap sample 0/200\n", "Bootstrap sample 50/200\n", "Bootstrap sample 100/200\n", "Bootstrap sample 150/200\n" ] } ], "source": [ "r_boot,fSNR_boot,_ = pcm.bootstrap_group_corr(D,Mflex_sim,fixed_effect=None,n_bootstr=200,verbose=True)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGOCAYAAABFQAMcAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMSxJREFUeJzt3X9wHeV97/HPkakkHCNBbCxhoeAG8PAryIkdG0MT8FSxmmQ0de8k9ZAEGwdoGQwBBHODwT/qZILSUhx3io3yA4eEDoNp5kJKy9gBBd8Mg7i+MYmmJTXGpEbGRMJybiVbGAmsvX8okvVj95zdc3b3eXb3/ZpxiNfnSM95ds/zfPf7/Nic4ziOAAAADCkzXQAAAJBtBCMAAMAoghEAAGAUwQgAADCKYAQAABhFMAIAAIwiGAEAAEYRjAAAAKNOM10AP4aHh/X222/rjDPOUC6XM10cAADgg+M4OnbsmObMmaOyMu/8RyKCkbffflv19fWmiwEAAIpw6NAhnXvuuZ7/nohg5IwzzpA08mGqqqoMlwYAAPjR39+v+vr6sX7cSyKCkdGhmaqqKoIRAAASptAUCyawAgAAowhGAACAUQQjAADAKIIRAABgFMEIAAAwimAEAAAYRTACAACMIhgBAABGEYwAAACjCEYAAIBRgYORX/ziF2pubtacOXOUy+X09NNPF3zP7t279YlPfEIVFRW64IIL9OijjxZRVAAAkEaBg5GBgQE1NDRo69atvl7/X//1X/r85z+vpUuX6te//rXuuOMO3Xjjjdq1a1fgwgIAgPQJ/KC8z372s/rsZz/r+/VtbW364z/+Yz344IOSpIsvvlgvvviivvOd76ipqSnor4/Mrle7te2FA9rfc1zzambolqUXqOnSWs/jlDVY2WzkVlZJsZc/SXUWh6D1EeT1aa3rqD5XUurLRDltrxvbyzdZznEcp+g353J66qmntHz5cs/XfPrTn9YnPvEJbdmyZezYD3/4Q91xxx3q6+vz9Xv6+/tVXV2tvr6+SJ7au+vVbv31Y3snHMvlpJs/fb4e/t9vTDne9pUFxk6qzWX1KpvJ+vLiWlZJk78MUZc/SXUWh6D1EeT1aa3rqD5XUurLRDltrxubyue3/458Amt3d7dqamomHKupqVF/f79OnDjh+p7BwUH19/dP+BOlbS8cmHLMcaQfvXTQ9fi23W9MOR4Xr7I+2nHQ9XicZfUqm8n68uJaVpfXRV3+JNVZHILWR5DXp7Wuo/pcSakvE+W0vW5sL5+bwMM0cWhtbdWmTZti+V3PPPOM/vPtnEbuiyd69/0PXI//5+H/p2eeeSb6wrnwKuuJIfNl9Sqbyfry4lVW19dGWP4k1VkcgtZHkNenta6j+lxJqS8T5bS9boopX3Nzc8Slyi/yzEhtba16enomHOvp6VFVVZVOP/101/esXbtWfX19Y38OHToUbRmnux8v96idczxeH4eoy9p5VHqgM6e7OnJ6oDOnzqOll81kfXnxKqubKMufpDqLQ9D6CPL6tNZ1VJ8rKfVlopy2143t5XMTeTCyZMkStbe3Tzj23HPPacmSJZ7vqaioUFVV1YQ/UVp2rqPcpCR9To6uPsf9+LJzi55mU7Ioy9p5VPrBvjJ1Hc9paDinruM5PbLPf0DiVTaT9eXFrawjAzXxlj9JdRaHoPUR5PVpreuoPldS6stEOW2vG9vL5ybwBNbjx4/rwIGR8aiPf/zj2rx5s5YuXaoPf/jD+shHPqK1a9fq8OHD+vGPfyxpZGnvZZddpjVr1uirX/2qfv7zn+trX/ua/u3f/s33apooJ7COpqw6j0rPvZXT794diR6Xnevo8pnex02KqqwPdI4EIJOdN8PR3Q3+LhMb68uLW1kdxV/+JNVZHILWR5DXp7Wuo/pcSakvE+W0vW6Cli+qYRq//XfgYGT37t1aunTplOOrVq3So48+quuvv14HDx7U7t27J7znzjvv1G9+8xude+65Wr9+va6//vrQP0wxbBjfs8VdHSMZkcnKyxw9uMTeiBoAUBrTwUjgCazXXHON8sUvbrurXnPNNfrVr34V9FchZrXTpa7jU4/bPM4IAEg+nk2DMUkcZwQAJB/BCMY0zJRuuMjReTMclZeN/PfGi+waBwUApI+V+4zAnIaZUsNMMiEAgPiQGQEAAEaRGQFC1HlU+tlbOXW/OzIheNm5jhoY5gKAvAhGgJCMbho3quu49Mi+kXk4BCSAGdwgJAPDNEBIfvbW1D1aHOX0nMtxANErdVdpxIfMCBCS7nfdj//O4ziAYIJmObxvEJiobxsyI0BIkvhwKiApislycIOQHAQjQEjYNA62KOXp27YqZhiUG4TkYJgGocvqhLHRTeOee0slPTwrq/WHcKR1InUxWY5l5zp6ZN9I0DKKGwQ7EYwgVGltCP0qddO4rNdfmpgKKtM6T6KYZ2eFdYOA6BGMIFRpbQjjQv3ZodRAwmRQmdZ5EsVmOdhVOhmYM4JQpbUhjAv1Z14Yy0FNLvNO6zwJnp2VbmRGEKpiUqk4hfozL2h2yi2LYjKoTPM8CbIc6UVmBKFiRUlpqD/zggQSXlmU6nL3nxFHUEkGAUlEZgShYsJYYfnmI5RSf6zCCUeQ7JRXFmXkfx1j2QkyCEgaghGEjobQm5+JjcXUH6twwhNkmMMri9I3lD+oJHAEJiIYAWIU1WoZW1bhpKGTDZKdypdF8QoqCRyBqQhGgBhFNbHRhlU4aepk/WanipksakvgCNiEYASYJMq7+6hWy8S9CsetjrLYyRYzx8eGwBGwDcEIME7Ud/dRLbuMczmnVx1N89hC4/BA6EWwStA5PizftksahhbTgKW9wDhRb1YV1bLLOJdzetVRmUcVfeAoFQ9qCwvLt+0RxgZ3CAeZEWCcOFLoUa02imsVk1cdDTuS5EiaHJWke6gmKJa/2yPf0KLkGM2YlJKxSWK2h2AEnpJ4QZeKFHphXnVU96GRIZkPXGIO5kNMZMvy9yx+x8fzCqwPD5idjF3KcHFSJ5IzTANXWU1fkkIvLF8dzfmQ+3sI5uyT1e/4eF7P8XEbcozr2UJSacPFJp+LVAqCEbjyuqB/vD/djRVbaReWr46iCOY6j0oPdOZ0V0dOD3Sm+/qLU1I7rTB5Xa/DHpdrXBm+UoaLk7pai2EauPK6oIeGc4lI+ZXClhS6zbzqKOz5EElNObuxbUgkqZ1WmLyu111v5YwO15YyXJzUoWaCEbjyuqCl9O8dgdIUu529W0edlr1LbAyqktpphc3tenVk9snHpSzVT+pTmxmmCSBL6WK39OV4Wbp7QrTyzV1Iy927jUMizI/yZnq4tpTfb7rsxSIz4pONdzZRGr2gf7x/ZGhmslLvnmxLWSdNmuovX/bD9N17qfU8+n6vLGMxQVVY554lxvmZHq4t5febLnsxCEZ8Sku6OIiGmdLKeeGn/LIW2IUtbfWXL/sRxfXnZXInP6/a0fOHi6/nyefJTdCgKuxzn8ROC+nEMI1PaUkXBxVFys/GlHUUohrWS1v9eS2vHH3ybRwpZ7ehoucPl1bPbudpvGKCqrSde2AUmRGfTKeLTQr77ikLgV2U2Yu01V+hCXdx3L27Bw7uHbzfevY6T5Kj82YUNySStnMPjCIz4hOTvcKT7044LaK8g01b/dkw4c47cJhqcj17ZcC8ztN5M6S7G4r7fGk798AoghGfbGgw0yILgV2Ud7BprL+GmSMd9INLnKI76lJ4dfIqUM/5VgJFcZ7SeO4BiWGaQJjsFY4szOKPclgvC/UXN6+hosY6R/v7vOs538T2uxuc0M9TmOc+TSuykoD6zo9gBEakPbCLeuOhtNdf3PJ38t717CcD5oz7bxhnLIxzn7YVWbajvgsjGAEiQPYieYrp5PNlwGzugGzfqiBtWQTb69sGBCPAJGFuKpX1hiZtncpk+TJguyzugGxelWNzEFcsm+vbFgQjhqW9sU6aNDaEphRTl0n7PuTLgP1ov/t7bOiAbN6qII1ZBJvr2xYEIwalveNLWscipbMhNCVoXSb1++CVAbO5A7L5YWppzCLYWN+T2+fyj3ar6dJaY+Vhaa9Bad5NMd+SR5ulsSE0xasu3zwu1x1p0/Z9sHkZrs1bFfjZSyVpDy21rb7d2ueb/2mvdr3abaZAIjNiVCkdn+1ZB6+O5acHpZ+9JWvLbfPdbNJ41aWUc816pC0QtH0Ss8k5Tfnar0JZhLRl0ExwbZ8dadvuN4xlR8iMGFTsbopJyDp4dSxH3stZXW6b72aTxq0ux5uc9Ujj7qKmN3OLk99sRaH2q1AWIW0ZNBO82ufXe47FW5BxyIwYVOw4oo3zGibf6VSVS73vFX6f6XJPZuvdrO2ZMDfj6/LN45Lbs17GZz1sHFeHP0GyFX7ar3xZhLRl0IpVSpvglbW8sOaMcAsZAMGIQcV2fLZ9Gd0aolNbPI1veCb/fYRtjYhN6VQpuWlp6VRdPtCZKzj8ZWsgiMKC3CCV2n4xlFp6m+Aa+OekNdecH0VxfSEYMSzsjZbC5Dfy9nri6dmVjqaf5ox1LAMfuGdLstSIFOPpg/ZlwoLym/WwLRAcL4nZqbgECTBKbb/IoJWeHXcL/O/7i4VaZnA1DcFIAsXxZQwSeXs1RH1D0oYFEx8qFrTcWe8AOo9Kve+V9ih7GyQ965Hk7FQcggQYpbZfSb+WwhBGdnxy4G8yEJEIRhIpji9jkMjbb0MUtNx0AF5ZpxFJyyjZnPUoxMZ5WjaZV+2o6/jkOnI0r3pq3YTRfjX84ZlBP3srp9+9K+16KydH2WkX0jhURTCSUFE37EEi7yB3OkHKTQfgfR6UsbS0abbN07LN/j73odr9fZLb4wFLbb+yfqOSxqEqlvbCVZBlllFt6EMH4H0ezq5UptLSpqVx2XGY4v6upnV5r9/l0bZtohYGMiNwFTTyjiJTk8ZUZFBe52H53OTeASVRGu9EwxT3dzWNNypBsz1JHvZ0QzACV3FNEitlJ8YsYLKeHTgP+cX9XU3jjUrUw9K2LwYgGIGnqCPvQncCdAAj0nYHlFRxfB9s7izyifu7GteKQq/zEcW5ijLbk4Q5NgQjMKbUnRiBtEhCZ1FInN/VqIOffOdDiuZcRZntScJiAIIRGJPGcV+gGEnoLGwTZfCT73y4/cYwzlWU2Z4ktLWspoExrFAARiShs8iSfOcjqnMV5QqZJLS1ZEZCkuTxXlOYoAqMSOOEzCTLdz4cRXeuosr2JKGtJRgJgdv44g/2SaflpDkfIjDxwgRVYIQNnQU3VKfkOx+Ogj/WwrQktLU5x3EC1+DWrVv1wAMPqLu7Ww0NDfrHf/xHLVq0yPP1W7Zs0cMPP6yuri7NmjVLX/jCF9Ta2qrKykpfv6+/v1/V1dXq6+tTVVVV0OLm9cwzz5T8M0aeSOq92U5OTqImoqUVjW20qN/SdB6VnvvD9ubnTB/ZYv21vnjqc/INlTS13cra+Z18PsZ33vn+Lamam5sj+bl+++/AmZEdO3aopaVFbW1tWrx4sbZs2aKmpia99tprmj179pTXP/7447rnnnu0fft2XXnlldq/f7+uv/565XI5bd68Oeivt5L3lt0jbJuIlrVGRUrHagWbUb+lG5+ij7s+C02gzeL5zTdkMvpvo23pj/bnMtOWRiXwBNbNmzfrpptu0urVq3XJJZeora1N06dP1/bt211f/9JLL+mqq67Sl770Jc2dO1fLli3Ttddeqz179pRceFt4TQ4az5aJaKONStfxnIaGRzI6j+zz3nY4LWzfPtrvNtC2sr1+kybu+iw0KZPzO1VW29KoBApGhoaGtHfvXjU2Np76AWVlamxsVEdHh+t7rrzySu3du3cs+Pjtb3+rZ599Vp/73Oc8f8/g4KD6+/sn/LHZsnMd5VwXfJ1iy0S0rDYqNq9WSEOjZnP9JlHc9VlotQXnd6qo2tKk35gUK9AwTW9vr06ePKmampoJx2tqarRv3z7X93zpS19Sb2+v/uRP/kSO4+iDDz7QzTffrHvvvdfz97S2tmrTpk1BimbU+MlBhwekDxxJlk5uymqjYvNqhTTsMWFr/SZhSNKtjKXUZzGfudAEWj/lSUJdhymKtjSLw2GjIl9Ns3v3bt1///3atm2bFi9erAMHDuj222/XN7/5Ta1fv971PWvXrlVLS8vY3/v7+1VfXx91UUsyebzX1lnLVeVS73vu/3ZXR3rHPd0aW8nRwAcj58vP542qsU1DgGjDapDJktCwe5XxT+scHTrurz7HX5fV5dKR93JTfl6hz1xotUWh85uEug5bFAF40BsTtzZp9OckLSgMFIzMmjVL06ZNU09Pz4TjPT09qq2tdX3P+vXrdd111+nGG2+UJH3sYx/TwMCA/uqv/kr33XefysqmjhRVVFSooqIiSNGsYusW5p1Hpd733FOIQ8Mjx9PaiIw2tj89OL6xzqn3PX+fN8rG1tasQhA2Lh1MQsbJq4yv9/mrz8nX5RGXGw2/n7nQhM185UlCXYctigA8yI2J15YS47PySWrPAwUj5eXlWrBggdrb27V8+XJJ0vDwsNrb23Xrrbe6vufdd9+dEnBMmzZNklTEqmKUwK3BcJO0RsRvxqJhpvSzt6Ye9/N5o2xsbcwqFMO2IDwJGad8ZfRTn36/02F85nzlyfc50jp8E0UAHuTGxP3cJzcoDDxM09LSolWrVmnhwoVatGiRtmzZooGBAa1evVqStHLlStXV1am1tVXSyNrlzZs36+Mf//jYMM369evV3Nw8FpQgHoWWII9nU4OdT9CMRbEdVJQd2/iszegQ2sxK92dgwL8kZJxKLaPf73TUn9nrc1SXp3v4JuwAPMiNSdra88DByIoVK3TkyBFt2LBB3d3dmj9/vnbu3Dk2qbWrq2tCJmTdunXK5XJat26dDh8+rLPPPlvNzc361re+Fd6ngC9eDYYbmxrsfIJmLIpt/OPo2MaP9fsdPsoqP3fbScg4+ZmLke9z+vlOx/GZvT6Hm6TcqZsQJNuStva8qB1Y42b7DqxxKTXd2XlUemRfbsokzhETG5GwHtAUtbs6cmPzXSZy9JEZU+vIrQ78fN5i3+eX1y6+581wdHeD91c0rSnwfPzsFjr+tbbvlOlVRr+7ov5gX05u6fncuO9AHJ/Z7XP8aL/797O8zNGDS5wJ783adVyqsNvzxO3ACjPCmEDpFXU7snf1TyHedwc51zoqdpw36gmaxQwDZXEFgxQsG2bbPBYvzrj/jv5/P5+zYabUWOfo+cOS6RsKt7r2uyQ4i9dxqdLWnhOMJERYEyi9GuckNNhu3JfsnuJWR8V2UFF2bMUMA5VyTURxJxpG5s7P+5MwMdWvfB2x38/553OluWfYtZJplJ+hsjhW4qQ185Km9pxgJCHibICT9MUdf3fw5nHJLV2dhE6qmPkNXtdE1/H8e8ZEcSda6s8M8v4kTEz1K19HHORz2poB8pNRjGrzsFL3XkG8CEZ8Mt1Bx9UAJzFlOtoQj8y7mPrvcXZSxV4nxQwDeV0TjnIaGvY+d1HciZb6M4O8PwkTU/3K1xGvnJeOz1koUAp7d9cw915BfAI/KC+Lon52iJ9nEbg9/ybMhmm0DI/sS+6za6Kuo0JKvU4aZkp3N4xM7Lu7oXCa3c8zkdzOXRR3oqX+zCDvHw3czpvhqLxs5L9JmXA9Wb5nwqTpc+ZT6Hsb9HsV594rCA+ZER+iHNP0m4mIcgKl26z9yZLwxQ2jjkrJgMW9C+Xkz/v+sPvcmcnnLoosW6k/cySV7v/9YQ5LmMx6Fsry2Dr8Eqawd3e1Ze8VBEMw4kOU8zVsWBng504iKV/cUuqo1CEqExMrx39ev8NUUQxzlLJfRufRiWP6p0Sf1TI9LGnjNvomFLu7qxtb9l6xgenpBUEQjPgQ5XwNG1YGFLqTyMoXt9TMRpTXSZibfEXRAeb7mYU6fK9g+OxKRd4p2/BMlSxkP0rpFIN+rwqtsCsvc7RqXvoDPtOBdlAEIz5EOWHOhpUBXmWIe9Mk00oNDKO6TqIYyouiA/T6mYU6fK967xtyPx7m3Z4NNwNpV2qnGPR7NfGhmCOvHv++LAQikh2BdhAEIz5EmUq1YWWAVxnSNlmu2K21/QaGUV0nUQzlxZm+LdThB6n3sO/2gp7zJKW9bVFqp1jM92r0ezCyK2w2h8CSFmgTjOQRR8Njw5ixDWWImp9OLIzAMIqMQ9iNStzp20IdfpB6D/tuL8jvTlra2xZhXL82blRoOxuy7kEQjHiIs+Gx4QtjQxmi5HdrbRuDsrAblbjTt35WjPit97ADsyC/O2lpb1skrVNMCxuy7kEQjHig4SksSSlrv52YW1Bm+nOG3ajEnb710+H7DYaj6Nj8/u6kpb1tYaJTNP2dtYGtN1deCEY8eDU8hwdGllBm+SKX3DNHP9gnzaqU+ofsq5tiOzEbUvNhNyom7lTDyryZvNvjDr84cXeKNnxnbZGkjDfBiAevhucD59ReDkm+yEu9c3BfjplT7x82riqmbqK8mym2E7MlQxZmo1JKh276jtPk3V7S0t42ibNTtOU7GwbT37c4EYx4cF+r7mjyg9iSeJGHcefgZ5fDIHUT9d1MsZ1YGlPzfutickM4r9rR84fN33F6dWxeDXdYDXrS0t5ZlZbvbNYyPAQjHtwansMD0gcu/WrSLvIw7hz87HIoFX6CbJhlKqSYu7Okp+a9OuJCdeHWEHo9lM+GYNyr4f7TunADqCSlvZMi7Lv/pH9nR6Upw+NH5h+Ul+8hdZMfXDbnQ+4/I2kXeRh3Dn4e0iaNPkG28MOtbL2bMf3wvVKU8uA+r2E4N1GdIz8PkBzl1XD/4nfJffBjFkTxENIkf2fHs7VNjEqmg5Fdr3YH+iKk5SLP96RQvyY/UfTsSkcq4gmyYZYpCsU8OTVIJxol7zurwh2x34eNSdGco6CdlFd5h4bdj6e1QU+aUq5RL2l52rGtbWJUMj1Ms+2FA1OO5UuDpWXMOKyJeJNT1uN3O/T7BNmwyxSFIKn5sMZ5w0hdl3Jn5T0MN3HeVFTnKGiK2qu85WXuAUlaG/Skieru361tStoqSJvbxChkOhjZ3+M+6SHfFyENY8ZRBVXFPEE26jLFLYxx3rACmnxj54WCHa+GsLHO0f6+6M9R0E7Kq7xXn+Po+cPZadCTJo75HUmdCJqWNtGvTAcj82pmqPOtvinHs3DXlC+oCuOuvJioPg2BXhh3emFNXPM6B/OqnYKNc/6G0L49PfKV97wzstOgJ00cd/9JngiahjbRr0wHI7csvUA3P/ZL7prGCesuopSoPslr68O40wsrde11Dnb5bJxNPnQvzGA2Sw160sRx95+1iaBJlXMcx/pvaX9/v6qrq9XX16eqqqrQfm4ul9PpF16h6iu+qD+a9RG939ulvo4ndeLA/wntdyRN7XWbVTFn3pTjg2+/pu7H7or8959+4RWa/T/WTTjmOMM68tT9OvH6y5H//lKdfuEVOvsv7lUudyqgc5xhHflf35pyXY1ce3956tp7+UmdeP3lyM9B/Z0/UVl55ZTjw0MndOg7Xwz0s6I8X3w3EQbTbVrShB0S+O2/M50ZkaQTr7+ct9H06jDS6o9mfSTQ8bBVX/GXU47lcmWqvuKLxuo9yDVw4vWXdeSp+wt2opM78Yo583T2X9yrI0/dr76Xn3QNaPo6ngzl87zf2+XaOL/f2xX4Z0V5vgp9N22TtbYiKaL+PiEcmV7aW8hoh1ExZ57KyivHOozTL7zCdNEi49UhFdNRFcN0MDRZMdfAiddfVvdjd+nQd76o7sfucr2bL9SJH3nqfg2+/ZqGh05o8O3XXDMrxep7+Uk5zsQlJsU2zradL1Oy2FYkRdTfJ4Qj08M0kvTMM894/tvIipCp4+vnzRjZBC2NOo9Kj+zLTRmrj2udvm11HlV57uoY2T9jsvKykU32ojayDDtX8ji9befLFOoBSdfc3BzJz2WYJgRZnPhkejmZbWvro7oGTG9ZnYan6Noki21Fkieawz4EI3mY7jBMMbn6wHQwNFlU10BaOnHbzpcpWWsril11NxrAvD0gleWkYUea8yECGRCM5JWWDiNpbFqKGdU1kKZO3KbzZUrW2opi9u6YHMCMbleTlE3IEC2CkTzS1GGYlOR0bpTXAJ14emStrfAalsr3lG73hy+OSMomZIgOwUgBdBilyfdo9/19yQhQuAbgR5auE69hqZGndLtnOwo9fDHN82tsYuvNIcEIIr04vdK5zx8+dZw0bXrY2tAhXG7DUpNNznZ4P3xxRJzza7J6neab6xPNWhr/2Gck44I+qj0ov4+iL/Wx4TAv6mtp/O95oDOnuzpyeqAz/J+PwkaHpc6b4ai8zFHO43lF47Mdy871fl2c82viuk5t5D3Xx3zbSzASkaQ0mFFfnLUB7nZI0yZbHA1dljsS2zTMlO5uGNkXp36G+2vGZzvGBzCn5UaCmNNyI3+Pax8jye4OOWo2L0FnmCYCSXpkddQXp3s615Fc0rtpXQaZFXE0dEl+Amua+V1NZMO8Gps75KjZvASdzEgEkhR5e2Uuwro4J6dzz5vh6DN1U9O1aV4GmRVRX0tStjsSm7l9z+PMdgQRx3VqK7ehMlvaXjIjEUhSgxnH/ghud0PnnZGdZZBZEce1ZPOdXdbZkPXwI2t7woxn8xJ0gpEIFNtgmpjhberiTErDBf/iuJay3JEgHDZ3yKOi7AtsbXt5UF6eB+UVq5iHzU3ZnfAP77FxnglgUlgP+QNsZKov4EF5KVRM5M3EPMAfW+/sgDBktS8gGIlI0AYzSfNMANOyumkV0i+rfQGraSyR5RneQBDsNYI0qy4PdjwtyIxYgol5yJpisxtZTWMjG7J6BROMWCIJM7yBsJSyMWBW09hJxrCaf/1D7sf7PI6nBcGIRZiYh6woJbvBXiPJkqQdqW2Q1eubOSMAYldKdsPmXSQxVZJ2pLZBVq9vMiNIDFK95oV1Dkq5+0v6kGbWrmOG1YJJ+vVdLIIRJAKpXvPCPAelTthO6pBmFq/jrA47lCKp13cpGKZBIpDqNS/Mc5CkB6uFKYvXcVaHHRAMmREkAqle88I+B1m8+8vidZzVYQcEQzCCRCDVax7noHRV5VLve1OPp70Osxh4IhiGaZAIpHrN4xyUpvOo1Pue23AMdQgQjCARsjrHwCacg9K4zReRpLMrRR0i8ximQWKQ6jWPc1A8r/kiad9ZE/CDYAQAYpC0OTdZ2w8FZhGMAEAMwn4YZpTBQhb3Q4FZzBkBgBiEOedmNFjoOp7T0HBOXcdzemRfTp1HwylrFvdDgVlkRgAgJmHNuSnlQYN+ZHE/FJhVVGZk69atmjt3riorK7V48WLt2bMn7+v/+7//W2vWrNE555yjiooKzZs3T88++2xRBUa6dR6VHujM6a6OnB7oDO9OD0iTqIOFWo95LLbOb0HyBc6M7NixQy0tLWpra9PixYu1ZcsWNTU16bXXXtPs2bOnvH5oaEif+cxnNHv2bP3kJz9RXV2d3nzzTZ155plhlB8WKnYsm3FqwJ9iJsMG+V6GPb8FKCRwZmTz5s266aabtHr1al1yySVqa2vT9OnTtX37dtfXb9++Xb///e/19NNP66qrrtLcuXN19dVXq6GhoeTCwz6ljGUzTm0PMlR2C7oBXdDvJXvKIG6BgpGhoSHt3btXjY2Np35AWZkaGxvV0dHh+p5/+Zd/0ZIlS7RmzRrV1NTosssu0/3336+TJ096/p7BwUH19/dP+INkKCWgYJzaDlFPjkTpggYLxXwvG2ZKdzc4enCJo7sbCEQQrUDDNL29vTp58qRqamomHK+pqdG+fftc3/Pb3/5WP//5z/XlL39Zzz77rA4cOKBbbrlF77//vjZu3Oj6ntbWVm3atClI0YrW3Nwcy+/Jiv+5Z6ekqYHmkaHT1Nz8Z3nf+4M3X1TnW31Tjl9cd5aam68Kq4go4AcPvShp4nlwlNP/HThL667nPNiiWdI6n68t5XsJxCHypb3Dw8OaPXu2vve972nBggVasWKF7rvvPrW1tXm+Z+3aterr6xv7c+jQoaiLiZDMq5nhevzCmjMKvveWpRcoN+lGLZeT1lxzfhhFg0/7e1wmI0h6vedYzCVBWEr5XgJxCBSMzJo1S9OmTVNPT8+E4z09PaqtrXV9zznnnKN58+Zp2rRpY8cuvvhidXd3a2jIfR/kiooKVVVVTfiDZCgloGi6tFZtX1mghvozNb18mhrqz9R3v7JAyy51v7YQDTqu9CHQh+0CBSPl5eVasGCB2tvbx44NDw+rvb1dS5YscX3PVVddpQMHDmh4eHjs2P79+3XOOeeovLy8yGLDVqUGFE2X1uqna67Sb77xZ/rpmqsIRAyg40qfJAX6u17t1p8/9KIuXr9Tf/7Qi9r1arfpIiEGOcdxAq3V2rFjh1atWqXvfve7WrRokbZs2aInn3xS+/btU01NjVauXKm6ujq1trZKkg4dOqRLL71Uq1at0m233abXX39dX/3qV/W1r31N9913n6/f2d/fr+rqavX19ZElAWKw69Vubdv9hl7vOaYLa87QmmvOt7LjQrrserVbf/3Y3gnHcjmp7SsL1MT1l0h+++/A+4ysWLFCR44c0YYNG9Td3a358+dr586dY5Nau7q6VFZ2KuFSX1+vXbt26c4779Tll1+uuro63X777fr6179exMcCEIemS2tp/BG7bS8cmHLMcaRtu9/geky5wJkRE8iMAED6Xbx+p068P3XVz/TyafrNN1j1k0R++28elAcAsAKTp7OLYAQAYAUmT2cXwQgAwApJWvWDcAWewAoAQFSYPJ1NZEYAAIBRBCMAAMAoghEAAGAUwQgAADCKYAQAABhFMAIAAIwiGAEAAEYRjAAAAKMIRgAAgFEEIwAAwCiCEQAAYBTBCAAAMIoH5QEAUMCuV7u17YUD2t9zXPNqZuiWpRfwQL8QkRkBACCPXa92668f26vOt/p04v2T6nyrTzf/017terXbdNFSg2AEAIA8tr1wYMoxx5G27X7DQGnSiWAEAIA89vccdz3+es+xmEuSXgQjAADkMa9mhuvxC2vOiLkk6UUwAgBAHrcsvUC53MRjuZy05przzRQohQhGAADIo+nSWrV9ZYEa6s/U9PJpaqg/U9/9ygItYzVNaFjaCwBAAU2X1rKUN0IEIwCQIOx3gTRimAYAEoL9LpBWBCMAkBDsd4G0IhgBgIRgvwukFcEIACQE+10grQhGACAh2O8CaUUwAgAJwX4XSCuW9gJAgrDfBdKIzAgAADCKYAQAABhFMAIAAIwiGAEAAEYRjAAAAKMIRgAAgFEEIwAAwCiCEQAAYBTBCAAAMIpgBAAAGEUwAgAAjCIYAQAARhGMAAAAowhGAACAUQQjAADAKIIRAABgFMEIAAAwimAEAAAYRTACAACMIhgBAABGEYwAAACjCEYAAIBRBCMAAMAoghEAAGAUwQgAADCKYAQAABhFMAIAAIwqKhjZunWr5s6dq8rKSi1evFh79uzx9b4nnnhCuVxOy5cvL+bXAgCAFAocjOzYsUMtLS3auHGjXnnlFTU0NKipqUnvvPNO3vcdPHhQd999tz71qU8VXVgAAJA+gYORzZs366abbtLq1at1ySWXqK2tTdOnT9f27ds933Py5El9+ctf1qZNm/TRj360pAIDAIB0CRSMDA0Nae/evWpsbDz1A8rK1NjYqI6ODs/3feMb39Ds2bN1ww03FF9SAACQSqcFeXFvb69OnjypmpqaCcdramq0b98+1/e8+OKLeuSRR/TrX//a9+8ZHBzU4ODg2N/7+/uDFBMAACRIpKtpjh07puuuu07f//73NWvWLN/va21tVXV19dif+vr6CEsJAABMCpQZmTVrlqZNm6aenp4Jx3t6elRbWzvl9W+88YYOHjyo5ubmsWPDw8Mjv/i00/Taa6/p/PPPn/K+tWvXqqWlZezv/f39BCQAAKRUoGCkvLxcCxYsUHt7+9jy3OHhYbW3t+vWW2+d8vqLLrpI//7v/z7h2Lp163Ts2DH9wz/8g2eAUVFRoYqKiiBFAwAACRUoGJGklpYWrVq1SgsXLtSiRYu0ZcsWDQwMaPXq1ZKklStXqq6uTq2traqsrNRll1024f1nnnmmJE05DgAAsilwMLJixQodOXJEGzZsUHd3t+bPn6+dO3eOTWrt6upSWRkbuwIAAH9yjuM4pgtRSH9/v6qrq9XX16eqqirTxQEAAD747b9JYQAAAKMIRgAAgFEEIwAAwCiCEQAAYBTBCAAAMIpgBAAAGEUwAgAAjCIYAQAARhGMAAAAowhGAACAUQQjAADAKIIRAABgFMEIAAAwimAEAAAYRTACAACMIhgBAABGEYwAAACjCEYAAIBRBCMAAMAoghEAAGAUwQgAADCKYAQAABhFMAIAAIwiGAEAAEYRjAAAAKMIRgAAgFEEIwAAwCiCEQAAYBTBCAAAMIpgBAAAGEUwAgAAjCIYAQAARhGMAAAAowhGAACAUQQjAADAKIIRAABgFMEIAAAwimAEAAAYRTACAACMIhgBAABGEYwAAACjCEYAAIBRBCMAAMAoghEAAGAUwQgAADCKYAQAABhFMAIAAIwiGAEAAEYRjAAAAKMIRgAAgFEEIwAAwCiCEQAAYBTBCAAAMIpgBAAAGEUwAgAAjCIYAQAARhGMAAAAowhGAACAUQQjAADAKIIRAABgVFHByNatWzV37lxVVlZq8eLF2rNnj+drv//97+tTn/qUzjrrLJ111llqbGzM+3oAAJAtgYORHTt2qKWlRRs3btQrr7yihoYGNTU16Z133nF9/e7du3XttdfqhRdeUEdHh+rr67Vs2TIdPny45MIDAIDkyzmO4wR5w+LFi/XJT35SDz30kCRpeHhY9fX1uu2223TPPfcUfP/Jkyd11lln6aGHHtLKlSt9/c7+/n5VV1err69PVVVVQYoLAAAM8dt/B8qMDA0Nae/evWpsbDz1A8rK1NjYqI6ODl8/491339X777+vD3/4w56vGRwcVH9//4Q/AAAgnQIFI729vTp58qRqamomHK+pqVF3d7evn/H1r39dc+bMmRDQTNba2qrq6uqxP/X19UGKCQAAEiTW1TTf/va39cQTT+ipp55SZWWl5+vWrl2rvr6+sT+HDh2KsZQAACBOpwV58axZszRt2jT19PRMON7T06Pa2tq87/37v/97ffvb39bzzz+vyy+/PO9rKyoqVFFREaRoAAAgoQJlRsrLy7VgwQK1t7ePHRseHlZ7e7uWLFni+b6/+7u/0ze/+U3t3LlTCxcuLL60AAAgdQJlRiSppaVFq1at0sKFC7Vo0SJt2bJFAwMDWr16tSRp5cqVqqurU2trqyTpb//2b7VhwwY9/vjjmjt37tjckhkzZmjGjBkhfhQAAJBEgYORFStW6MiRI9qwYYO6u7s1f/587dy5c2xSa1dXl8rKTiVcHn74YQ0NDekLX/jChJ+zceNG/c3f/E1ppQcAAIkXeJ8RE9hnBACA5IlknxEAAICwEYwAAACjCEYAAIBRBCMAAMAoghEAAGAUwQgAADCKYAQAABhFMAIAAIwiGAEAAEYRjAAAAKMIRgAAgFEEIwAAwCiCEQAAYBTBCAAAMIpgBAAAGEUwAgAAjCIYAQAARhGMAAAAowhGAACAUQQjAADAKIIRAABgFMEIAAAwimAEAAAYRTACAACMIhgBAABGEYwAAACjCEYAAIBRBCMAAMAoghEAAGAUwQgAADCKYAQAABhFMAIAAIwiGAEAAEYRjAAAAKMIRgAAgFEEIwAAwCiCEQAAYBTBCAAAMIpgBAAAGEUwAgAAjCIYAQAARhGMAAAAowhGAACAUQQjAADAKIIRAABgFMEIAAAwimAEAAAYRTACAACMIhgBAABGEYwAAACjCEYAAIBRBCMAAMAoghEAAGAUwQgAADCKYAQAABhFMAIAAIwiGAEAAEYRjAAAAKMIRgAAgFFFBSNbt27V3LlzVVlZqcWLF2vPnj15X//P//zPuuiii1RZWamPfexjevbZZ4sqLAAASJ/AwciOHTvU0tKijRs36pVXXlFDQ4Oampr0zjvvuL7+pZde0rXXXqsbbrhBv/rVr7R8+XItX75c//Ef/1Fy4QEAQPLlHMdxgrxh8eLF+uQnP6mHHnpIkjQ8PKz6+nrddtttuueee6a8fsWKFRoYGNC//uu/jh274oorNH/+fLW1tfn6nf39/aqurlZfX5+qqqqCFBcAABjit/8OlBkZGhrS3r171djYeOoHlJWpsbFRHR0dru/p6OiY8HpJampq8nw9AADIltOCvLi3t1cnT55UTU3NhOM1NTXat2+f63u6u7tdX9/d3e35ewYHBzU4ODj2976+PkkjERYAAEiG0X670CBMoGAkLq2trdq0adOU4/X19QZKAwAASnHs2DFVV1d7/nugYGTWrFmaNm2aenp6Jhzv6elRbW2t63tqa2sDvV6S1q5dq5aWlrG/Dw8P6/e//71mzpypXC4XpMgALNff36/6+nodOnSIOWFAyjiOo2PHjmnOnDl5XxcoGCkvL9eCBQvU3t6u5cuXSxoJFNrb23Xrrbe6vmfJkiVqb2/XHXfcMXbsueee05IlSzx/T0VFhSoqKiYcO/PMM4MUFUDCVFVVEYwAKZQvIzIq8DBNS0uLVq1apYULF2rRokXasmWLBgYGtHr1aknSypUrVVdXp9bWVknS7bffrquvvloPPvigPv/5z+uJJ57QL3/5S33ve98L+qsBAEAKBQ5GVqxYoSNHjmjDhg3q7u7W/PnztXPnzrFJql1dXSorO7VI58orr9Tjjz+udevW6d5779WFF16op59+Wpdddll4nwIAACRW4H1GACBMg4ODam1t1dq1a6cMzwLIBoIRAABgFA/KAwAARhGMAAAAowhGAACAUQQjAADAKIIRAABgFMEIAAAwimAEAAAYRTACAACMIhgBAABGEYwAAACjCEYAAIBR/x+xhbQ8KVxNoQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ax=sns.stripplot(y=r_boot,dodge=True,legend=False)\n", "ax.set_ylim([0,1.05])\n", "x_coord = ax.get_xticks()\n", "CI_low = np.percentile(r_boot, 2.5)\n", "CI_high = np.percentile(r_boot, 97.5)\n", "pcm.plot_CI(ax,x_coord[0],center=r_group,low=CI_low,high=CI_high,width=0.2,facecolor='black',alpha=0.3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Figure shows the distribution of the group correlation estimates obtained from the bootstrap procedure. The black line indicates the group estimate obtained from the original data, and the shaded area indicates the 95% confidence interval obtained from the bootstrap distribution. As we can see, the true value of 0.7 is well within the confidence interval, so we cannot reject the hypothesis that r = 0.7. However, if we had hypothesized that r = 1, we would have been able to reject this hypothesis, as 1 is outside the confidence interval. Note the the lower bound of the confidence interval tends to be a bit too high, so testing whether r is larger than a specific value needs to be done with caution (Diedrichsen et al., 2026)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Paired-sample inference via bootstrap\n", "We can extend this bootstrap approach to test whether the correlation is higher in one brain area than another. In this case, we fit a separate models for each brain area. For the boostrap procedure, we resample the datasets (participants) with replacement, and for each resampled dataset, we calculate the group correlation estimate for both brain areas. This gives us a distribution of the difference in group correlation estimates between the two brain areas. " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "n_subj = 30\n", "rng = np.random.default_rng(seed=102)\n", "# Data set for brain area 1 with low signal\n", "signal = np.random.gamma(shape=2,scale=0.03,size=(n_subj,))\n", "D1 = pcm.sim.make_dataset(model=Mtrue_sim, \\\n", " theta=[0,0], \\\n", " cond_vec=cond_vec, \\\n", " part_vec=part_vec, \\\n", " n_sim=n_subj, \\\n", " signal=signal,\\\n", " rng=rng)\n", "# Data set for brain area 2 with higher signal, but lower correlation\n", "Mtrue_sim = pcm.CorrelationModel('corr', num_items=1, corr=0.5, cond_effect=False)\n", "signal = np.random.gamma(shape=2,scale=0.08,size=(n_subj,))\n", "D2 = pcm.sim.make_dataset(model=Mtrue_sim, \\\n", " theta=[0,0], \\\n", " cond_vec=cond_vec, \\\n", " part_vec=part_vec, \\\n", " n_sim=n_subj, \\\n", " signal=signal,\\\n", " rng=rng)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated group correlation for the different areas are: 0.8293 , 0.4975\n", "Estimated fSNR for the different areas are: 0.2947 , 1.0843\n" ] } ], "source": [ "# Get the group estimates for the two areas and compare the estimated correlation parameters:\n", "T_gr1, theta_gr1 = pcm.fit_model_group(D1, Mflex_sim, fixed_effect=None, fit_scale=True, verbose=False)\n", "theta_g1,_= pcm.group_to_individ_param(theta_gr1[0],Mflex_sim,n_subj)\n", "r_group1 = Mflex_sim.get_correlation(theta_g1).mean()\n", "fSNR_group1 = Mflex_sim.get_fSNR(theta_g1, n_part, separate=False).mean()\n", "T_gr2, theta_gr2 = pcm.fit_model_group(D2, Mflex_sim, fixed_effect=None, fit_scale=True, verbose=False)\n", "theta_g2,_= pcm.group_to_individ_param(theta_gr2[0],Mflex_sim,n_subj)\n", "r_group2 = Mflex_sim.get_correlation(theta_g2).mean()\n", "fSNR_group2 = Mflex_sim.get_fSNR(theta_g2, n_part, separate=False).mean()\n", "print(f\"Estimated group correlation for the different areas are: {r_group1:.4f} , {r_group2:.4f}\")\n", "print(f\"Estimated fSNR for the different areas are: {fSNR_group1:.4f} , {fSNR_group2:.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the fSNR is much lower for the first dataset, so if we calculated Pearson correlations, we would have concluded that the correlation is higher in the second dataset. However, when we take into account the measurement noise, we see that the estimated group correlation is actually higher for the first dataset.\n", "\n", "To test this hypothesis we are running the bootstrap for each dataset. The bootstrap for the second dataset will be done with the same bootstrap indices than for the first dataset. We can use the distribution of the difference in correlation of those paired bootstrap samples to test whether the correlation of D1 > D2. \n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bootstrapping 200 samples: \n", "Bootstrap sample 0/200\n", "Bootstrap sample 50/200\n", "Bootstrap sample 100/200\n", "Bootstrap sample 150/200\n", "Bootstrapping 200 samples: \n", "Bootstrap sample 0/200\n", "Bootstrap sample 50/200\n", "Bootstrap sample 100/200\n", "Bootstrap sample 150/200\n" ] } ], "source": [ "r_boot1,fSNR_boot1,boot_indx = pcm.bootstrap_group_corr(D1,Mflex_sim,fixed_effect=None,n_bootstr=200,verbose=True)\n", "r_boot2,fSNR_boot2,boot_indx = pcm.bootstrap_group_corr(D2,Mflex_sim,fixed_effect=None,n_bootstr=200,verbose=True,boot_indx=boot_indx)\n", "diff_boot = r_boot1 - r_boot2" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "95% confidence interval for the difference in correlation is: [0.1322, 0.4859]\n" ] } ], "source": [ "CI_low = np.percentile(diff_boot, 2.5)\n", "CI_high = np.percentile(diff_boot, 97.5)\n", "print(f\"95% confidence interval for the difference in correlation is: [{CI_low:.4f}, {CI_high:.4f}]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given these results, we can conclude that the correlation in D1 is significantly different from D2. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Multi-pattern setting: Estimating the correlation between two sets of patterns across two conditions\n", "In the second part of this notebook, we will simulate data from a hypothetical experiment, in which participants observed 3 hand gestures or executed the same 3 hand gestures. Thus, we have 3 items (i.e., the hand gestures) in each of 2 conditions (i.e., either observe or execute). \n", "\n", "We are interested in the average correlation between the patterns associated with observing and executing action A, observing and executing action B, and observing and executing action C, while accounting for overall differences in the average patterns of observing and executing.\n", "\n", "For the condition effect, we have two options: We can either model the condition effect as a random effect with an associated pattern variance. We will take this approach when generating the data. We can also model the condition effect as a fixed effect, which is equivalent to removing the mean pattern for each condition from the data - this approach is taken here when fitting the model. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulating data\n", "First, we create our true model (`Mtrue`): one where the all actions are equally strongly encoded in each condition, but where the strength of encoding can differ between conditions (i.e., between observation or execution). \n", "\n", ">For example, we could expect the difference between actions to be smaller during observation than during execution (simply due to overall levels of brain activation). \n", "\n", "Next, we also model the covariance between items within each condition with a condition effect (i.e., by setting `condEffect` to `True`). Finally, we set the ground-truth correlation to be 0.7. " ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The true model has 4 hyper parameters\n" ] } ], "source": [ "Mtrue = pcm.CorrelationModel('corr', num_items=3, corr=0.7, cond_effect=True, within_cov=None)\n", "print(f\"The true model has {Mtrue.n_param} hyper parameters\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These four parameters are concerned with the condition effect and item effect for observation and execution, respectively. Visualizing the components of the second moment matrix (also known as variance-covariance, or simply covariance matrix) helps to understand this:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAACXCAYAAABJNBKHAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAADDpJREFUeJzt3U9oXPUaxvFnkpiJeifTamiaMdHaRSham2A1oXgVwcHgdXF1lUUXEkQ3o1jiQrIxuoorEbRUESQbJXEjBZEUCcRSaWlJKCiC10q5DsRJLFxmpl1M28zvLrwdzU3G5kze82dOvx840JyenvN2+uTXp2emMwnnnBMAAICBlrAHAAAA8UGxAAAAZigWAADADMUCAACYoVgAAAAzFAsAAGCGYgEAAMy0BXmxarWq5eVlpVIpJRKJIC+NBjjnVC6Xlclk1NJi10HJQXPxIwdkoLmwFsBLBgItFsvLy+rr6wvykjCQz+fV29trdj5y0Jwsc0AGmhNrAbaSgUCLRSqVkiT9e2mPOv8Wr2dhnu9/KOwRzF3XNZ3SV7U/Nys3zvd3/UNtus303LDnRw4s14I4fu9FDWsBvGQg0GJx41ZX599a1JmKV7FoS8Twm+J/b/ZufYvyxvnadFs8H7e48SEHlmsBGQoAawE8ZCBef7sDAIBQUSwAAICZhorF0aNHtWfPHnV0dGh4eFhnz561ngsRRwYgkQOQAWzkuVjMzs5qfHxck5OTWlpa0sDAgEZGRrS6uurHfIggMgCJHIAMYHOei8W7776rl156SWNjY3rggQf04Ycf6o477tAnn3zix3yIIDIAiRyADGBznorF1atXtbi4qGw2+8cJWlqUzWZ1+vTpDcdXKhWVSqV1G5qb1wxI5CCOWAvAWoB6PBWLS5cuaW1tTd3d3ev2d3d3q1AobDh+ampK6XS6tvFGKM3PawYkchBHrAVgLUA9vv6vkImJCRWLxdqWz+f9vBwiihyADEAiB7cKT2+Q1dXVpdbWVq2srKzbv7Kyot27d284PplMKplMbm9CRIrXDEjkII5YC8BagHo83bFob2/XwYMHNT8/X9tXrVY1Pz+vQ4cOmQ+H6CEDkMgByADq8/yW3uPj43rhhRf0yCOPaGhoSO+9956uXLmisbExP+ZDBJEBSOQAZACb81wsRkdH9dtvv+nNN99UoVDQ4OCg5ubmNryAB/FFBiCRA5ABbC7hnHNBXaxUKimdTus//9obuw8hG8kMhj2CuevumhZ0XMViUZ2dnWbnvZGDJ/VPPnioCfiRA8u1II7fe1HDWgAvGYjX3+4AACBUFAsAAGDG82ssAMDK8/0PbfsW+Inl8yaz8JQKYIM7FgAAwAzFAgAAmKFYAAAAMxQLAABghmIBAADMUCwAAIAZigUAADBDsQAAAGYoFgAAwAzFAgAAmKFYAAAAMxQLAABghmIBAADMUCwAAIAZigUAADBDsQAAAGYoFgAAwExb2AMAwHaMZAZNznNi+bzJeazmAZoVdywAAIAZigUAADBDsQAAAGYoFgAAwIynYjE1NaVHH31UqVRKu3bt0nPPPacff/zRr9kQQWQAEjkAGUB9norFN998o1wupzNnzujrr7/WtWvX9PTTT+vKlSt+zYeIIQOQyAHIAOrz9N9N5+bm1n09PT2tXbt2aXFxUU888YTpYIgmMgCJHIAMoL5tvY9FsViUJN11112b/nylUlGlUql9XSqVtnM5RNDNMiCRg1sBawFYC3BDwy/erFarOnLkiB577DHt379/02OmpqaUTqdrW19fX8ODInq2kgGJHMQdawFYC/BnDReLXC6n77//XjMzM3WPmZiYULFYrG35fL7RyyGCtpIBiRzEHWsBWAvwZw09FfLKK6/oyy+/1MmTJ9Xb21v3uGQyqWQy2fBwiK6tZkAiB3HGWgDWAvw/T8XCOadXX31VX3zxhRYWFnT//ff7NRciigxAIgcgA6jPU7HI5XL67LPPdPz4caVSKRUKBUlSOp3W7bff7suAiBYyAIkcgAygPk+vsTh27JiKxaKefPJJ9fT01LbZ2Vm/5kPEkAFI5ABkAPV5fioEtzYyAIkcgAygPj4rBAAAmKFYAAAAM9t6581GPd//kNoSt4VxacBXJ5bPhz2CuVK5qp39YU/hv5HMYNgjIEas1oJmzCV3LAAAgBmKBQAAMEOxAAAAZigWAADADMUCAACYoVgAAAAzFAsAAGCGYgEAAMxQLAAAgBmKBQAAMEOxAAAAZigWAADADMUCAACYoVgAAAAzFAsAAGCGYgEAAMxQLAAAgJm2sAcAACBuRjKDJuc5sXze5DxW82wFdywAAIAZigUAADBDsQAAAGYoFgAAwMy2isU777yjRCKhI0eOGI2DZkMGQAYgkQP8oeFice7cOX300Uc6cOCA5TxoImQAZAASOcB6DRWLy5cv6/Dhw/r444+1c+dO65nQBMgAyAAkcoCNGioWuVxOzz77rLLZ7F8eV6lUVCqV1m2Ih61mQCIHcUUGIJEDbOT5DbJmZma0tLSkc+fO3fTYqakpvf322w0NhujykgGJHMQRGYBEDrA5T3cs8vm8XnvtNX366afq6Oi46fETExMqFou1LZ/PNzwoosFrBiRyEDdkABI5QH2e7lgsLi5qdXVVDz/8cG3f2tqaTp48qQ8++ECVSkWtra21n0smk0omk3bTInReMyCRg7ghA5DIAerzVCyeeuopfffdd+v2jY2Nad++fXrjjTc2hAjxQwZABiCRA9TnqVikUint379/3b4777xTd99994b9iCcyADIAiRygPt55EwAAmNn2x6YvLCwYjIFmRgZABiCRA/yOOxYAAMAMxQIAAJjZ9lMhAADAHyOZQZPznFg+v61fXypXtbN/a8dyxwIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMAMxQIAAJihWAAAADMUCwAAYIZiAQAAzFAsAACAGYoFAAAwQ7EAAABmKBYAAMBMW5AXc85Jkq7rmuSCvDIacV3XJP3x52Ylzjkolathj2CudPn335NlDuKcgThiLWh+212bvKwDgRaLcrksSTqlr4K8LLapXC4rnU6bnk+KZw529oc9gX8scxDnDMQZa0HzslqbtpKBhLOuoH+hWq1qeXlZqVRKiURi02NKpZL6+vqUz+fV2dkZ1Gi3nK08zs45lctlZTIZtbTYPWtGDqIjrByQgehgLYB1BgK9Y9HS0qLe3t4tHdvZ2UmIAnCzx9nyXyc3kIPoCToHZCB6WAtglQFevAkAAMxQLAAAgJnIFYtkMqnJyUklk8mwR4m1qD/OUZ8vLqL8OEd5tjiJ+uMc9fniwPoxDvTFmwAAIN4id8cCAAA0L4oFAAAwQ7EAAABmKBYAAMAMxQIAAJiJXLE4evSo9uzZo46ODg0PD+vs2bNhjxQrb731lhKJxLpt3759YY+1DhnwVzNkQCIHfmuGHJABf/mVgUgVi9nZWY2Pj2tyclJLS0saGBjQyMiIVldXwx4tVh588EH9+uuvte3UqVNhj1RDBoIR5QxI5CAoUc4BGQiGLxlwETI0NORyuVzt67W1NZfJZNzU1FSIU8XL5OSkGxgYCHuMusiA/6KeAefIQRCingMy4D+/MhCZOxZXr17V4uKistlsbV9LS4uy2axOnz4d4mTx89NPPymTyWjv3r06fPiwfvnll7BHkkQGghTVDEjkIEhRzQEZCI4fGYhMsbh06ZLW1tbU3d29bn93d7cKhUJIU8XP8PCwpqenNTc3p2PHjunixYt6/PHHVS6Xwx6NDAQkyhmQyEFQopwDMhAMvzIQ6MemI3zPPPNM7ccHDhzQ8PCw7rvvPn3++ed68cUXQ5wMQSEDkMgB/MtAZO5YdHV1qbW1VSsrK+v2r6ysaPfu3SFNFX87duxQf3+/Lly4EPYoZCAkUcqARA7CEqUckIFwWGUgMsWivb1dBw8e1Pz8fG1ftVrV/Py8Dh06FOJk8Xb58mX9/PPP6unpCXsUMhCSKGVAIgdhiVIOyEA4zDJg/nLQbZiZmXHJZNJNT0+7H374wb388stux44drlAohD1abLz++utuYWHBXbx40X377bcum826rq4ut7q6GvZozjkyEISoZ8A5chCEqOeADPjPrwxEqlg459z777/v7r33Xtfe3u6GhobcmTNnwh4pVkZHR11PT49rb29399xzjxsdHXUXLlwIe6x1yIC/miEDzpEDvzVDDsiAv/zKQMI552xuogAAgFtdZF5jAQAAmh/FAgAAmKFYAAAAMxQLAABghmIBAADMUCwAAIAZigUAADBDsQAAAGYoFgAAwAzFAgAAmKFYAAAAM/8FDiS319Omrp0AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "H = Mtrue.n_param\n", "for i in range(H):\n", " plt.subplot(1, H, i+1)\n", " plt.imshow(Mtrue.Gc[i,:,:])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">The first two components plotted above reflect the condition effect and model the covariance between items within each condition (observation, execution). The second two components reflect the item effect and model the item-specific variance for each item (3 hand gestures) in each condition." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To Simulate a dataset, we need to simulate an experimental design. Let's assume we measure the 6 trial types (3 items x 2 conditions) in 8 imaging runs and submit the beta-values from each run to the model as $\\mathbf{Y}$. \n", "\n", "We then generate a dataset where there is a strong overall effect for both observation (exp(0)) and execution (exp(1)). In comparison, the item-specific effects for observation (exp(-1.5)) and execution (exp(-1)) are pretty weak (this is a rather typical for real fMRI data).\n", "\n", ">Note that all hyper parameters are log(variances) — this ensures that the variances positive and the math easy. " ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "# Create the design. In this case it's 8 runs, 6 trial types\n", "n_part = 8\n", "cond_vec, part_vec = pcm.sim.make_design(n_cond=6, n_part=n_part)\n", "#print(cond_vec)\n", "#print(part_vec)\n", "\n", "# Starting from the true model above, generate 20 datasets/participants\n", "theta_true = [0,1,-1.5,-1]\n", "D = pcm.sim.make_dataset(model=Mtrue, \\\n", " theta=theta_true, \\\n", " cond_vec=cond_vec, \\\n", " part_vec=part_vec, \\\n", " n_sim=n_subj, \\\n", " signal=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inspecting the data\n", "As a quick check, let's plot the predicted second moment matrix of our true model (using the simulation parameters) and the crossvalidated estimate from the first dataset. " ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, 'dataset')" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAEjCAYAAABuGEhQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHhRJREFUeJzt3XtwlPXd/vFrk5gNhCQKJIE8BMPBiomCPzkVqQJCwZSCtFqp2hoQta2JQqkdh+mUgK2GTmcoKhjRUWitDB54wMNoGA4C2oLG8KAclAcVMYgkHMomBNjA7vf3Rx9WtxtINnyTew/v18zOZO98d+/Ljfvh2ntPLmOMEQAAgAUJTgcAAACxg2IBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZigXPasGGDXC6XNmzYENg2ZcoU5eXlOZbpPzWVEUDrzJkzRy6Xy+kYiHIUC7SLRx99VKtWrXI6BoA28OSTT2rp0qVOx5AkHThwQHPmzNG2bducjhK3KBYIyzPPPKPdu3eHfTmKBRC7Iq1YzJ07l2LhoCSnA8A+v9+vxsZGpaSkWL/uiy66yPp1AgBiB0csItjZ5zs/+eQT3XrrrUpPT1eXLl00ffp0nTp1KrDO5XKppKREL7zwggoKCuR2u1VRUSFJ+uqrr3TXXXcpOztbbrdbBQUFeu6550L2tX//fk2aNEmpqanKysrSr3/9a3m93pB1Tb3Gwu/367HHHtNVV12llJQUZWZm6sYbb9QHH3wQyNfQ0KC//vWvcrlccrlcmjJlSuDytjMCaN67776rwYMHKyUlRX369NHixYtD1ixZskQ33HCDsrKy5Ha7lZ+fr/Ly8qA1eXl52rlzpzZu3Bi4f48cOVKSdPToUT344IO66qqr1KlTJ6Wnp6uwsFAffvhhyL6eeOIJFRQUqGPHjrrkkks0aNAgLVu2LGhNc7Niw4YNGjx4sCRp6tSpgTyRcjQlXnDEIgrceuutysvLU1lZmbZs2aLHH39c//rXv/S3v/0tsGb9+vV66aWXVFJSoq5duyovL081NTX67ne/GygemZmZeuuttzRt2jTV1dVpxowZkqSTJ09q9OjR+vLLL/XAAw8oJydHzz//vNavX9+ifNOmTdPSpUtVWFiou+++W2fOnNE777yjLVu2aNCgQXr++ed19913a8iQIbr33nslSX369JGkdssI4Bvbt2/X2LFjlZmZqTlz5ujMmTMqLS1VdnZ20Lry8nIVFBRo4sSJSkpK0uuvv6777rtPfr9fxcXFkqQFCxbo/vvvV6dOnfS73/1OkgLX8/nnn2vVqlX6yU9+ol69eqmmpkaLFy/WiBEjtGvXLuXk5Ej691OsDzzwgG655ZbAA6ePPvpI7733nm6//XZJLZsVV1xxhR5++GHNnj1b9957r6677jpJ0rXXXtsutyv+j0HEKi0tNZLMxIkTg7bfd999RpL58MMPjTHGSDIJCQlm586dQeumTZtmunfvbg4fPhy0/ac//anJyMgwJ06cMMYYs2DBAiPJvPTSS4E1DQ0Npm/fvkaSefvttwPbi4qKzKWXXho4v379eiPJPPDAAyH5/X5/4OfU1FRTVFQUsqYtMgI4v0mTJpmUlBSzb9++wLZdu3aZxMRE8+1/Fs7e/75t3Lhxpnfv3kHbCgoKzIgRI0LWnjp1yvh8vqBte/fuNW632zz88MOBbTfddJMpKCg4b+aWzorKykojySxZsuS814e2w1MhUeDsI4Oz7r//fknSm2++Gdg2YsQI5efnB84bY7RixQpNmDBBxhgdPnw4cBo3bpw8Ho+2bt0auJ7u3bvrlltuCVy+Y8eOgaML57NixQq5XC6VlpaG/K65t621V0YA3/D5fFq9erUmTZqknj17BrZfccUVGjduXNDaDh06BH72eDw6fPiwRowYoc8//1wej6fZfbndbiUkJAT2e+TIEXXq1EmXX3554L4tSRdffLH279+vysrKJq8nnFkB5/FUSBS47LLLgs736dNHCQkJ+uKLLwLbevXqFbTm0KFDOnbsmJ5++mk9/fTTTV5vbW2tJGnfvn3q27dvSBG4/PLLm8322WefKScnR507d27Jf4ojGQF849ChQzp58mTIXJH+fX/69gOWf/zjHyotLdXmzZt14sSJoLUej0cZGRnn3dfZ1189+eST2rt3r3w+X+B3Xbp0Cfz80EMPae3atRoyZIj69u2rsWPH6vbbb9fw4cMDmVs6K+A8ikUUaupIwLcfWUj/vkNL0s9+9jMVFRU1eT39+/e3Hy4M0ZARiFefffaZRo8erX79+mn+/PnKzc1VcnKy3nzzTf3lL38J3H/P59FHH9Xvf/973XXXXfrDH/6gzp07KyEhQTNmzAi6/BVXXKHdu3frjTfeUEVFhVasWKEnn3xSs2fP1ty5c5kVUYZiEQX27NkTdETi008/ld/vP+8nYGZmZiotLU0+n09jxow57/Vfeuml2rFjh4wxQaWlJZ9X0adPH61evVpHjx4971GLpspQe2UE8I3MzEx16NBBe/bsCfndt+9Pr7/+urxer1577bWgp0zefvvtkMud62nPV155RaNGjdKzzz4btP3YsWPq2rVr0LbU1FRNnjxZkydPVmNjo3784x/rkUce0axZs8KaFXxyqPN4jUUUWLRoUdD5J554QpJUWFh4zsskJibq5ptv1ooVK7Rjx46Q3x86dCjw8w9+8AMdOHBAr7zySmDbiRMnznnI8dtuvvlmGWM0d+7ckN8ZYwI/p6am6tixY45kBPCNxMREjRs3TqtWrdKXX34Z2P7xxx9r9erVQeuk4Puxx+PRkiVLQq6zqfv32ev49uUl6eWXX9ZXX30VtO3IkSNB55OTk5Wfny9jjE6fPh3WrEhNTZWkJvOgfXDEIgrs3btXEydO1I033qjNmzfr73//u26//XYNGDDgvJebN2+e3n77bQ0dOlT33HOP8vPzdfToUW3dulVr167V0aNHJUn33HOPFi5cqDvvvFNVVVXq3r27nn/+eXXs2LHZbKNGjdLPf/5zPf7449qzZ49uvPFG+f1+vfPOOxo1apRKSkokSQMHDtTatWs1f/585eTkqFevXho6dGi7ZAQQbO7cuaqoqNB1112n++67T2fOnAl8jsRHH30kSRo7dqySk5M1YcIE/eIXv9Dx48f1zDPPKCsrS19//XXQ9Q0cOFDl5eX64x//qL59+yorK0s33HCDfvjDH+rhhx/W1KlTde2112r79u164YUX1Lt376DLjx07Vt26ddPw4cOVnZ2tjz/+WAsXLtT48eOVlpYmqeXzrE+fPrr44ov11FNPKS0tTampqRo6dGjI69DQhhx5Lwpa5OzbTXft2mVuueUWk5aWZi655BJTUlJiTp48GVgnyRQXFzd5HTU1Naa4uNjk5uaaiy66yHTr1s2MHj3aPP3000Hr9u3bZyZOnGg6duxounbtaqZPn24qKiqafbupMcacOXPG/PnPfzb9+vUzycnJJjMz0xQWFpqqqqrAmk8++cRcf/31pkOHDkZS0FtPbWcE0LyNGzeagQMHmuTkZNO7d2/z1FNPBWbOWa+99prp37+/SUlJMXl5eeZPf/qTee6554wks3fv3sC6gwcPmvHjx5u0tDQjKfDW01OnTpnf/OY3pnv37qZDhw5m+PDhZvPmzWbEiBFBb09dvHixuf76602XLl2M2+02ffr0Mb/97W+Nx+MJytzSWfHqq6+a/Px8k5SUxFtPHeAy5j+OUyFizJkzR3PnztWhQ4dCno8EACAS8RoLAABgDcUCAABYQ7EAAADW8BoLAABgDUcsAACANRQLAABgTbt/QJbf79eBAweUlpbGR68CDjDGqL6+Xjk5OYFvnox0zA3AeS2dHe1eLA4cOKDc3Nz23i2A/1BdXa0ePXo4HaNFmBtA5GhudrR7sTj78ax9fzVbie6U9t79ebmP8jrWlsqs2Ot0hCb5ag81vyjOndFpvas3A/fFaHA2a+8Zs5UQYXMj8aTTCUKdSXU6QdN6/XWf0xGa5K8/7nSEEK7EyDuaeMY0aqPnxWZnR7sXi7OHMRPdKRFXLBKTKRYtlZSQ7HSEJrlcFzkdIfL93//m0fSUwtmsCZE4N5r/9vB2ZyLrJgpISnA7HaFJflej0xFCuFyRVyzOam52RG5yAAAQdSgWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArGlVsVi0aJHy8vKUkpKioUOH6v3337edC0AMYnYAsS/sYvHiiy9q5syZKi0t1datWzVgwACNGzdOtbW1bZEPQIxgdgDxIexiMX/+fN1zzz2aOnWq8vPz9dRTT6ljx4567rnn2iIfgBjB7ADiQ1jForGxUVVVVRozZsw3V5CQoDFjxmjz5s1NXsbr9aquri7oBCC+hDs7mBtA9AqrWBw+fFg+n0/Z2dlB27Ozs3Xw4MEmL1NWVqaMjIzAKTc3t/VpAUSlcGcHcwOIXm3+rpBZs2bJ4/EETtXV1W29SwBRjrkBRK+kcBZ37dpViYmJqqmpCdpeU1Ojbt26NXkZt9stt9vd+oQAol64s4O5AUSvsI5YJCcna+DAgVq3bl1gm9/v17p16zRs2DDr4QDEBmYHED/COmIhSTNnzlRRUZEGDRqkIUOGaMGCBWpoaNDUqVPbIh+AGMHsAOJD2MVi8uTJOnTokGbPnq2DBw/q6quvVkVFRciLsgDg25gdQHwIu1hIUklJiUpKSmxnARDjmB1A7OO7QgAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFjTqu8KscF91Cgx2Ti1+yZ5u7icjhDCfSSybqOzfDW1TkdoUmJ2ltMRQkTqbRWNEk9KiX6nUwQ7lRl599HEU5E3yyRJCZH5WNb07el0hBCu+hNORwjl80rHml8WmX9lAAAQlSgWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArKFYAAAAaygWAADAmrCLxaZNmzRhwgTl5OTI5XJp1apVbRALQCxhbgDxI+xi0dDQoAEDBmjRokVtkQdADGJuAPEjKdwLFBYWqrCwsC2yAIhRzA0gfoRdLMLl9Xrl9XoD5+vq6tp6lwCiHHMDiF5t/uLNsrIyZWRkBE65ubltvUsAUY65AUSvNi8Ws2bNksfjCZyqq6vbepcAohxzA4hebf5UiNvtltvtbuvdAIghzA0gevE5FgAAwJqwj1gcP35cn376aeD83r17tW3bNnXu3Fk9e/a0Gg5AbGBuAPEj7GLxwQcfaNSoUYHzM2fOlCQVFRVp6dKl1oIBiB3MDSB+hF0sRo4cKWNMW2QBEKOYG0D84DUWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArAn7u0JimftI5H2XgbeLy+kIUcVXU+t0hBCJ2VlORwhi/I1S5N1MLXImVTIpTqcIlngq8u6jjV19TkdokmlsdDpCk8z/7HQ6QghXj/9yOkIIl79lfz+OWAAAAGsoFgAAwBqKBQAAsIZiAQAArKFYAAAAaygWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwBqKBQAAsCasYlFWVqbBgwcrLS1NWVlZmjRpknbv3t1W2QDECGYHED/CKhYbN25UcXGxtmzZojVr1uj06dMaO3asGhoa2iofgBjA7ADiR1I4iysqKoLOL126VFlZWaqqqtL1119vNRiA2MHsAOJHWMXiP3k8HklS586dz7nG6/XK6/UGztfV1V3ILgHEgOZmB3MDiF6tfvGm3+/XjBkzNHz4cF155ZXnXFdWVqaMjIzAKTc3t7W7BBADWjI7mBtA9Gp1sSguLtaOHTu0fPny866bNWuWPB5P4FRdXd3aXQKIAS2ZHcwNIHq16qmQkpISvfHGG9q0aZN69Ohx3rVut1tut7tV4QDElpbODuYGEL3CKhbGGN1///1auXKlNmzYoF69erVVLgAxhNkBxI+wikVxcbGWLVumV199VWlpaTp48KAkKSMjQx06dGiTgACiH7MDiB9hvcaivLxcHo9HI0eOVPfu3QOnF198sa3yAYgBzA4gfoT9VAgAhIvZAcQPvisEAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANWF9V4hNmRV7lZSQ7NTum+SrqXU6QtRYfWCb0xGa9FpDR6cjNOF/nQ4Q5ES9T+v+n9MpWqfXX/cpKcHtdIxgCZH3+Mw0NjodoUlPV/630xGatLK+wOkIIS5zf+x0hBAn6n1ae3Xz6yLvHgEAAKIWxQIAAFhDsQAAANZQLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFgTVrEoLy9X//79lZ6ervT0dA0bNkxvvfVWW2UDECOYHUD8CKtY9OjRQ/PmzVNVVZU++OAD3XDDDbrpppu0c+fOtsoHIAYwO4D4kRTO4gkTJgSdf+SRR1ReXq4tW7aooKDAajAAsYPZAcSPsIrFt/l8Pr388stqaGjQsGHDzrnO6/XK6/UGztfV1bV2lwBiQEtmB3MDiF5hv3hz+/bt6tSpk9xut375y19q5cqVys/PP+f6srIyZWRkBE65ubkXFBhAdApndjA3gOgVdrG4/PLLtW3bNr333nv61a9+paKiIu3ateuc62fNmiWPxxM4VVdXX1BgANEpnNnB3ACiV9hPhSQnJ6tv376SpIEDB6qyslKPPfaYFi9e3OR6t9stt9t9YSkBRL1wZgdzA4heF/w5Fn6/P+i5UABoCWYHEJvCOmIxa9YsFRYWqmfPnqqvr9eyZcu0YcMGrV69uq3yAYgBzA4gfoRVLGpra3XnnXfq66+/VkZGhvr376/Vq1fr+9//flvlAxADmB1A/AirWDz77LNtlQNADGN2APGD7woBAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANRQLAABgDcUCAABYE9aXkNnkqz0kl+sip3bfpMTsLKcjhPDV1DodoUmvNXR0OkKTJqaecDpCiEi9raKRv/64/K5Gp2MEMX17Oh0hhPmfnU5HaNLK+gKnIzTpB50i7/Z652RvpyOEOOk706J1HLEAAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDUUCwAAYA3FAgAAWEOxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1F1Qs5s2bJ5fLpRkzZliKAyDWMTeA2NbqYlFZWanFixerf//+NvMAiGHMDSD2tapYHD9+XHfccYeeeeYZXXLJJbYzAYhBzA0gPrSqWBQXF2v8+PEaM2ZMs2u9Xq/q6uqCTgDiD3MDiA9J4V5g+fLl2rp1qyorK1u0vqysTHPnzg07GIDYwdwA4kdYRyyqq6s1ffp0vfDCC0pJSWnRZWbNmiWPxxM4VVdXtyoogOjE3ADiS1hHLKqqqlRbW6trrrkmsM3n82nTpk1auHChvF6vEhMTgy7jdrvldrvtpAUQdZgbQHwJq1iMHj1a27dvD9o2depU9evXTw899FDIcAAA5gYQX8IqFmlpabryyiuDtqWmpqpLly4h2wFAYm4A8YZP3gQAANaE/a6Q/7RhwwYLMQDEE+YGELs4YgEAAKyhWAAAAGsoFgAAwBqKBQAAsIZiAQAArKFYAAAAaygWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMCaC/6ukFjiq6l1OkKIxOwspyOcw/86HaBJrzV0dDpCiImpJ5yOEKTO73c6Qqu5EhPkckXW4yFXfWT9fSXJ1eO/nI7QpMvcHzsdoUnvnOztdIQQU9Ij79+jOpdfJS1YF1n3UAAAENUoFgAAwBqKBQAAsIZiAQAArKFYAAAAaygWAADAGooFAACwhmIBAACsoVgAAABrKBYAAMAaigUAALCGYgEAAKyhWAAAAGsoFgAAwJqwisWcOXPkcrmCTv369WurbABiBLMDiB9J4V6goKBAa9eu/eYKksK+CgBxiNkBxIew79lJSUnq1q1bW2QBEMOYHUB8CPs1Fnv27FFOTo569+6tO+64Q19++eV513u9XtXV1QWdAMSfcGYHcwOIXmEVi6FDh2rp0qWqqKhQeXm59u7dq+uuu0719fXnvExZWZkyMjICp9zc3AsODSC6hDs7mBtA9HIZY0xrL3zs2DFdeumlmj9/vqZNm9bkGq/XK6/XGzhfV1en3NxcjdRNSnJd1Npdx43E7CynIzTpl+++63SEqDEx9YTTEYLU1ft1yXc+l8fjUXp6uiMZmpsd55oboy/+uZJcye0ZtXldOzudIITrVKPTEZp0/4a1zS9ywMEzGU5HCDElvdbpCCFaOjsu6NVTF198sb7zne/o008/Pecat9stt9t9IbsBEGOamx3MDSB6XdDnWBw/flyfffaZunfvbisPgDjA7ABiV1jF4sEHH9TGjRv1xRdf6J///Kd+9KMfKTExUbfddltb5QMQA5gdQPwI66mQ/fv367bbbtORI0eUmZmp733ve9qyZYsyMzPbKh+AGMDsAOJHWMVi+fLlbZUDQAxjdgDxg+8KAQAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDUUCwAAYE1Y3xVigzFGknRGpyXT3nuPPsbf6HSEJp2o9zkdIWrU+f1ORwhSd/zfec7eF6NBYG6YCLw/+LxOJwjhYm6E5aTvjNMRQtS5ImtuSC2fHS7TztNl//79ys3Nbc9dAmhCdXW1evTo4XSMFmFuAJGjudnR7sXC7/frwIEDSktLk8vlavX11NXVKTc3V9XV1UpPT7eYMDZxe7VcrN9WxhjV19crJydHCQnR8Wyorbkhxf7f1yZuq5aLh9uqpbOj3Z8KSUhIsPooKT09PWb/iG2B26vlYvm2ysjIcDpCWGzPDSm2/762cVu1XKzfVi2ZHdHxcAUAAEQFigUAALAmaouF2+1WaWmp3G6301GiArdXy3FbxTb+vi3HbdVy3FbfaPcXbwIAgNgVtUcsAABA5KFYAAAAaygWAADAGooFAACwhmIBAACsidpisWjRIuXl5SklJUVDhw7V+++/73SkiFNWVqbBgwcrLS1NWVlZmjRpknbv3u10rKgwb948uVwuzZgxw+kosIi50TLMjtZjdkRpsXjxxRc1c+ZMlZaWauvWrRowYIDGjRun2tpap6NFlI0bN6q4uFhbtmzRmjVrdPr0aY0dO1YNDQ1OR4tolZWVWrx4sfr37+90FFjE3Gg5ZkfrMDv+j4lCQ4YMMcXFxYHzPp/P5OTkmLKyMgdTRb7a2lojyWzcuNHpKBGrvr7eXHbZZWbNmjVmxIgRZvr06U5HgiXMjdZjdjSP2fGNqDti0djYqKqqKo0ZMyawLSEhQWPGjNHmzZsdTBb5PB6PJKlz584OJ4lcxcXFGj9+fND/X4h+zI0Lw+xoHrPjG+3+7aYX6vDhw/L5fMrOzg7anp2drU8++cShVJHP7/drxowZGj58uK688kqn40Sk5cuXa+vWraqsrHQ6CixjbrQes6N5zI5gUVcs0DrFxcXasWOH3n33XaejRKTq6mpNnz5da9asUUpKitNxgIjB7Dg/ZkeoqCsWXbt2VWJiompqaoK219TUqFu3bg6limwlJSV64403tGnTJvXo0cPpOBGpqqpKtbW1uuaaawLbfD6fNm3apIULF8rr9SoxMdHBhLgQzI3WYXY0j9kRKupeY5GcnKyBAwdq3bp1gW1+v1/r1q3TsGHDHEwWeYwxKikp0cqVK7V+/Xr16tXL6UgRa/To0dq+fbu2bdsWOA0aNEh33HGHtm3bFneDIdYwN8LD7Gg5ZkeoqDtiIUkzZ85UUVGRBg0apCFDhmjBggVqaGjQ1KlTnY4WUYqLi7Vs2TK9+uqrSktL08GDByVJGRkZ6tChg8PpIktaWlrI88epqanq0qULzyvHCOZGyzE7Wo7ZESoqi8XkyZN16NAhzZ49WwcPHtTVV1+tioqKkBdmxbvy8nJJ0siRI4O2L1myRFOmTGn/QICDmBstx+zAhXAZY4zTIQAAQGyIutdYAACAyEWxAAAA1lAsAACANRQLAABgDcUCAABYQ7EAAADWUCwAAIA1FAsAAGANxQIAAFhDsQAAANZQLAAAgDX/H+H6W0PBqOauAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get the predicted G-matrix from the true model\n", "G,_ = Mtrue.predict([0,1,-1.5,-1])\n", "\n", "# The estimated G-matrix from the first dataset\n", "trial_type = D[1].obs_descriptors['cond_vec']\n", "G_hat,_ = pcm.est_G_crossval(D[1].measurements, trial_type, part_vec)\n", "\n", "# Visualize the second moment (G) matrices\n", "plt.subplot(1,2,1)\n", "plt.imshow(G)\n", "plt.title('predicted')\n", "\n", "plt.subplot(1,2,2)\n", "plt.imshow(G_hat)\n", "plt.title('dataset')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fitting models\n", "We now use a flexible correlation model (`Mflex`). In this case we do not want to model the condition as random effect, as we are removing them as a fixed effect when fitting the model." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The true model has 3 hyper parameters\n" ] } ], "source": [ "# Now make the flexible model without a condition effect (given that we are removing the condition effect as a fixed effect when fitting the model).\n", "Mflex = pcm.CorrelationModel(\"flex\", num_items=3, corr=None, cond_effect=False)\n", "print(f\"The true model has {Mflex.n_param} hyper parameters\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now fit the model to the datasets in one go. The resulting dataframe `T` has the log-likelihoods for each model (columns) / dataset (rows). The second return argument `theta` contains the parameters for each model fit. " ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "# If we want to remove the condition effect as a fixed effect, we need an indicator for that\n", "c_vec,item_vec= pcm.cond_to_item(cond_vec)\n", "X = pcm.matrix.indicator(c_vec)\n" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
variablelikelihoodnoiseiterations
modelflexflexflex
0-867.4813831.0136835.0
1-789.7482010.9148025.0
2-867.8862921.0243305.0
3-834.8341550.9870005.0
4-832.1244850.9824135.0
\n", "
" ], "text/plain": [ "variable likelihood noise iterations\n", "model flex flex flex\n", "0 -867.481383 1.013683 5.0\n", "1 -789.748201 0.914802 5.0\n", "2 -867.886292 1.024330 5.0\n", "3 -834.834155 0.987000 5.0\n", "4 -832.124485 0.982413 5.0" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "T, theta = pcm.fit_model_individ(D, Mflex, fixed_effect=X, fit_scale=False, verbose=False)\n", "T.head()" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
variablelikelihoodnoisescaleiterations
modelflexflexflexflex
0-870.6843701.0151430.05
1-790.2296140.9138110.05
2-868.2718461.0232430.05
3-836.5382270.9851930.05
4-833.7753460.9780060.05
\n", "
" ], "text/plain": [ "variable likelihood noise scale iterations\n", "model flex flex flex flex\n", "0 -870.684370 1.015143 0.0 5\n", "1 -790.229614 0.913811 0.0 5\n", "2 -868.271846 1.023243 0.0 5\n", "3 -836.538227 0.985193 0.0 5\n", "4 -833.775346 0.978006 0.0 5" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Now the group fit\n", "T_gr, theta_gr = pcm.fit_model_group(D, Mflex, fixed_effect=X, fit_scale=False, verbose=False)\n", "T_gr.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interpreting the model fits \n", "Again, the key diagnostic plot is to look at the ML estimates of the parameters against the overall SNR. Note that the first two parameters are the main condition effects " ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPI9JREFUeJzt3X90VPWd//HXZCAZccnwIyQzYL4QsGIjvwSbNIhHjwYT7MnC/nABRZAF3GaxR01tJa2QRq3xR0vRSknLguByKliLVQpNbaOhawlml8iRiKJgFNRMwo8yCXGT6Mz9/pHN4JDJjxuSzJ3J83HOHMmd971+Prn58cr9fO7n2gzDMAQAAGBhMeFuAAAAQFcILAAAwPIILAAAwPIILAAAwPIILAAAwPIILAAAwPIILAAAwPIILAAAwPIGhbsBvcHv9+uzzz7T0KFDZbPZwt0cAADQDYZhqKGhQaNHj1ZMTOfXUKIisHz22WdKTk4OdzMAAEAPnDhxQpdddlmnNVERWIYOHSqptcPx8fFhbg0AAOiO+vp6JScnB36PdyYqAkvbMFB8fDyBBQCACNOd6RxMugUAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJbXo8Cyfv16jRs3Tg6HQ+np6aqoqOi0ft26dZo4caIuueQSJScn67777lNTU1Pg/R/96Eey2WxBryuvvLInTQMAAFHI9MJxO3bsUF5enoqLi5Wenq5169YpKytLR44cUWJiYrv6X//611q1apU2b96smTNn6v3339edd94pm82mtWvXBuquuuoq/fnPfz7fsEFRsaYdAACW5/Mbqqg+o7qGJiUOdSgtZYTsMdZ6Np/pVLB27VqtWLFCS5culSQVFxdr9+7d2rx5s1atWtWuft++fbr22mt12223SZLGjRunhQsX6s033wxuyKBBcrlcPekDAADooZKqGhXuOqwa7/mRD7fToYKcVGVPcoexZcFMDQm1tLTowIEDyszMPH+AmBhlZmaqvLw85D4zZ87UgQMHAsNGH374ofbs2aNbbrklqO6DDz7Q6NGjNX78eN1+++06fvx4h+1obm5WfX190AsAAJhTUlWj3G2VQWFFkjzeJuVuq1RJVU2YWtaeqcBy6tQp+Xw+JSUlBW1PSkqSx+MJuc9tt92mhx56SLNmzdLgwYM1YcIE3XDDDfrBD34QqElPT9eWLVtUUlKiDRs2qLq6Wtddd50aGhpCHrOoqEhOpzPw4knNAACY4/MbKtx1WEaI99q2Fe46LJ8/VEX/6/O7hMrKyvToo4/qF7/4hSorK7Vz507t3r1bDz/8cKBmzpw5uvXWWzVlyhRlZWVpz549Onv2rF544YWQx8zPz5fX6w28Tpw40dfdAAAgqlRUn2l3ZeWrDEk13iZVVJ/pv0Z1wtQcloSEBNntdtXW1gZtr62t7XD+yerVq3XHHXdo+fLlkqTJkyersbFRd911l374wx8qJqZ9Zho2bJiuuOIKHT16NOQx4+LiFBcXZ6bpAADgK+oaOg4rPanra6ausMTGxmrGjBkqLS0NbPP7/SotLVVGRkbIfT7//PN2ocRut0uSDCP0ZaZz587p2LFjcrutM9kHAIBokjjU0at1fc30XUJ5eXlasmSJrrnmGqWlpWndunVqbGwM3DW0ePFijRkzRkVFRZKknJwcrV27VldffbXS09N19OhRrV69Wjk5OYHgcv/99ysnJ0djx47VZ599poKCAtntdi1cuLAXuwoAANqkpYyQ2+mQx9sUch6LTZLL2XqLsxWYDizz58/XyZMntWbNGnk8Hk2bNk0lJSWBibjHjx8PuqLy4IMPymaz6cEHH9Snn36qUaNGKScnRz/+8Y8DNZ988okWLlyo06dPa9SoUZo1a5b279+vUaNG9UIXAQDAhewxNhXkpCp3W6VsUlBoaVuBpSAn1TLrsdiMjsZlIkh9fb2cTqe8Xq/i4+PD3RwAACJGONdhMfP7m+VkAQAYwLInuTU71RV9K90CAIDoYo+xKWPCyHA3o1M8rRkAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFgegQUAAFjeoHA3AAAAWIPPb6ii+ozqGpqUONShtJQRssfYwt0sSQQWAAAgqaSqRoW7DqvG2xTY5nY6VJCTquxJ7jC2rBVDQgAADHAlVTXK3VYZFFYkyeNtUu62SpVU1YSpZef1KLCsX79e48aNk8PhUHp6uioqKjqtX7dunSZOnKhLLrlEycnJuu+++9TUFPxJMXtMAABw8Xx+Q4W7DssI8V7btsJdh+Xzh6roP6YDy44dO5SXl6eCggJVVlZq6tSpysrKUl1dXcj6X//611q1apUKCgr07rvvatOmTdqxY4d+8IMf9PiYAACgd1RUn2l3ZeWrDEk13iZVVJ/pv0aFYDqwrF27VitWrNDSpUuVmpqq4uJiDRkyRJs3bw5Zv2/fPl177bW67bbbNG7cON18881auHBh0BUUs8cEAAC9o66h47DSk7q+YiqwtLS06MCBA8rMzDx/gJgYZWZmqry8POQ+M2fO1IEDBwIB5cMPP9SePXt0yy239PiYzc3Nqq+vD3oBAADzEoc6erWur5i6S+jUqVPy+XxKSkoK2p6UlKT33nsv5D633XabTp06pVmzZskwDH355Zf69re/HRgS6skxi4qKVFhYaKbpAAAghLSUEXI7HfJ4m0LOY7FJcjlbb3EOpz6/S6isrEyPPvqofvGLX6iyslI7d+7U7t279fDDD/f4mPn5+fJ6vYHXiRMnerHFAAAMHPYYmwpyUiW1hpOvavu4ICc17OuxmLrCkpCQILvdrtra2qDttbW1crlcIfdZvXq17rjjDi1fvlySNHnyZDU2Nuquu+7SD3/4wx4dMy4uTnFxcWaaDgAAOpA9ya0Ni6a3W4fFZaF1WEwFltjYWM2YMUOlpaWaN2+eJMnv96u0tFR33313yH0+//xzxcQEX8ix2+2SJMMwenRMAADQu7InuTU71RU9K93m5eVpyZIluuaaa5SWlqZ169apsbFRS5culSQtXrxYY8aMUVFRkSQpJydHa9eu1dVXX6309HQdPXpUq1evVk5OTiC4dHVMAADQ9+wxNmVMGBnuZoRkOrDMnz9fJ0+e1Jo1a+TxeDRt2jSVlJQEJs0eP3486IrKgw8+KJvNpgcffFCffvqpRo0apZycHP34xz/u9jEBAMDAZjMMI7xL1/WC+vp6OZ1Oeb1excfHh7s5AACgG8z8/ubhhwCikpWfOgvAPAILgKhj9afOAjCPpzUDiCqR8NRZAOYRWABEjUh56iwA8wgsAKJGpDx1FoB5BBYAUSNSnjoLwDwCC4CoESlPnQVgHoEFQNRoe+psRzcv29R6t1C4nzoLwDwCC4CoESlPnQVgHoEFQFRpe+qsyxk87ONyOrRh0XTWYQEiFAvHAYg6Vn/qLADzCCwAopKVnzoLwDyGhAAAgOURWAAAgOURWAAAgOUxhwWwAJ/fYIIoAHSCwAKEWUlVjQp3HQ56Bo7b6VBBTiq34ALA/2FICAijkqoa5W6rbPfAPo+3SbnbKlVSVROmlgGAtRBYgDDx+Q0V7josI8R7bdsKdx2Wzx+qAgAGFgILECYV1WfaXVn5KkNSjbdJFdVn+q9RAGBRBBYgTOoaOg4rPakDgGhGYAHCJHGoo+siE3UAEM0ILECYpKWMkNvpaPdU4TY2td4tlJYyoj+bBQCWRGABwsQeY1NBTqoktQstbR8X5KSyHgsAiMAChFX2JLc2LJoulzN42MfldGjDoumswwIA/4eF44Awy57k1uxUFyvdAkAnCCyABdhjbMqYMDLczQAAy2JICAAAWB6BBQAAWB6BBQAAWB5zWAAMOD6/wSRnIML06ArL+vXrNW7cODkcDqWnp6uioqLD2htuuEE2m63d61vf+lag5s4772z3fnZ2dk+aBgCdKqmq0azHX9PCjft1z/aDWrhxv2Y9/hpPxgYsznRg2bFjh/Ly8lRQUKDKykpNnTpVWVlZqqurC1m/c+dO1dTUBF5VVVWy2+269dZbg+qys7OD6p5//vme9QgAOlBSVaPcbZXtHjrp8TYpd1sloQWwMNOBZe3atVqxYoWWLl2q1NRUFRcXa8iQIdq8eXPI+hEjRsjlcgVef/rTnzRkyJB2gSUuLi6obvjw4T3rEQCE4PMbKtx1WEaI99q2Fe46LJ8/VAWAcDMVWFpaWnTgwAFlZmaeP0BMjDIzM1VeXt6tY2zatEkLFizQpZdeGrS9rKxMiYmJmjhxonJzc3X69OkOj9Hc3Kz6+vqgFwB0pqL6TLsrK19lSKrxNqmi+kz/NQpAt5kKLKdOnZLP51NSUlLQ9qSkJHk8ni73r6ioUFVVlZYvXx60PTs7W88995xKS0v1+OOPa+/evZozZ458Pl/I4xQVFcnpdAZeycnJZroBYACqa+g4rPSkLhx8fkPlx07r5YOfqvzYaa4GYUDp17uENm3apMmTJystLS1o+4IFCwL/njx5sqZMmaIJEyaorKxMN910U7vj5OfnKy8vL/BxfX09oQVApxKHOrouMlHX30qqalS463DQVSK306GCnFSeOYUBwdQVloSEBNntdtXW1gZtr62tlcvl6nTfxsZGbd++XcuWLevy/zN+/HglJCTo6NGjId+Pi4tTfHx80AsAOpOWMkJup6Pdk7Hb2NQaANJSRvRns7qFycKAycASGxurGTNmqLS0NLDN7/ertLRUGRkZne77m9/8Rs3NzVq0aFGX/59PPvlEp0+fltvNXw0Aeoc9xqaCnFRJahda2j4uyEm13HosTBYGWpm+SygvL08bN27U1q1b9e677yo3N1eNjY1aunSpJGnx4sXKz89vt9+mTZs0b948jRwZ/IC3c+fO6Xvf+57279+vjz76SKWlpZo7d64uv/xyZWVl9bBbANBe9iS3NiyaLpczeNjH5XRow6LplhxaYbIw0Mr0HJb58+fr5MmTWrNmjTwej6ZNm6aSkpLARNzjx48rJiY4Bx05ckRvvPGGXn311XbHs9vtevvtt7V161adPXtWo0eP1s0336yHH35YcXFxPewWAISWPcmt2amuiFnpNhomCwO9wWYYRsRfR6yvr5fT6ZTX62U+C4CoUn7stBZu3N9l3fMrvqmMCSO7rAOsxMzvbx5+CAAWFsmThYHeRGABAAuL1MnCQG8jsACAxUXiZGGgt/XrwnEAgJ6JtMnCQG8jsABAhLDH2JhYiwGLISEAAGB5XGEBeonPb3C5HgD6CIEF6AU8mA4A+hZDQsBF4sF0AND3CCzAReDBdADQPwgswEXgwXQA0D8ILMBF4MF0ANA/CCzARUgc6ui6yEQdACA0AgtwEXgwHQD0DwILcBF4MB0A9A8CC3CReDAdAPQ9Fo4DegEPpkM0YdVmWBGBBeglPJgO0YBVm2FVDAkBACSxajOsjcACAGDVZlgegQUAwKrNsDzmsACwNCaA9g8rrNrMuUZnCCwALIsJoP0n3Ks2c67RFYaEAFgSE0D7VzhXbeZcozsILAAshwmg/S9cqzZzrtFdBBYAlsME0PAIx6rNnGt0F3NYAFiOFSaADlT9vWoz5xrdRWABYDnhngA60PXnqs2ca3QXQ0IALCecE0DRvzjX6C4CCwDLCdcEUPQ/zjW6i8ACwJLCMQEU4cG5RnfYDMMwfa/Y+vXr9eSTT8rj8Wjq1Kn6+c9/rrS0tJC1N9xwg/bu3dtu+y233KLdu3dLkgzDUEFBgTZu3KizZ8/q2muv1YYNG/S1r32tW+2pr6+X0+mU1+tVfHy82e4AsDBWPx04ONcDj5nf36Yn3e7YsUN5eXkqLi5Wenq61q1bp6ysLB05ckSJiYnt6nfu3KmWlpbAx6dPn9bUqVN16623BrY98cQTevrpp7V161alpKRo9erVysrK0uHDh+VwMNEKGMj6cwIowotzjc6YvsKSnp6ub3zjG3rmmWckSX6/X8nJyfrOd76jVatWdbn/unXrtGbNGtXU1OjSSy+VYRgaPXq0vvvd7+r++++XJHm9XiUlJWnLli1asGBBl8fkCgsAAJHHzO9vU3NYWlpadODAAWVmZp4/QEyMMjMzVV5e3q1jbNq0SQsWLNCll14qSaqurpbH4wk6ptPpVHp6eofHbG5uVn19fdALAABEL1OB5dSpU/L5fEpKSgranpSUJI/H0+X+FRUVqqqq0vLlywPb2vYzc8yioiI5nc7AKzk52Uw3AABAhOnXu4Q2bdqkyZMndzhBt7vy8/Pl9XoDrxMnTvRSCwEAgBWZCiwJCQmy2+2qra0N2l5bWyuXy9Xpvo2Njdq+fbuWLVsWtL1tPzPHjIuLU3x8fNALQHTz+Q2VHzutlw9+qvJjp3kYHjDAmAossbGxmjFjhkpLSwPb/H6/SktLlZGR0em+v/nNb9Tc3KxFixYFbU9JSZHL5Qo6Zn19vd58880ujwlgYCipqtGsx1/Two37dc/2g1q4cb9mPf6aSqpqwt00AP3E9JBQXl6eNm7cqK1bt+rdd99Vbm6uGhsbtXTpUknS4sWLlZ+f326/TZs2ad68eRo5MviWNZvNpnvvvVePPPKIXnnlFR06dEiLFy/W6NGjNW/evJ71CkDUKKmqUe62ynZP9PV4m5S7rZLQAgwQptdhmT9/vk6ePKk1a9bI4/Fo2rRpKikpCUyaPX78uGJignPQkSNH9MYbb+jVV18Neczvf//7amxs1F133aWzZ89q1qxZKikpYQ0WYIDz+Q0V7jqsUIM/hlqXbi/cdVizU10sMIaIxqJ5XevRSrdWwzosQHQqP3ZaCzfu77Lu+RXfZMExRKySqhoV7jocdBXR7XSoICc16h9L0GfrsABAf6praOq6yEQdYDUMeXYfgQWAZSUO7d6wcHfrACvpashTah3y5I64VgQWAJaVljJCbqdDHY3k29R66TwtZUR/NgvoFRXVZ9pdWfkqQ1KNt0kV1Wf6r1EWRmABEMRK653YY2wqyEmVpHahpe3jgpxUJiciIjHkaY7pu4QARC8rTv7LnuTWhkXT27XLNUAmJSJ6MeRpDoEFgKTzk/8uvJ7SNvlvw6LpYQ0ts1Nd3PaJqNI25OnxNoWcx2JTazBnyLMVQ0IAImLynz3GpowJIzV32hhlTBhJWEHEY8jTHAILACb/AWHSNuTpcgYP+7icjrBe1bQihoQAMPkPCCOGPLuHwAKAyX9AmLUNeaJjDAkBYL0TAJZHYAHA5D8AlkdgASCJyX8ArI05LAACmPwHwKoILACCMPkPgBUxJAQAACyPwAIAACyPwAIAACyPwAIAACyPwAIAACyPwAIAACyPwAIAACyPwAIAACyPheMAYIDx+Q1WM0bEIbAAwABSUlWjwl2HVeNtCmxzOx0qyEnleVGwNIaEAGCAKKmqUe62yqCwIkkeb5Nyt1WqpKomTC0DukZgAYABwOc3VLjrsIwQ77VtK9x1WD5/qAog/AgsADAAVFSfaXdl5asMSTXeJlVUn+m/RgEmEFgAYACoa+g4rPSkDuhvBBYAGAAShzp6tQ7obwQWABgA0lJGyO10qKObl21qvVsoLWVEfzYL6LYeBZb169dr3LhxcjgcSk9PV0VFRaf1Z8+e1cqVK+V2uxUXF6crrrhCe/bsCbz/ox/9SDabLeh15ZVX9qRpAIAQ7DE2FeSkSlK70NL2cUFOKuuxwLJMB5YdO3YoLy9PBQUFqqys1NSpU5WVlaW6urqQ9S0tLZo9e7Y++ugjvfjiizpy5Ig2btyoMWPGBNVdddVVqqmpCbzeeOONnvUIABBS9iS3NiyaLpczeNjH5XRow6LprMMCSzO9cNzatWu1YsUKLV26VJJUXFys3bt3a/PmzVq1alW7+s2bN+vMmTPat2+fBg8eLEkaN25c+4YMGiSXy2W2OQAAE7InuTU71cVKt4g4pq6wtLS06MCBA8rMzDx/gJgYZWZmqry8POQ+r7zyijIyMrRy5UolJSVp0qRJevTRR+Xz+YLqPvjgA40ePVrjx4/X7bffruPHj3fYjubmZtXX1we9AADdY4+xKWPCSM2dNkYZE0YSVhARTAWWU6dOyefzKSkpKWh7UlKSPB5PyH0+/PBDvfjii/L5fNqzZ49Wr16tn/70p3rkkUcCNenp6dqyZYtKSkq0YcMGVVdX67rrrlNDQ0PIYxYVFcnpdAZeycnJZroBAAAiTJ8/S8jv9ysxMVG/+tWvZLfbNWPGDH366ad68sknVVBQIEmaM2dOoH7KlClKT0/X2LFj9cILL2jZsmXtjpmfn6+8vLzAx/X19YQWAACimKnAkpCQILvdrtra2qDttbW1Hc4/cbvdGjx4sOx2e2Db17/+dXk8HrW0tCg2NrbdPsOGDdMVV1yho0ePhjxmXFyc4uLi2r/R2Ch95f8TYLdLDkdwXUdiYqRLLulZ7eefS0YHy1rbbNKQIT2r/d//lfz+jttx6aU9q21qki4Ymutx7ZAhre2WpOZm6csve6f2kktaP8+S1NIiffFF79Q6HOe/VszUfvFFa31H4uKkQYPM1375ZevnoiOxsdL/zQEzVevztZ67jgwe3Fpvttbvb/1a643aQYNaPxdS6/fE55/3Tq2Z73t+RoSu5WeE+Vp+RrT+28zPiO4yTEpLSzPuvvvuwMc+n88YM2aMUVRUFLI+Pz/fGDt2rOHz+QLb1q1bZ7jd7g7/Hw0NDcbw4cONp556qltt8nq9hiTD2/rt3f51yy3BOwwZErpOMozrrw+uTUjouPaaa4Jrx47tuDY1Nbg2NbXj2rFjg2uvuabj2oSE4Nrrr++4dsiQ4Npbbum49sIvjX/+585rz507X7tkSee1dXXna//93zuvra4+X3v//Z3XVlWdry0o6Ly2ouJ87RNPdF77+uvna595pvPa3//+fO2zz3Ze+8IL52tfeKHz2mefPV/7+993XvvMM+drX3+989onnjhfW1HReW1BwfnaqqrOa++//3xtdXXntf/+7+dr6+o6r12y5HztuXOd1/7zPxtBOqvlZ0Tri58R518VFcaXPr+x7+gpo+reH3Zey8+I1lcPfkYEfn97vUZXTN/WnJeXp40bN2rr1q169913lZubq8bGxsBdQ4sXL1Z+fn6gPjc3V2fOnNE999yj999/X7t379ajjz6qlStXBmruv/9+7d27Vx999JH27dunf/iHf5DdbtfChQvNNg8AgIu279gpzXr8NS3cuF8vH/ws3M2BJJthGIbZnZ555hk9+eST8ng8mjZtmp5++mmlp6dLkm644QaNGzdOW7ZsCdSXl5frvvvu08GDBzVmzBgtW7ZMDzzwQGCYaMGCBfrLX/6i06dPa9SoUZo1a5Z+/OMfa8KECd1qT319vZxOp7yffab4+Pj2BVzuDV3L5V7ztVzubf03Q0I9q+VnROu/Lf4z4k/vePTtF9+RL6b1+36w7wsN8vkCC+w9tWCaZl/1lWkQUf4zwjdocOtt8Gcb5YqVrhnXwW3wPfgZEfj97fWG/v39FT0KLFZjpsMDkc9vsOYCAHSDz29o1uOvdfhka5taF9p744EbB8TP0ZKqGhXuOhz0+XA7HSrISe2VhQbN/P7u87uEEF59/cWG0AiJQGSqqD7TYViRJENSjbdJFdVnlDFhZP81LAxKqmqUu61SF17V8HiblLutst9XRyawRDGrfbENFIREIHLVNXQyPNKDukjl8xsq3HW43e8PqTW02SQV7jqs2amufvtjjKc1R6muvtik1i82nz/iRwQtpS0kXvgXWltILKmqCVPLAHRH4lBH10Um6iKVmStN/YXAEqWs+MUW7QiJQORLSxkht9PR7onWbWxqvWKaljKiP5vV76x4pYnAEqWs+MUW7QiJQOSzx9hUkJMqSe1CS9vHBTmpUT8nzYpXmggsUcqKX2zRjpAIRIfsSW5tWDRdLmfwz0eX0zFg5v5Z8UoTk26jVNsXm8fbFHKIou3WvGi/rNmfCIlA9Mie5NbsVNeAvduv7UpT7rZK2aSg3yPhutLEFZYoxWXN/mfFv0gA9Jw9xqaMCSM1d9oYZUwYOeB+XlrtShMLx0U5brHtX213CUmh/yIZKJeTAUSPvlxXipVuEYRFzPoXIREAuofAAoQZIREAusbS/ECYtY19AwB6B5NuAQCA5XGFBYApDHcBCAcCC4BuY0IxgHBhSAhAt/BgRwDhRGAB0CUe7Agg3AgsALrEgx0BhBuBBUCXeLAjgHAjsADoEg92BBBuBBYAXeLBjgDCjcACoEs8/RtAuBFYAHSL1R41D2BgYeE4AN2WPcmt2akuVroF0O8ILABM4cGOAMKBISEAAGB5BBYAAGB5BBYAAGB5BBYAAGB5BBYAAGB5BBYAAGB5PQos69ev17hx4+RwOJSenq6KiopO68+ePauVK1fK7XYrLi5OV1xxhfbs2XNRxwQQfXx+Q+XHTuvlg5+q/Nhp+fxGuJsEwCJMr8OyY8cO5eXlqbi4WOnp6Vq3bp2ysrJ05MgRJSYmtqtvaWnR7NmzlZiYqBdffFFjxozRxx9/rGHDhvX4mACiT0lVjQp3HVaN9/wTn91OhwpyUllFF4BshmGY+hMmPT1d3/jGN/TMM89Ikvx+v5KTk/Wd73xHq1ataldfXFysJ598Uu+9954GDx7cK8e8UH19vZxOp7xer+Lj4810B4AFlFTVKHdbpS78YdS2fi5L/wPRyczvb1NDQi0tLTpw4IAyMzPPHyAmRpmZmSovLw+5zyuvvKKMjAytXLlSSUlJmjRpkh599FH5fL4eH7O5uVn19fVBLwCRyec3VLjrcLuwIimwrXDXYYaHgAHOVGA5deqUfD6fkpKSgrYnJSXJ4/GE3OfDDz/Uiy++KJ/Ppz179mj16tX66U9/qkceeaTHxywqKpLT6Qy8kpOTzXQDgIVUVJ8JGga6kCGpxtukiuoz/dcoAJbT53cJ+f1+JSYm6le/+pVmzJih+fPn64c//KGKi4t7fMz8/Hx5vd7A68SJE73YYgD9qa6h47DSkzoA0cnUpNuEhATZ7XbV1tYGba+trZXL5Qq5j9vt1uDBg2W32wPbvv71r8vj8ailpaVHx4yLi1NcXJyZpgOwqMShjl6tAxCdTF1hiY2N1YwZM1RaWhrY5vf7VVpaqoyMjJD7XHvttTp69Kj8fn9g2/vvvy+3263Y2NgeHRNA9EhLGSG30xGYYHshm1rvFkpLGdGfzQJgMaaHhPLy8rRx40Zt3bpV7777rnJzc9XY2KilS5dKkhYvXqz8/PxAfW5urs6cOaN77rlH77//vnbv3q1HH31UK1eu7PYxAUQve4xNBTmpktQutLR9XJCTKntMR5EGwEBgeh2W+fPn6+TJk1qzZo08Ho+mTZumkpKSwKTZ48ePKybmfA5KTk7WH//4R913332aMmWKxowZo3vuuUcPPPBAt48JILplT3Jrw6Lp7dZhcbEOC4D/Y3odFitiHRYgOvj8hiqqz6iuoUmJQ1uHgbiy0j187hCJzPz+Nn2FBQD6ij3GpowJI8PdjIjDKsEYCHj4IQBEsLZVgi9cy8bjbVLutkqVVNWEqWVA7yKwAECEYpVgDCQEFgCIUKwSjIGEOSywBCYMAuaxSjAGEgILwo4Jg0DPsEowBhKGhBBWTBgEeo5VgjGQEFgQNkwYBC4OqwRjICGwIGyYMAhcvLZVgl3O4GEfl9OhDYumM6yKqMEcFoQNEwaB3pE9ya3ZqS4mriOqEVgQNkwYBHoPqwQj2jEkhLBhwiAAoLsILAgbJgwCALqLwIKwYsIgAKA7mMOCsGPCIACgKwQWWAITBgEAnWFICAAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB53CQFh4PMb3MYNACYQWIB+VlJVo8Jdh4OeVO12OlSQk8pCeQDQAYaEgH5UUlWj3G2VQWFFkjzeJuVuq1RJVU2YWgYA1kZgAfqJz2+ocNdhGSHea9tWuOuwfP5QFUDP+fyGyo+d1ssHP1X5sdN8jSEiMSQE9JOK6jPtrqx8lSGpxtukiuozrPqLXsMQJKIFV1iAflLX0HFY6Ukd0BWGIBFNuMIC9JPEoY6ui0zU9SXuYrKu7p6broYgbWodgpyd6uLcIiIQWIB+kpYyQm6nQx5vU8hfIjZJLmfrL6BwYgjBusycG4YgEW0YEgL6iT3GpoKcVEmt4eSr2j4uyEkN61+7DCFYl9lzwxAkog2BBehH2ZPc2rBoulzO4GEfl9OhDYumh/UKRiTdxTTQ7nrpybmJpCFIoDt6FFjWr1+vcePGyeFwKD09XRUVFR3WbtmyRTabLejlcAR/g9x5553tarKzs3vSNMDysie59cYDN+r5Fd/UUwum6fkV39QbD9wY9uEWM0MI4VRSVaNZj7+mhRv3657tB7Vw437Nevy1qL7605Nz0zYE2dH1Optah5PCPQQJdJfpwLJjxw7l5eWpoKBAlZWVmjp1qrKyslRXV9fhPvHx8aqpqQm8Pv7443Y12dnZQTXPP/+82aYBEcMeY1PGhJGaO22MMiaMtMSkx0gYQhioQ1Y9OTeRMAQJmGE6sKxdu1YrVqzQ0qVLlZqaquLiYg0ZMkSbN2/ucB+bzSaXyxV4JSUltauJi4sLqhk+fLjZpgG4CFYfQoikIave1tNzY+UhSMAsU3cJtbS06MCBA8rPzw9si4mJUWZmpsrLyzvc79y5cxo7dqz8fr+mT5+uRx99VFdddVVQTVlZmRITEzV8+HDdeOONeuSRRzRyZOiZ683NzWpubg58XF9fb6YbAEKw+l1MA/mul4s5N9mT3Jqd6uI2dUQ8U1dYTp06JZ/P1+4KSVJSkjweT8h9Jk6cqM2bN+vll1/Wtm3b5Pf7NXPmTH3yySeBmuzsbD333HMqLS3V448/rr1792rOnDny+Xwhj1lUVCSn0xl4JScnm+kGgBCsPoQQCUNWfeViz40VhyABs/r8LqGMjAwtXrxY06ZN0/XXX6+dO3dq1KhR+uUvfxmoWbBggf7+7/9ekydP1rx58/T73/9e//3f/62ysrKQx8zPz5fX6w28Tpw40dfdAAYEKw8hWH3Iqq9Z+dwA/cHUkFBCQoLsdrtqa2uDttfW1srlcnXrGIMHD9bVV1+to0ePdlgzfvx4JSQk6OjRo7rpppvavR8XF6e4uDgzTQfQTVYdQrD6kFV/sOq5AfqDqSsssbGxmjFjhkpLSwPb/H6/SktLlZGR0a1j+Hw+HTp0SG53x38NfPLJJzp9+nSnNQD6jhWHEKw+ZNVfrHhugP5gekgoLy9PGzdu1NatW/Xuu+8qNzdXjY2NWrp0qSRp8eLFQZNyH3roIb366qv68MMPVVlZqUWLFunjjz/W8uXLJbVOyP3e976n/fv366OPPlJpaanmzp2ryy+/XFlZWb3UTQDRgGERYOAy/Syh+fPn6+TJk1qzZo08Ho+mTZumkpKSwETc48ePKybmfA7629/+phUrVsjj8Wj48OGaMWOG9u3bp9TU1r+U7Ha73n77bW3dulVnz57V6NGjdfPNN+vhhx9m2AdAOwyLAAOTzTCMiF+0oL6+Xk6nU16vV/Hx8eFuDgAA6AYzv795lhAAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALA8AgsAALC8HgWW9evXa9y4cXI4HEpPT1dFRUWHtVu2bJHNZgt6ORyOoBrDMLRmzRq53W5dcsklyszM1AcffNCTpgEAgChkOrDs2LFDeXl5KigoUGVlpaZOnaqsrCzV1dV1uE98fLxqamoCr48//jjo/SeeeEJPP/20iouL9eabb+rSSy9VVlaWmpqazPcIAABEHdOBZe3atVqxYoWWLl2q1NRUFRcXa8iQIdq8eXOH+9hsNrlcrsArKSkp8J5hGFq3bp0efPBBzZ07V1OmTNFzzz2nzz77TL/73e961CkAABBdTAWWlpYWHThwQJmZmecPEBOjzMxMlZeXd7jfuXPnNHbsWCUnJ2vu3Ll65513Au9VV1fL4/EEHdPpdCo9Pb3DYzY3N6u+vj7oBQAX8vkNlR87rZcPfqryY6fl8xvhbhKAHhpkpvjUqVPy+XxBV0gkKSkpSe+9917IfSZOnKjNmzdrypQp8nq9+slPfqKZM2fqnXfe0WWXXSaPxxM4xoXHbHvvQkVFRSosLDTTdAADTElVjQp3HVaN9/zQstvpUEFOqrInucPYMgA90ed3CWVkZGjx4sWaNm2arr/+eu3cuVOjRo3SL3/5yx4fMz8/X16vN/A6ceJEL7YYQKQrqapR7rbKoLAiSR5vk3K3VaqkqiZMLQPQU6YCS0JCgux2u2pra4O219bWyuVydesYgwcP1tVXX62jR49KUmA/M8eMi4tTfHx80AsApNZhoMJdhxVq8KdtW+GuwwwPARHGVGCJjY3VjBkzVFpaGtjm9/tVWlqqjIyMbh3D5/Pp0KFDcrtbL8mmpKTI5XIFHbO+vl5vvvlmt48JAG0qqs+0u7LyVYakGm+TKqrP9F+jAFw0U3NYJCkvL09LlizRNddco7S0NK1bt06NjY1aunSpJGnx4sUaM2aMioqKJEkPPfSQvvnNb+ryyy/X2bNn9eSTT+rjjz/W8uXLJbXeQXTvvffqkUce0de+9jWlpKRo9erVGj16tObNm9d7PQUwINQ1dG85hO7WAbAG04Fl/vz5OnnypNasWSOPx6Np06appKQkMGn2+PHjiok5f+Hmb3/7m1asWCGPx6Phw4drxowZ2rdvn1JTUwM13//+99XY2Ki77rpLZ8+e1axZs1RSUtJugTkA6Eri0O793OhuHQBrsBmGEfEDufX19XI6nfJ6vcxnAQY4n9/QrMdfk8fbFHIei02Sy+nQGw/cKHuMrb+bB+ArzPz+5llCAKKKPcamgpzWK7gXxpG2jwtyUgkrQIQhsACIOtmT3NqwaLpczuBhH5fToQ2LprMOCxCBTM9hAYBIkD3JrdmpLlVUn1FdQ5MShzqUljKCKytAhCKwAIha9hibMiaMDHczAPQChoQAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlEVgAAIDlDQp3A6zM5zdUUX1GdQ1NShzqUFrKCNljbOFuFgAAAw6BpQMlVTUq3HVYNd6mwDa306GCnFRlT3KHsWUAAAw8DAmFUFJVo9xtlUFhRZI83iblbqtUSVVNmFoGAMDARGC5gM9vqHDXYRkh3mvbVrjrsHz+UBUAAKAvEFguUFF9pt2Vla8yJNV4m1RRfab/GgUAwABHYLlAXUPHYaUndQAA4OIRWC6QONTRq3UAAODi9SiwrF+/XuPGjZPD4VB6eroqKiq6td/27dtls9k0b968oO133nmnbDZb0Cs7O7snTbtoaSkj5HY61NHNyza13i2UljKiP5sFAMCAZjqw7NixQ3l5eSooKFBlZaWmTp2qrKws1dXVdbrfRx99pPvvv1/XXXddyPezs7NVU1MTeD3//PNmm9Yr7DE2FeSkSlK70NL2cUFOKuuxAADQj0wHlrVr12rFihVaunSpUlNTVVxcrCFDhmjz5s0d7uPz+XT77bersLBQ48ePD1kTFxcnl8sVeA0fPtxs03pN9iS3NiyaLpczeNjH5XRow6LpUbcOi89vqPzYab188FOVHzvNHVAAAMsxtXBcS0uLDhw4oPz8/MC2mJgYZWZmqry8vMP9HnroISUmJmrZsmX6r//6r5A1ZWVlSkxM1PDhw3XjjTfqkUce0ciRI800r1dlT3Jrdqor6le6ZYE8AEAkMBVYTp06JZ/Pp6SkpKDtSUlJeu+990Lu88Ybb2jTpk06ePBgh8fNzs7WP/7jPyolJUXHjh3TD37wA82ZM0fl5eWy2+3t6pubm9Xc3Bz4uL6+3kw3us0eY1PGhPCFpr7WtkDehddT2hbIi8arSQCAyNSnS/M3NDTojjvu0MaNG5WQkNBh3YIFCwL/njx5sqZMmaIJEyaorKxMN910U7v6oqIiFRYW9kmbB4quFsizqXWBvNmprqi7qgQAiDym5rAkJCTIbrertrY2aHttba1cLle7+mPHjumjjz5STk6OBg0apEGDBum5557TK6+8okGDBunYsWMh/z/jx49XQkKCjh49GvL9/Px8eb3ewOvEiRNmugGxQB4AILKYusISGxurGTNmqLS0NHBrst/vV2lpqe6+++529VdeeaUOHToUtO3BBx9UQ0ODnnrqKSUnJ4f8/3zyySc6ffq03O7QwxFxcXGKi4sz03RcgAXyAACRxPSQUF5enpYsWaJrrrlGaWlpWrdunRobG7V06VJJ0uLFizVmzBgVFRXJ4XBo0qRJQfsPGzZMkgLbz507p8LCQv3TP/2TXC6Xjh07pu9///u6/PLLlZWVdZHdQ0dYIA8AEElMB5b58+fr5MmTWrNmjTwej6ZNm6aSkpLARNzjx48rJqb7I012u11vv/22tm7dqrNnz2r06NG6+eab9fDDD3MVpQ+1LZDn8TaFnMdiU+tt3CyQBwCwApthGBG/6EZ9fb2cTqe8Xq/i4+PD3ZyI0XaXkKSg0NI2xZa7hAAAfcnM72+eJTSADbQF8gAAkatPb2uG9Q2UBfIAAJGNwIKoXyAPABD5GBICAACWR2ABAACWR2ABAACWR2ABAACWR2ABAACWR2ABAACWR2ABAACWR2ABAACWR2ABAACWFxUr3bY9v7G+vj7MLQEAAN3V9nu7O89hjorA0tDQIElKTk4Oc0sAAIBZDQ0NcjqdndbYjO7EGovz+/367LPPNHToUNls0f3Qvvr6eiUnJ+vEiRNdPoo7GtH/gdv/gdx3if7T/+jsv2EYamho0OjRoxUT0/kslai4whITE6PLLrss3M3oV/Hx8VH1RWsW/R+4/R/IfZfoP/2Pvv53dWWlDZNuAQCA5RFYAACA5RFYIkxcXJwKCgoUFxcX7qaEBf0fuP0fyH2X6D/9H9j9l6Jk0i0AAIhuXGEBAACWR2ABAACWR2ABAACWR2ABAACWR2AJo7/85S/KycnR6NGjZbPZ9Lvf/a7T+rKyMtlstnYvj8cTVLd+/XqNGzdODodD6enpqqio6MNe9JzZ/t95550h+3/VVVcFan70ox+1e//KK6/s456YV1RUpG984xsaOnSoEhMTNW/ePB05cqTL/X7zm9/oyiuvlMPh0OTJk7Vnz56g9w3D0Jo1a+R2u3XJJZcoMzNTH3zwQV91o8d60v+NGzfquuuu0/DhwzV8+HBlZma2+9oO9TWSnZ3dl13pkZ70f8uWLe365nA4gmqi+fzfcMMNIb//v/WtbwVqIuH8b9iwQVOmTAksAJeRkaE//OEPne4TLd/3F4vAEkaNjY2aOnWq1q9fb2q/I0eOqKamJvBKTEwMvLdjxw7l5eWpoKBAlZWVmjp1qrKyslRXV9fbzb9oZvv/1FNPBfX7xIkTGjFihG699daguquuuiqo7o033uiL5l+UvXv3auXKldq/f7/+9Kc/6YsvvtDNN9+sxsbGDvfZt2+fFi5cqGXLlumtt97SvHnzNG/ePFVVVQVqnnjiCT399NMqLi7Wm2++qUsvvVRZWVlqamrqj251W0/6X1ZWpoULF+r1119XeXm5kpOTdfPNN+vTTz8NqsvOzg46/88//3xfd8e0nvRfal3l9Kt9+/jjj4Pej+bzv3PnzqC+V1VVyW63t/v+t/r5v+yyy/TYY4/pwIED+p//+R/deOONmjt3rt55552Q9dH0fX/RDFiCJOOll17qtOb11183JBl/+9vfOqxJS0szVq5cGfjY5/MZo0ePNoqKinqppX2jO/2/0EsvvWTYbDbjo48+CmwrKCgwpk6d2ruN6wd1dXWGJGPv3r0d1vzLv/yL8a1vfStoW3p6uvFv//ZvhmEYht/vN1wul/Hkk08G3j979qwRFxdnPP/8833T8F7Snf5f6MsvvzSGDh1qbN26NbBtyZIlxty5c/ughX2rO/1/9tlnDafT2eH7A+38/+xnPzOGDh1qnDt3LrAtUs//8OHDjf/4j/8I+V40f9+bxRWWCDRt2jS53W7Nnj1bf/3rXwPbW1padODAAWVmZga2xcTEKDMzU+Xl5eFoap/atGmTMjMzNXbs2KDtH3zwgUaPHq3x48fr9ttv1/Hjx8PUwu7zer2SpBEjRnRYU15eHnRuJSkrKytwbqurq+XxeIJqnE6n0tPTLX/+u9P/C33++ef64osv2u1TVlamxMRETZw4Ubm5uTp9+nSvtrUvdLf/586d09ixY5WcnNzur/KBdv43bdqkBQsW6NJLLw3aHknn3+fzafv27WpsbFRGRkbImmj+vjeLwBJB3G63iouL9dvf/la//e1vlZycrBtuuEGVlZWSpFOnTsnn8ykpKSlov6SkpHbzXCLdZ599pj/84Q9avnx50Pb09HRt2bJFJSUl2rBhg6qrq3XdddepoaEhTC3tmt/v17333qtrr71WkyZN6rDO4/F0em7b/htp57+7/b/QAw88oNGjRwf9oM7OztZzzz2n0tJSPf7449q7d6/mzJkjn8/XF03vFd3t/8SJE7V582a9/PLL2rZtm/x+v2bOnKlPPvlE0sA6/xUVFaqqqmr3/R8p5//QoUP6u7/7O8XFxenb3/62XnrpJaWmpoasjdbv+56Iiqc1DxQTJ07UxIkTAx/PnDlTx44d089+9jP953/+Zxhb1v+2bt2qYcOGad68eUHb58yZE/j3lClTlJ6errFjx+qFF17QsmXL+rmV3bNy5UpVVVVZcq5Nf+hJ/x977DFt375dZWVlQRNPFyxYEPj35MmTNWXKFE2YMEFlZWW66aaberXdvaW7/c/IyAj6K3zmzJn6+te/rl/+8pd6+OGH+7qZfaYn53/Tpk2aPHmy0tLSgrZHyvmfOHGiDh48KK/XqxdffFFLlizR3r17OwwtaMUVlgiXlpamo0ePSpISEhJkt9tVW1sbVFNbWyuXyxWO5vUJwzC0efNm3XHHHYqNje20dtiwYbriiisCnyOrufvuu/X73/9er7/+ui677LJOa10uV6fntu2/kXT+zfS/zU9+8hM99thjevXVVzVlypROa8ePH6+EhISoOP8XGjx4sK6++upA3wbK+W9sbNT27du79QeIVc9/bGysLr/8cs2YMUNFRUWaOnWqnnrqqZC10fh931MElgh38OBBud1uSa3fBDNmzFBpaWngfb/fr9LS0g7HRyPR3r17dfTo0W79wDp37pyOHTsW+BxZhWEYuvvuu/XSSy/ptddeU0pKSpf7ZGRkBJ1bSfrTn/4UOLcpKSlyuVxBNfX19XrzzTctd/570n+p9W6Ihx9+WCUlJbrmmmu6rP/kk090+vTpqDj/F/L5fDp06FCgbwPh/Eutt/g2Nzdr0aJFXdZa9fxfyO/3q7m5OeR70fR9f9HCOuV3gGtoaDDeeust46233jIkGWvXrjXeeust4+OPPzYMwzBWrVpl3HHHHYH6n/3sZ8bvfvc744MPPjAOHTpk3HPPPUZMTIzx5z//OVCzfft2Iy4uztiyZYtx+PBh46677jKGDRtmeDyefu9fV8z2v82iRYuM9PT0kMf87ne/a5SVlRnV1dXGX//6VyMzM9NISEgw6urq+rQvZuXm5hpOp9MoKyszampqAq/PP/88UHPHHXcYq1atCnz817/+1Rg0aJDxk5/8xHj33XeNgoICY/DgwcahQ4cCNY899pgxbNgw4+WXXzbefvttY+7cuUZKSorxv//7v/3av670pP+PPfaYERsba7z44otB+zQ0NBiG0fr1dP/99xvl5eVGdXW18ec//9mYPn268bWvfc1oamrq9z52pif9LywsNP74xz8ax44dMw4cOGAsWLDAcDgcxjvvvBOoiebz32bWrFnG/Pnz222PlPO/atUqY+/evUZ1dbXx9ttvG6tWrTJsNpvx6quvGoYR3d/3F4vAEkZttylf+FqyZIlhGK236F1//fWB+scff9yYMGGC4XA4jBEjRhg33HCD8dprr7U77s9//nPj//2//2fExsYaaWlpxv79+/upR+aY7b9htN6ud8kllxi/+tWvQh5z/vz5htvtNmJjY40xY8YY8+fPN44ePdrHPTEvVL8lGc8++2yg5vrrrw98Ltq88MILxhVXXGHExsYaV111lbF79+6g9/1+v7F69WojKSnJiIuLM2666SbjyJEj/dAjc3rS/7Fjx4bcp6CgwDAMw/j888+Nm2++2Rg1apQxePBgY+zYscaKFSssGdZ70v9777038H2dlJRk3HLLLUZlZWXQcaP5/BuGYbz33nuGpMAv96+KlPP/r//6r8bYsWON2NhYY9SoUcZNN90U1J9o/r6/WDbDMIy+unoDAADQG5jDAgAALI/AAgAALI/AAgAALI/AAgAALI/AAgAALI/AAgAALI/AAgAALI/AAgAALI/AAgAALI/AAgAALI/AAgAALI/AAgAALO//A5mfbL0H8kBdAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot individual parameter estimates against the individual noise level\n", "r_indiv = Mflex.get_correlation(theta[0])\n", "fSNR_indiv = Mflex.get_fSNR(theta[0], n_part, separate=False)\n", "plt.scatter(fSNR_indiv, r_indiv)\n", "\n", "# Add the group estimate as a horizontal line\n", "theta_g,_= pcm.group_to_individ_param(theta_gr[0],Mflex,n_subj)\n", "r_group = Mflex.get_correlation(theta_g)\n", "plt.axhline(r_group[0], color='red', linestyle='--')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inference via bootstrap\n", "Again, we can use the bootstrap to get confidence intervals for the group-level correlation parameter\n", "In this case we need to specify the fixed (condition) effects. For this example we draw 200 bootstrap samples-in practice we recommend 5000 or more. " ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bootstrapping 200 samples: \n", "Bootstrap sample 0/200\n", "Bootstrap sample 50/200\n", "Bootstrap sample 100/200\n", "Bootstrap sample 150/200\n" ] } ], "source": [ "\n", "r_boot,fSNR_boot,_ = pcm.bootstrap_group_corr(D,Mflex,fixed_effect=X,n_bootstr=200,verbose=True)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "95% confidence interval for the group-level correlation parameter: [0.609, 0.685]\n" ] } ], "source": [ "# From the bootstrap samples, we can get an approximate confidence interval using the bootstrap method.\n", "CI_low = np.percentile(r_boot, 2.5)\n", "CI_high = np.percentile(r_boot, 97.5)\n", "print(f\"95% confidence interval for the group-level correlation parameter: [{CI_low:.3f}, {CI_high:.3f}]\")" ] } ], "metadata": { "kernelspec": { "display_name": "base", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.5" } }, "nbformat": 4, "nbformat_minor": 4 }