{
"metadata": {
"name": "",
"signature": "sha256:06aaaf7085ddd47dd0b96d93bf3c322d5569eeaac8e2cb804e1a926e01aa02e8"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The Generalised Linear Model and ANOVA\n",
"\n",
"[Tom Wallis](http://www.tomwallis.info) and [Philipp Berens](http://philippberens.wordpress.com/)\n",
"\n",
"The University of T\u00fcbingen"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy as np\n",
"import pandas as pd\n",
"import patsy\n",
"import seaborn as sns\n",
"import statsmodels.formula.api as smf\n",
"import statsmodels.stats as sms\n",
"import matplotlib.pyplot as plt\n",
"\n",
"%matplotlib inline\n",
"\n",
"sns.set_style(\"white\")"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 21
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introduction\n",
"\n",
"The analysis of variance (ANOVA) is typically used to \"*learn the relative importance of different sources of variation in a dataset*\" (see Gelman & Hill, Chapter 22). Recall that when we think about regression models, we ask questions like: \n",
"\n",
"**What is our prediction for distance flown given paper and design?**\n",
"\n",
"The ANOVA for the same model fit focuses on a slightly different question:\n",
"\n",
"**How much of the variance in distance flown is associated with changes in paper and design?**\n",
"\n",
"ANOVA can be thought of as *breaking up* the total variance in the response into (a) how much variance is associated with each independent variable in our design and (b) the remaining *error* or *residual* variance. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# for this demonstration we will use the *planes* example because it's balanced.\n",
"\n",
"# read in the data, rename some values:\n",
"planes = pd.read_csv('../Datasets/planes.txt', sep='\\t')\n",
"\n",
"# column names to lowercase:\n",
"planes.columns = [c.lower() for c in planes.columns]\n",
"\n",
"# turn \"paper\" variable into an object, not int64:\n",
"planes['paper'] = planes['paper'].astype(object)\n",
"\n",
"planes.replace(to_replace={'paper': {1: '80',\n",
" 2: '50'},\n",
" 'angle': {1: 'horz',\n",
" 2: '45deg '},\n",
" 'plane': {1: 'sleek',\n",
" 2: 'simple'}}, \n",
" inplace=True)\n",
"\n",
"# convert distance into meters for ease of SSR:\n",
"planes['distance'] = planes['distance'] / 1000\n",
"\n",
"# rename \"plane\" to \"design\":\n",
"planes = planes.rename(columns={'plane': 'design'})\n",
"\n",
"# planes is a balanced dataset (same number of cases in each cell):\n",
"planes.groupby(['paper', 'design']).size()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 23,
"text": [
"paper design\n",
"50 simple 4\n",
" sleek 4\n",
"80 simple 4\n",
" sleek 4\n",
"dtype: int64"
]
}
],
"prompt_number": 23
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ANOVA as model comparison between nested models\n",
"\n",
"Recall for the planes data we are interested in the effects of `design` and `paper` on distance. This is a $2 \\times 2$ factorial design. We can think about all the possible models here as a series of nested models. That is, each simpler model is a subset of the most complex model (the model with all main effects and interactions). We can visualise this in a diagram (adapted from [1]):\n",
"\n",
"
\n",
"\n",
"In the diagram, the arrows indicate a differencing of the sum of squared residuals (SSR). The SSR for the upper model (with more parameters, and thus a lower SSR) is subtracted from the SSR of the lower model. This difference score gives the change in SS associated with the new predictor.\n",
"\n",
"Notice that there are actually two comparisons here that speak to the \"main effects\" of `paper` (in red) and `design` (in blue). These correspond to different hypotheses. Let's consider the test for `paper`. The first comparison (red text, right side of figure) tests whether paper has an effect *ignoring* design. This is usually called the **main effect** of paper. The second comparison (red text, left side of figure) tests whether adding paper changes the residuals *after discounting the variance due to design*. This corresponds to the SS for the **simple effects** of paper (in an interaction model, we also add the SS from the interaction term -- see [2], p.433).\n",
"\n",
"The `planes` dataset we use here is an example of a *balanced design*. Balanced designs have an equal number of observations in each cell (four planes were thrown in each combination of `paper` and `design`). It turns out that in balanced designs, the main effects and the sum of simple effects give the same answer. That's because they're just different ways to divide up the total sum of squares. When we have unbalanced data however, it matters which comparison you do (see Types of Sums of Squares, below). "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Orthogonal contrasts\n",
"\n",
"When we do ANOVA, we are breaking up the total variance from the grand mean into subgroups of variance. To do this in a way that makes sense, we need to code categorical variables in a way that is independent of the other columns. These are called *orthogonal contrasts*. \n",
"\n",
"What is an orthogonal contrast? Intuitively, in a design matrix of orthogonal contrasts, the overall effect of one column is independent of the other columns. Mathematically, two columns of a design matrix are orthogonal if their product sums to zero.\n",
"\n",
"Throughout this course so far, recall that we have coded categorical variables using *treatment* or *dummy* coding. In this case, the Intercept term is the mean of one of the subgroups in the design (the *reference* level) and all other coefficients give deflections from this reference level. Are treatment contrasts orthogonal?\n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"X = patsy.dmatrix('~ design * paper', planes)\n",
"X"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 24,
"text": [
"DesignMatrix with shape (16, 4)\n",
" Intercept design[T.sleek] paper[T.80] design[T.sleek]:paper[T.80]\n",
" 1 1 1 1\n",
" 1 1 1 1\n",
" 1 0 1 0\n",
" 1 0 1 0\n",
" 1 1 1 1\n",
" 1 1 1 1\n",
" 1 0 1 0\n",
" 1 0 1 0\n",
" 1 1 0 0\n",
" 1 1 0 0\n",
" 1 0 0 0\n",
" 1 0 0 0\n",
" 1 1 0 0\n",
" 1 1 0 0\n",
" 1 0 0 0\n",
" 1 0 0 0\n",
" Terms:\n",
" 'Intercept' (column 0)\n",
" 'design' (column 1)\n",
" 'paper' (column 2)\n",
" 'design:paper' (column 3)"
]
}
],
"prompt_number": 24
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# example for coefficients of design and paper:\n",
"np.sum(X[:, 1] * X[:, 2])"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 25,
"text": [
"array(4.0)"
]
}
],
"prompt_number": 25
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So in the *treatment* contrasts regime, we're not using orthogonal contrasts (the sum of the product of the columns is not zero). Therefore ANOVAs done with these design matrices become hard or impossible to interpret.\n",
"\n",
"The fix for this is easy: we can tell Patsy to code our categorical variables using orthogonal contrasts. The most common of these is Sum coding, in which the values for each category sum to one. This means that the coefficients in the regression represent *deviations* from the overall mean across the category levels. There are a number of other contrasts available in Patsy, which can test different hypotheses about the variables. See [here](http://statsmodels.sourceforge.net/devel/contrasts.html), [here](http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm) and [here](http://talklab.psy.gla.ac.uk/tvw/catpred/) for further reading.\n",
"\n",
"What does the design matrix look like for sum coding?"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"X = patsy.dmatrix('~ C(design, Sum) * C(paper, Sum)', planes)\n",
"print(X)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[[ 1. -1. -1. 1.]\n",
" [ 1. -1. -1. 1.]\n",
" [ 1. 1. -1. -1.]\n",
" [ 1. 1. -1. -1.]\n",
" [ 1. -1. -1. 1.]\n",
" [ 1. -1. -1. 1.]\n",
" [ 1. 1. -1. -1.]\n",
" [ 1. 1. -1. -1.]\n",
" [ 1. -1. 1. -1.]\n",
" [ 1. -1. 1. -1.]\n",
" [ 1. 1. 1. 1.]\n",
" [ 1. 1. 1. 1.]\n",
" [ 1. -1. 1. -1.]\n",
" [ 1. -1. 1. -1.]\n",
" [ 1. 1. 1. 1.]\n",
" [ 1. 1. 1. 1.]]\n"
]
}
],
"prompt_number": 26
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Instead of our levels being either 0 or 1, they're now either -1 or 1. Now, the sum over our contrast columns is zero (note the intercept is not a contrast, so does not sum to zero), and their cross products are all zero:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print('sum of columns = ' + str(X.sum(axis=0)))\n",
"print('sum of products of design and paper = ' + str((X[:, 1] * X[:, 2]).sum()))\n",
"print('sum of products of design and interaction = ' + str((X[:, 1] * X[:, 3]).sum()))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"sum of columns = [ 16. 0. 0. 0.]\n",
"sum of products of design and paper = 0.0\n",
"sum of products of design and interaction = 0.0\n"
]
}
],
"prompt_number": 27
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using these orthogonal contrasts, let's fit the model and look at the ANOVA table:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fit_full = smf.ols('distance ~ C(design, Sum) * C(paper, Sum)', planes).fit()\n",
"print(fit_full.summary())"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: distance R-squared: 0.727\n",
"Model: OLS Adj. R-squared: 0.659\n",
"Method: Least Squares F-statistic: 10.66\n",
"Date: Tue, 24 Feb 2015 Prob (F-statistic): 0.00106\n",
"Time: 20:51:51 Log-Likelihood: -18.584\n",
"No. Observations: 16 AIC: 45.17\n",
"Df Residuals: 12 BIC: 48.26\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"================================================================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"----------------------------------------------------------------------------------------------------------------\n",
"Intercept 3.6853 0.223 16.514 0.000 3.199 4.171\n",
"C(design, Sum)[S.simple] -0.1552 0.223 -0.696 0.500 -0.641 0.331\n",
"C(paper, Sum)[S.50] 0.3277 0.223 1.469 0.168 -0.158 0.814\n",
"C(design, Sum)[S.simple]:C(paper, Sum)[S.50] -1.2090 0.223 -5.418 0.000 -1.695 -0.723\n",
"==============================================================================\n",
"Omnibus: 0.751 Durbin-Watson: 2.909\n",
"Prob(Omnibus): 0.687 Jarque-Bera (JB): 0.745\n",
"Skew: 0.359 Prob(JB): 0.689\n",
"Kurtosis: 2.224 Cond. No. 1.00\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"prompt_number": 28
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What do the coefficients mean? The intercept coefficient is now the grand mean of distance (recall that in the Treatment contrast regime it was the mean distance when the design was \"simple\" and the paper weight was \"50\"):"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"np.round(planes['distance'].mean(), decimals=4)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 29,
"text": [
"3.6852"
]
}
],
"prompt_number": 29
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The coefficient values for the contrasts tell us the average deviations from this grand mean for each category:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# example for design:\n",
"print(fit_full.params[0] + fit_full.params[1]) # model prediction for simple\n",
"print(fit_full.params[0] + -1*fit_full.params[1]) # model prediction for sleek\n",
"print('\\n\\n Actual Means:')\n",
"print(planes.groupby('design').distance.mean()) # actual means"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"3.53\n",
"3.8405\n",
"\n",
"\n",
" Actual Means:\n",
"design\n",
"simple 3.5300\n",
"sleek 3.8405\n",
"Name: distance, dtype: float64\n"
]
}
],
"prompt_number": 30
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ANOVA table for this model"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"anova_table = sms.anova.anova_lm(fit_full) # type I SS by default; same answer in this case\n",
"print(anova_table) "
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
" df sum_sq mean_sq F PR(>F)\n",
"C(design, Sum) 1 0.385641 0.385641 0.484016 0.499861\n",
"C(paper, Sum) 1 1.718721 1.718721 2.157158 0.167628\n",
"C(design, Sum):C(paper, Sum) 1 23.386896 23.386896 29.352777 0.000156\n",
"Residual 12 9.561029 0.796752 NaN NaN\n"
]
}
],
"prompt_number": 31
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Relating the ANOVA table to nested model comparisons\n",
"\n",
"Recall that above we talked about ANOVA as a comparison of nested models. I'll reproduce the diagram here:\n",
"\n",
"\n",
"\n",
"We can produce all the sums of squares from the ANOVA table by computing the relevant change in residual (error) sum of squares for each model.\n",
"\n",
"### Fit all models in the diagram"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# fit all the models in the tree\n",
"null = smf.ols('distance ~ 1', planes).fit()\n",
"design = smf.ols('distance ~ C(design, Sum)', planes).fit()\n",
"paper = smf.ols('distance ~ C(paper, Sum)', planes).fit()\n",
"additive = smf.ols('distance ~ C(design, Sum) + C(paper, Sum)', planes).fit()\n",
"full = smf.ols('distance ~ C(design, Sum) * C(paper, Sum)', planes).fit()"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 32
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The null model represents the total sum of squared residuals"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print(null.summary())"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: distance R-squared: 0.000\n",
"Model: OLS Adj. R-squared: 0.000\n",
"Method: Least Squares F-statistic: nan\n",
"Date: Tue, 24 Feb 2015 Prob (F-statistic): nan\n",
"Time: 20:51:51 Log-Likelihood: -28.977\n",
"No. Observations: 16 AIC: 59.95\n",
"Df Residuals: 15 BIC: 60.73\n",
"Df Model: 0 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"------------------------------------------------------------------------------\n",
"Intercept 3.6853 0.382 9.643 0.000 2.871 4.500\n",
"==============================================================================\n",
"Omnibus: 0.701 Durbin-Watson: 1.723\n",
"Prob(Omnibus): 0.704 Jarque-Bera (JB): 0.714\n",
"Skew: 0.351 Prob(JB): 0.700\n",
"Kurtosis: 2.239 Cond. No. 1.00\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"prompt_number": 33
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This model explains 0% of the variance in the data. Because it's predicting every data point to be the mean, the residual variance is equal to the total variance:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"null.mse_total # mean squared error total"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 34,
"text": [
"2.336819133333333"
]
}
],
"prompt_number": 34
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"null.mse_resid"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 35,
"text": [
"2.336819133333333"
]
}
],
"prompt_number": 35
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# this is the same as the variance in distance:\n",
"planes['distance'].var()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 36,
"text": [
"2.3368191333333317"
]
}
],
"prompt_number": 36
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# let's calculate the sum of squared residuals manually:\n",
"res = planes['distance'] - null.predict() # absolute residuals (difference between data and prediction)\n",
"res = res ** 2 # squared residuals\n",
"res.sum()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 37,
"text": [
"35.052286999999993"
]
}
],
"prompt_number": 37
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# the ssr is also contained in the model object:\n",
"null.ssr"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 38,
"text": [
"35.052286999999993"
]
}
],
"prompt_number": 38
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Main effect of design\n",
"\n",
"Recall that the SS for each model term can be calculated by comparing the change in the residual sum of squares for nested models. We do that here:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"## Main effect of design, ignoring paper:\n",
"print(null.ssr - design.ssr)\n",
"\n",
"## sum of simple effects of design, eliminating paper:\n",
"print(paper.ssr - additive.ssr)\n",
"\n",
"## main effect from ANOVA table:\n",
"print(anova_table['sum_sq'][0])"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.385641\n",
"0.385641\n",
"0.385641\n"
]
}
],
"prompt_number": 39
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Main effect of paper"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"## Main effect of paper, ignoring design:\n",
"print(null.ssr - paper.ssr)\n",
"\n",
"## sum of simple effects of paper, eliminating design:\n",
"print(design.ssr - additive.ssr)\n",
"\n",
"## main effect from ANOVA table:\n",
"print(anova_table['sum_sq'][1])"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"1.718721\n",
"1.718721\n",
"1.718721\n"
]
}
],
"prompt_number": 40
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Interaction effect"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# interaction term:\n",
"print(additive.ssr - full.ssr)\n",
"\n",
"## interaction effect from ANOVA table:\n",
"print(anova_table['sum_sq'][2])"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"23.386896\n",
"23.386896\n"
]
}
],
"prompt_number": 41
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Partitioning of residuals, variance explained\n",
"\n",
"As we add terms to the model, the residual sum of squares gets lower. It turns out that for any model, its $R^2$ (the coefficient of determination) is given by the proportion of the variance explained, relative to the null model:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# e.g. for full model:\n",
"print(full.summary())"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: distance R-squared: 0.727\n",
"Model: OLS Adj. R-squared: 0.659\n",
"Method: Least Squares F-statistic: 10.66\n",
"Date: Tue, 24 Feb 2015 Prob (F-statistic): 0.00106\n",
"Time: 20:51:51 Log-Likelihood: -18.584\n",
"No. Observations: 16 AIC: 45.17\n",
"Df Residuals: 12 BIC: 48.26\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"================================================================================================================\n",
" coef std err t P>|t| [95.0% Conf. Int.]\n",
"----------------------------------------------------------------------------------------------------------------\n",
"Intercept 3.6853 0.223 16.514 0.000 3.199 4.171\n",
"C(design, Sum)[S.simple] -0.1552 0.223 -0.696 0.500 -0.641 0.331\n",
"C(paper, Sum)[S.50] 0.3277 0.223 1.469 0.168 -0.158 0.814\n",
"C(design, Sum)[S.simple]:C(paper, Sum)[S.50] -1.2090 0.223 -5.418 0.000 -1.695 -0.723\n",
"==============================================================================\n",
"Omnibus: 0.751 Durbin-Watson: 2.909\n",
"Prob(Omnibus): 0.687 Jarque-Bera (JB): 0.745\n",
"Skew: 0.359 Prob(JB): 0.689\n",
"Kurtosis: 2.224 Cond. No. 1.00\n",
"==============================================================================\n",
"\n",
"Warnings:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"prompt_number": 42
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"diff = null.ssr - full.ssr\n",
"# express as a proportion of the total (null model) residual:\n",
"prop = diff / null.ssr\n",
"print(prop)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.727235230044\n"
]
}
],
"prompt_number": 43
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# same as R-squared:\n",
"print(full.rsquared)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.727235230044\n"
]
}
],
"prompt_number": 44
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# F-tests and t-tests\n",
"\n",
"As you may know from your basic statistics courses, the sum of squares for each term is divided by how many parameters were estimated for that term (the degrees of freedom added) to give the mean squared column, and then the ratio between the change and the amount of remaining error gives the F-statistic. This can then be tested against the F-distribution (with `model, error` degrees of freedom) for significance."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print(anova_table)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
" df sum_sq mean_sq F PR(>F)\n",
"C(design, Sum) 1 0.385641 0.385641 0.484016 0.499861\n",
"C(paper, Sum) 1 1.718721 1.718721 2.157158 0.167628\n",
"C(design, Sum):C(paper, Sum) 1 23.386896 23.386896 29.352777 0.000156\n",
"Residual 12 9.561029 0.796752 NaN NaN\n"
]
}
],
"prompt_number": 45
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the regression table for the full model above. How are the t-statistics related to the F-values for the model terms? \n",
"\n",
"In the $2 \\times 2$ design here, the t-tests are equivalent to the F-tests for the factors... their values are the square roots of the F-values (i.e. $F = t^2$)"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"full.tvalues[1:] **2"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 46,
"text": [
"C(design, Sum)[S.simple] 0.484016\n",
"C(paper, Sum)[S.50] 2.157158\n",
"C(design, Sum)[S.simple]:C(paper, Sum)[S.50] 29.352777\n",
"dtype: float64"
]
}
],
"prompt_number": 46
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print(anova_table['F'])"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"C(design, Sum) 0.484016\n",
"C(paper, Sum) 2.157158\n",
"C(design, Sum):C(paper, Sum) 29.352777\n",
"Residual NaN\n",
"Name: F, dtype: float64\n"
]
}
],
"prompt_number": 47
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# One way to visualise variance partitioning\n",
"\n",
"Recall that when we do ANOVA, we're partitioning the total variance into the variance associated with each factor, and that this is equivalent to nested model comparison. We can visualise this as a pie chart:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# plot the explained variances as a pie chart:\n",
"labels = 'design', 'paper', 'interaction', 'residual'\n",
"sizes = anova_table['sum_sq']\n",
"colors = ['yellowgreen', 'gold', 'lightskyblue', 'lightcoral']\n",
"plt.pie(sizes, labels=labels,\n",
" autopct='%1.0f%%',\n",
" colors=colors,\n",
" startangle=90)\n",
"# Set aspect ratio to be equal so that pie is drawn as a circle.\n",
"plt.axis('equal')\n",
"plt.show() "
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {
"png": {
"height": 253,
"width": 349
}
},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAroAAAH6CAYAAADoTe/oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAWJQAAFiUBSVIk8AAAIABJREFUeJzs3Xd4VFX+BvD3TJ+U6ZNM2qRzE0gghF4CoQVQRKrSm4od\nELuufS1rw7I/+6KuvXfXspZVxN4LXhUBaVKSkADpyf39cScYIKEmuVPez/PkGTJz595vYpy8OXPO\n9whFUUBEREREFG50WhdARERERNQRGHSJiIiIKCwx6BIRERFRWGLQJSIiIqKwxKBLRERERGGJQZeI\niIiIwhKDLhERERGFJQZdIiIiIgpLDLpEREREFJYYdImIiIgoLDHoEhEREVFYYtAlIiIiorDEoEtE\nREREYYlBl4iCgiRJWZIkNUmS9GAHXmNe4BqLOuoaREQUPBh0iSjYKB147q8BXAngkw68BhERBQmD\n1gUQEXUWWZa/BfCt1nUQEVHn4IguEREREYUljugSUaeTJCkPwDUAhkD9g/sVAHe3cpwJwLkAZgNI\nB1AJ4G0Al8myvGafY6cDOAtAV6ivbasALJdl+Z4Wx8wDsBzAObIs397i/uEALgfQE0AtgOcB/BPA\ndwCukmX5qsBx7wNIBTAYwE0ARgOwAPgCwOWyLP/vyL8rRETU3jiiS0SdSpKkAgAfATgGwBsAHgVQ\nDOC5fY4zAvgPgGsBVAC4M3D8ZACfS5LUrcWx0wA8BsANNcjeA8AJ4C5Jkv7WShlKi+dOAvAWgHwA\nTwN4CsAUAC/te2zg3zEAPgwc/yCAFwEMAvCmJEldD+d7QUREHYsjukTU2W4HEAVgjCzL7wCAJElX\nAHgPgK/FcUsADAPwD1mWL26+U5KkOwCshBpo+wXuPg/ALgC9ZFneHTjuagA/Qx3l/XtrhUiSFA3g\nLgA7APSTZXl14P4bAXzVylME1DD9IYCpsiw3Bo7/AWognw3g4laeR0REGuCILhF1GkmSkgAUAXiz\nOeQCgCzLZQCu2OfwkwCUA7i05Z2yLH8J4BkAfSRJyg3crYManvNbHLcTQB+oUx7aMhpAHIB/Nofc\nwHPXA7j1AM+7pTnkBvwncJt6gOcQEVEn44guEXWm7oHbz1p5bGXzPwIjrV0A/AngckmS9j22eeS3\nJ9S5uPcEPlZKkvQ9gNcDHytkWT5Qu7I+h1LPPhQAv+xzX0Xg1nyAaxERUSdj0CWizuQM3O5s5bGy\nFv+2B259UBeJtUZpPp8sy/dJkrQVwCKoI8b5AC4EsFGSpKWyLD/Txjk8gds/W3lsUxvPAdQFa/vW\nAqhTG4iIKEgw6BKFIUmSHgIwB0A81LfgjwPQAHVu6aWyLP/U4thUABcBKAGQGDhOBnC/LMv3tjhu\nHtR5sWMA9AdwKgAbgO8BXCvL8mut1DEVwDlQR3IbAfweeMje4phiAO8CWAo1KM4BMDbw8AeyLBcf\nytcsy/KLAF6UJMkOdW7veAAzADwhSdJPsiz/2MrTKgO3tlYea+0+IiIKIQy6ROHtDahzUB8AkAJg\nEoBiSZKGyrL8nSRJaQA+B2CF2lJrPYBkqF0H7pYkySDL8v/tc87roLbwehRqeJ0K4GVJkk6SZfmh\n5oMCi8H+BmAN1IAMACdCDbNToO5Q1tJlgdufAKyDGo7zJEmyyLJc0/LAQCuxbAAPA9gGdTHaTlmW\nl8myXAG1E8KLkiT9DuBqAAMAtBZ0vwjc9gPw/j6P9QMREYU0LkYjCm8OAD1kWT5PluUToQZNO9TO\nB4A6kusCMF6W5TmyLF8qy/JcqO2+AGB6K+fsDqBYluWFsiyfDqAv1DmqywKjqZAkqS/UkPsegG6y\nLC+SZXkRgFwAuwF0C4wQ71urAuALWZaPg9q6ywXgBkmS9kwJCLTwugfqCHCpLMtVABYAuEqSpH0X\nnjV/vq6N789LUKdMLAqE/uZrJAO4oI3nEBFRiGDQJQpv1wQ6GgAAZFl+HmoP26GSJCUCeATAAlmW\n3235JFmWPwdQA3U0eF+Py7L8WYtjf4e6uYIdwLGBuxcEbs9vORobqOWawKfLJUl6Dmr7LxG4Xks3\nAPgU6rzbLyRJulWSpOWB+6IBnCrL8q7AsRdB7W/7lSRJ90mS9A9Jkt4BMB/Ae7Isv93aNycQks+E\nOhf4S0mSHpAk6V4AX0Pt4gCoo9YtcR4uEVGI4NQFovD2Xiv3fQZ1g4Pusiy/AeAjSZJcAAoAZAGQ\noM7BNQPQH+I5Pw/cdgfwOIBegc+nSJI0fp9jk6GO3P4GdeFYdOD+F9FiBFmW5RpJkoYBOB/ANACn\nQ+13+yGA62VZ/rDFsU9KklQJdQrDeKihew3UlmU3tri2gr03gIAsy09JkrQbahuz6VBHnJ8IXOcp\nAFUHej4REQUvBl2i8LaxlfuaOwzYJUlyAlgGddGWAWqIWwN1cVhPtD56ecBzBm4dgduL2qhLAfCt\nLMtTWyxG2yLL8l7vMgVGg6/BX6PAbZJlubml2IGOeRjqvF4AgCRJsQBssiy/CuDVlsdKkjQ/8M/1\nLZ4/rI3zrgXfISMiCjoMukThzYr9W3k1h9BSqAvKxgK4G+o0hu9b7Cw2+wDn3FfzObcHbndB7d5g\n3WdjhWAjAfhMkqSHZFlesOdOSbJCndJQD2CFVsUREdHRYdAlCm99Abyzz30DoAa4L6GG3M9lWT6z\n5QGBhVlmtD6i2xfAK62cE1DnzwLAtwB6ACjEX9Mams/dB8BEAK/Lsqx1iPwS6lSOeS06UEQBGAfA\nD7UVW2s9domIKATwrTai8Ha1JEkxzZ9IkjQFakeFl6Au/lIAOCVJMrY4xgp1cRkAGLG/0yVJym5x\nvAR1wdgGAG8F7n4ocLssMD2g+dgYAPdCndKg+etPYNe0EqgtyOKhjuLOALAWwFRZlm/QrjoiIjpa\nQlG4roIo3LTYMKIC6hSF16AuApsAdY7tQFmWN0iS9DTUnrbfAngbaueC46CG4EaovXdjZFlWWmwY\n0bzd7bNQw+pkqKO/E2RZfrNFDbfhrwD8OtTdxCYG6ri7eRS5xRzd22RZXtr+3w0iIopUmo+oEFGH\nOgHAN1DbfQ2E2pu2nyzLGwKPnwTgNqhzbBcBGAl1sVYh1M0mLPirp26zqwHcD7W7wSQAKwEMbRly\nAUCW5SUAZkNdzDULavDeBLXl11nt+DUSERG1iiO6RGGoxYhuXsvtfo/ynPOgjuieJcvyXe1xTiIi\noo7EEV0iIiIiCksMukREREQUlhh0icJTR+zgxV3BiIgopHCOLhERERGFJY7oEhEREVFYYtAlIiIi\norDEoEtEREREYYlBl4iIiIjCkkHrAoiIwpEQQg91Zzkz1NdaAXVwwQZ1G+WmwEdji39XKYrSqEnB\nRERhiF0XiIhaIYQwAHAD8ALw2mKQ6HEizWyCv7JG10NEm3cZjMJsMAqT3iCMOr0w6vXCoNPDIIQw\n6I1CZzLrhNEsdHqDEDq9gBCA/Hate0y6VC6EaGpSFKWhqQlNTU1KQ1MTdtXVNdU3NtY3NDXV1Tc2\n1tc33zY21tc1NtY1KcpuAOVNilJe19i49c+dO+WahoY/AGwEsFVRlCYtv2dERMGGQZeIIpIQQgDw\nAOiS4EUPpw19TUak2WJgi45CVJQFlgQvjEnxMCfFIdrngcXrArwuYNI5ourkB3pGHcl1b5/6R/X7\nM0+3GvX6w35ubUMDdtbWYmdtLSpqavDnzp1N63bs2Pl7WVnVhoqK2qq6ut3V9fVVu+vrd9c3Nm5r\naGpavb2q6qcd1dVrAKwDsEHhiz4RRRBOXSCisCaEsAHIdjuQH+dCX4MBXewxcBwzBDZ/IqJ75sLW\nLQsx2alqiBXi4Oc0Gzu87NavazDAbDDAEx3dfJcOgD3wsZedtbXYXFmJTZWV+KOiYvcv27bt/LW0\ndPfA1NQdFTU15fVNTd9vrKhYUVVf/z2ANYqiNHTil0JE1CkYdIkoLAgh3GYT+qcmothkRLfYaLhs\nMbDPnYCYghzE5neBvUsakBQP6I5yGa7JgEOIw9qKNZsR6/Wii9cLANGBDwBAY1MT1pWXj/x527ZF\nX2/aVPbT1q27BqelVVTU1JTXNTbKW3btWlFRU/M5gF85AkxEoYxBl4hCjhDCDKAgMwVjoqwY7HHC\nd9JkOEcOhKcwF+b0ZMDYgaOuRoOCpqYm6I42MWtEr9Mhw+1GhtutPyYnxwt1HjIURcHGysphP2/b\ntvDTP/7Y/vWmTRUFiYnbquvrv1tXXv5ybWPjJ4qilGtcPhHRIWPQJaKgJ4TwuewYkeDFREcs0iaX\nwF3cB87BvWDPywYMnfxK5ogFdpU3wOY2de6FO5gQAsl2O5Ltdt3IrKw4AHGKomSvr6gY+Nn69fPf\nW716W1F6ellZVdXmnbW1766vqHgLwA+c9kBEwYpBl4iCSmCRWFqCF2PcDoxzO5By+jS4jx2KuEGF\nMDhsWlcIxLuhbN9YF3ZBtzVCCPgdDvgdDvOU/PxkAMl1DQ3df9iyZfSKtWsv+viPP8p7Jydv31lb\n+/kv27c/AuALtkgjomDBoEtEmhNCOHweTI5z4cQR/ZFc2A2usUWI698DwmrRurr9+TwQv2+pQ4bW\nhWjEZDCgMCkJhUlJzkWDBjkBZKwpK+v7xi+/zHrnt9+2dU9IWF9eXf3KhoqK5xRF+UPreokocjHo\nEpEmhBD+tCTMddlx7MxxSJg1HonD+sFgDoFB0gQvDF9trde6jKCS7nLh9P797af372+vb2zM+nLj\nxuKXf/rpgkFpadvKqqp+2VBR8ciuurp3FEXZpXWtRBQ5GHSJqFMEpiT0kNKx0GXHgEWzED9rPBJ6\ndTv6LgidLd4NY+VPdY0ADr8ZbgQw6vXo7/eL/n6/D4CvvKoq/73ffx/3yqpVfxYmJW3dUV399pry\n8vsURVmnda1EFN4YdImowwghjBYzhmYkY+HQPsgbNRDx046BK9OvdWVHx2UHqsvr6sGge0icUVGY\nlJdnnpSXl6ooSuoPW7b0efTrrxf0Sk7eXF5V9caa8vJ7GXqJqCMw6BJRuxJCxHqcOD7Bi9ljipAx\naRQSJo5EtMepdWXtx2UHairruN3uERBCIN/nwz/GjvUpiuL7YcuWno9+/fX83snJm8uqq99cU1Z2\nr6Ioa7Wuk4jCA4MuER01IUR0UjzmxLkwc+oYJM84Fgklg2CKsmpdWcdw2oD63WwscLT2Db0/qqF3\nXq/k5M3lDL1E1A4YdInoiAghhMmIAV3ScOG4YuSdPh0powfBqI+AN/MdNqC+pjHEZhYHNyEE8nw+\n3NAi9D72zTfzCpOSNm3bvfuRDRUVD3AhGxEdLgZdIjosQghfph9L+vfAsZNLkDx/IhzuMJqWcCj0\nekCv49a4HaU59F4/ZoyvqanJ9+Hatd3v/fTTRV3j479bXVp6XV1j4+fcmpiIDgWDLhEdlBBCb4mx\nz/e5as6aPxHus2cjqWcuhNZ1ackYASPXwUCn02FoRoZhaEZGellVVfrDX3455NWff16f5nQ+s27H\njv/jlsREdCAMukTUJiGEz+3PujitsKik9/i5/i+fvlG3/LpfgnALh85nMnBEsbO5oqJwTlGRc8ng\nwc7PNmzods8nn5zcLT5+1bry8ht319e/z1FeItoXgy4R7UUIIcxRsSOdiakXFY6b3aVozjnJvqw8\nAMD67z6q+vCLX1DUW+Mig4CJr56aEUKgX0qKvl9KSmplTU3qE99+O+C5H37YkOFyvbimvPwmRVEq\ntK6RiIIDX6qJCAAghDA7EvxnpuT1WdDz2FkphcfNtpmjYvY6ZvCc86POu+71qk8f3hKlUZlBw2RU\nInrqRrCwWSw4tV8/+6n9+tm/3Lgx99YPP5wpeb0rf9m+/W+KoqzRuj4i0haDLlGEE0LEuJIzLkgv\nLDpx6Pzz07oMHG1SNzHbn8efhdKGNF1NzRZYInwCg9UM1FY3wWxl84Vg0SspSffYtGlp63fsSLt1\nxYribvHxP68pL7+8qq7uI61rIyJtMOgSRSghhMvtz7oyu//IY4Yv/FtaWsHAQ1peVTjhVN3ld39d\nf+M5dcaOrjGYeV1QSjfXIjEjTJsFh7AUhwPLxo1L3Flbm/jAZ5/17J6Q8PvmnTuv27579wucx0sU\nWQT/nyeKLEKIBG+adIM3TRoy8rTLUxO6dD+st+Ab6mrx4IKeNauf+y2ix3QvvxPVO6RMa/5g+2E9\n7/apf1S/P/N0a0Q0HA4S9Y2NePLbb3cu/+KLdVt27frnpsrK5Yqi1GtdFxF1PI7oEkUIIURmXHru\nTQXHTO89YuFlKR5/1hGdx2Ayw5XVv+nTb39Dvx7tXGQI8bmhW7O1Tusy6BAY9XrMLiyMnVlQkPef\nX3657Y6PPjonzel8eN2OHbcoisL/iERhjEGXKMzpDcbu3jTppn5TFuYPP+WSBJs38ajPWTT3gqil\nN79Z9dHybRG7KM3ngXHHj3UKENn9hEOJTqfDsTk5lmNzcqR3V6++8ob335+b4nDcsaGi4l5FUbin\nM1EYYtAlClNGs7XA48+6ffCsJTlD550XF2V3tdu549Jz8GdNqq6ubhtMpnY7bUjxOKHbXVpfByBC\nvwOhbXhmpmlYRob06qpVN9+6YsUZCTbb9X/u3Pk45/AShRcuFyYKM0KIxPiMri/1n3rqG6c9+L8h\nYxdf164ht1nBcSeLv99nbGj3E4cIlx2orajjKGAIE0LguK5drf89+eRul48YcU93n+9rb3T0ONFW\n2xEiCjkMukRhQggR40nN/mf+yMkfL7jrtfHHnntTvDk6tsOuVzBmuvmx9/wRu6DHaQdqdzVw9C8M\n6HU6zCgoiHnnlFN6nDd06CN58fGfOqzWIq3rIqKjx6kLRCFOCKF3JPgXp/cqOv3YpTelJ+UWdspy\nfqPFCkdab3yzajUKcjvjisHFZQfqqho48hdGjHo9TuvXzzG3sLDPnStXvtAtPn7V6rKyRTX19V9r\nXRsRHRmO6BKFsBhX3LjEnJ7fjDvvlr+fct9/szor5DYrmnehdfGtnurOvGawsFoANHJANxxZjUZc\nMHSo+40FCwaf0qfPG5LX+4IQwqt1XUR0+Bh0iUKQ3mDM82XnrRg67/xHznh4RV63YcdbtZhW6MvK\nw/rKFNEQoTN1TQYuXApnNosFfx89Ou7pmTMnlGRnf5bucv1NCMF3QolCCIMuUQgRQvjiMnKf6zdl\n4dun/ev9QUWzlzj0Rm03KOt+7Hz8Y7k+IqOuUc+cGwnSnE48PXNm2o1jx/6te0LC157o6OFa10RE\nh4ZBlygECCFMHn/WLd2GHf/J/DtfnjT+wtt8ltjD25GroxSOm2156K3UiFyUZuTYXkQp6dLF/PZJ\nJ+Wd3KfP0zlxca8LIY6+KTURdSgGXaIgF2V3DUjMKfhq8hX3nTXrlmdSnYlpWpe0F5M1GjHJBcqP\nv2ldSeczGbhZRKQxGwy4eNgw9/OzZo09Nifn4wyX6zohBHspEwUpBl2iICWEiPamdnm413FzXjjt\nwQ+6pRcWBe0v08FzLrQuvsUVcYvSTEYFTU1NWpdBGkiy2/HYtGn+28ePP7dnYuI3cTExx2hdExHt\nj0GXKAjFuuPHpuT3/XL6Px6feey5N8UbzRatSzqgpNye4vfSFBFpmc9pAyq2R+T0ZAoYkp5ueuuk\nk3LP6N//UcnrfUUI0f67sxDREWPQJQoiQghnXHrOC/1POP2RU//1vpQo9ejUdmFHo1vJLCz7ty6i\ndgrzuaGUbqrVugzSmFGvxzlFRc6nZ8wYNyg19bNku/0ErWsiIhWDLlGQsMcnz0ovLPp8zm0vTBix\n8FK33hBaK516T1hguefV1Dqt6+hMPi9E2Z8RuQ6PWpHmcuGVuXMzT+vf/54u6uiuU+uaiCIdgy6R\nxoQQvviMrm8Xzz//zpPvezvT48/SuqQjYo6OhdWXh1/Xal1J50nwwLhjK4Mu/UWn0+HsgQObR3c/\nT7TZpmpdE1EkY9Al0ogQQriS0hdLg8Z8vODu10cOOPEMh04X2v9LDp5zgXnxLc4arevoLHFuGCq3\n13GSLu0nvXl0t1+/eyWv92WO7hJpI7R/qxKFKCFEYnxm1w9Lzrzq73PveCnNHpekdUntIiWvr27V\nn8mIlEVpLjtQXV7LoEut0ul0WDx4sPPJ6dOPG+D3f55ks03RuiaiSMOgS9TJ7PFJEzL7Dltx0j1v\nDioYOz1Gi617O4oQArkjpuHuJ0VELEpz2YHayvoIifV0pDLcbrw2b17mQnV090UhhE3rmogiBYMu\nUScRQhg9qdn39R4/7/4Fd/0nPdYdr3VJHaLPpFMsd7zgj4hFaU4bUFcVEZmejpJOp8OSwYNdj0+b\nNr5XUtJndoult9Y1EUUCBl2iTiCE8Mdn5X086bJ75o4640pPqM/FPRBrrAMmTzes3aB1JR3PHgs0\n1DaGz5A8dbgsj0e8Om+edGxOzkvpLtfFIpze0iEKQuH725YoSDgS/NO6DCj54JR73+qV0WtI0O5u\n1p4Gzj7PtPhme9jvlKbTAcaQ6XRMwcJqNOLuiRMTLy4uvriLx/MfTmUg6jgMukQdRAhh8qZ2ebjv\npJP/OffOl1OjnR6tS+o0aQWD9N9uiIyd0gw6RdG6BgpNJ/boEfvE9OklhYmJnMpA1EEYdIk6gBAi\nw5ed/+nUq5dPH37yxe5wnqrQGiEEugydjAdfRNhHXZOBOZeOXKbbLV6fP186Nifn5QyX6xJOZSBq\nX5H125eoEzgT/HNzhhz73sL7/1vg797PqHU9Wuk35TTLzU/6w35/XGNobWBHQciiTmVIuHjYsIuy\nPZ43hBB2rWsiChcMukTtRAhh8aZJTwyYduayOcue91ttkd0fPtrhhs6Rg41btK6kY5mNCkfgqF2c\n0L177JPTp48KdGXoo3U9ROGAQZeoHQghEnxZeSunXffvKUPmLHXy3UfVwFnnmZbcHBvWO6VFWYDq\n3dwzgtpHptstXps3r8vI7OwX/Q7HfK3rIQp1DLpER8kcHdsztceADxbc/Z+eSbmFfCO7hYzeQ/Vf\nrEnRuowOFeeGUropItoGUyexGI341+TJiTN79rwx0+1exnm7REeOQZfoKDgS/FOy+498+aS738gK\n1w0gjoYQAhkDxiuPvoKwXbHl8wClmxl0qX0JIXBRcbHnypEjT8lyu18VQli1rokoFDHoEh0hd0rW\n5d1Lpv7fjBufTDZa+DuoLQOmnWW9/tGUsJ2+4HNDX76lXusyKEyN79o1+sGpU0d3jYv7SAiRoHU9\nRKGGQZfoMAkhDN7ULo8VL7hg6TFLboiLtNZhhyvGFYem2C7YWqp1JR3D54WxYltd2I5Yk/byfT79\nC7Nn9+yXkvIB++0SHR7+hiY6DEIImzc9592Jl909pc+E+WwBdIgGTD/HuPSWmLAc1fU4oNtdVssh\nXepQ8bGxeGnOnKzhmZkvpjqds7SuJ1hIkrRWkqTyQzx2gyRJazqhpjRJkpokSXqho69FB8egS3SI\nhBCpiVLByjnLnh8cKVv5tpes/iMNH8nhuSjNZQdqK+rZdoE6nMVoxINTpyad0L37sky3+2YuUgMA\nLANw/WEc35nvvvCdniDAoEt0CKKdnkEZfYrfO/neN7t5/Fn85XKYdDodUvuMVZ57S+tK2p/TBtTu\nYs6lziGEwN+GD/dcNnz4qdlu9wtCiIjdlAYAZFm+XZblG7Wug4IXgy7RQbiS0hZIA0c/veCfr6VH\n+iYQR2PgjMXWqx5Mrta6jvbmtAP11Y3844c61cS8vJg7xo8fm+12vyGEsGhdD1GwYs9PogPwpnW5\nvtf4uScPP+VSD98lPDo2bwLqLJmibMcGuBxaV9N+LGZAKHyHkjrfgNRU0/KpU4ee/Oyz7wkhRiuK\nUql1TQAgSdL7AFIBnAngLgBeAK/KsnyiJEmFAC4HUATACkAGcI8sy/fuc454ANcBGAogCUAZgHcA\nXCXL8uoWx60FYJdl2dniPheAqwBMBOAE8CmAc1upcx6A5QDOkWX59la+hiEAHLIsVwbuMwI4A8B0\nADmB+jcD+A+Ay2RZ3n443yfqHBzRJWqFEEJ4/Nl3Fs1eevqIhX9jyG0n/U5coj9/WVSt1nW0N6Oe\nSZe0ke/z6R+fPr1f17i4/wkhPFrXE6AAcAN4EsAHAB4E8IEkSWMBrARQDOAlAHdAzSF3S5K0J+hK\nkmSBGh5nA/gcwK0AVkANmCslSdr3rTWlxXNjAHwINWT/DOAeADYA7wFoawFxW///7nv/E1DnBNcC\nuDdw7hoApwJ4vY1zkMY4oku0DyGEcPuz7h1+8sUnFo6bbdO6nnAiDR5jvP9+f436+yd8GA3MuaSd\nTLdbPDtrVsH0J574QAhRoijKBo1LEgBiANwiy/L5ACBJUhSAtQDKAfSTZfmPwP0XA3gKwCmSJL0o\ny/J/AIwEUAB19Paq5pNKknQugBsBTANwdxvXPh9ALoBrZFm+IvA8HYCHAMwCsO1IviBJkvoDmATg\nUVmW57S4Xw/gKwC9JUnKlmX51yM5P3UcjugStRAIuQ+VnHHVdIbc9qfT65Hcc6Ty2vtaV9K+TBwy\nII0l2mx4ftas3D7Jye+a9PpsresJeK7Fv8cD8AC4qTnkAoAsywqASwKfzg/cNmeTHpIkmVuc4y4A\nflmW2wq5gDrquwPANS2u0QRgKYCmI/kiAtYDmAt12sUesiw3Avgo8Kn3KM5PHYQvz0QBQgid25/1\n+DFLbhjXtXh8tNb1hKvBs5ZaL730+epjizeFzXZyJoPCuS2kOXd0NJ6fPTt7xhNPvOmwWqfuqK7+\nUsNyFAAte9b2Ctz2liTpylaOb4I6igsAbwP4HcAEAFskSfov1KkMr8qyvLGtC0qSZAWQBeB9WZb3\naoUiy/J2SZJ+AXBEC/cC131EkiRDYJ6xBCATQE8AIwKH6Y/k3NSxGHSJAAgh9B5/9jPjzr91jDRo\ndNgEsGBkj09GlT5dVO7aBFuM1tW0D5MBaGhogsHAN8lIW7FmM56ZOTN97tNPvxAfEzNvy65d72pY\nTssuK81LUKe1cawCdeEYZFmuDkwVuBTAVKhTBiYBaJIk6XkAp8qy3NomEc1zd3e2cY0yAImHXv7e\nJEk6FepSZNAZAAAgAElEQVSIbvNWzOUAPgGwCkA/qFM2KMjwVZkinhDC4EnNfmnCJf8cy5DbOfqe\nsEh/4W3WsFmU5nZA2bGVvXQpOFiMRjw2bVrK4PT0R5Lt9uO0ridgV+B2uCzLulY+9LIs73nrX5bl\n7bIsnyPLcjLUkd4LoQbKKWh7fm5z+G1r0dm+79Q1T65vLQtFtfxEkqSpgetuBXA8gBRZlt2yLB8L\n4Js2rkdBgEGXIpoQwuRJ7fLa5MvvHZXZdxh7UXaS3KHHGd/82h82K7ji3UDZ5rDJ7RQGDHo9Hpg0\nKbFvSso9vtjYUVrXA+DbwG2ffR+QJMkuSdItkiTNDHxeIknSnZIkZQCALMvfybJ8E4C+AHYDGNza\nBWRZroYahgsDnRtaXiMWakuwluoCt9H7HCsAZOxz7IzmW1mWX9lnCkVu4JYjukGIQZcilhDC4k3r\n8sbUq5cPS+s5mFv6diK9wYCEvGLlnY+1rqR9JHihK9tSr3UZRHvR6XS4f9KkxN7Jycu9MTFDNC7n\nBQCVAC6UJGnfxXI3AzgH6pxXQO2beyb2733rg9q7dt0BrvMQ1I4PNzTfEQiu1wIw73PsqsDtMYHO\nDM1OB+Da59iaFjXsIUnSHKj9dhUAEb1LXbDiHF2KSEKIKG+a9MaJ1/57QFJuT/5/oIHBc861XnDl\ny1VfDtgcdfCjg5vPA+OKdXUHP5Cokxn0ejw0dWryzCeffNQVFTWlrKrqs0669F6jm7IsV0iSdDKA\nxwF8LUnSC1A3WxgKdZT3M6iBF1D71Z4N4HRJkvKhzoO1QZ220Ih9Oh/s4zaoHR4WSZLUN/Dc/gC6\nB663py5Zlr+RJOlLAAMArJAk6YPAccOgbjLRr8V5HwFwIoAXJEl6Auo84L4ABgL4L9SWaMHSx5ha\n4IguRRwhhDUuPee/M258YiBDrnaciWmoVNJ0u3ZrXcnRi3NBv6u0tlHrOohaY9Tr8ciJJ6YUJCQ8\nE2s29+iESypoZRMGWZafhTr6+Q6AsQDOgjr6ejWAkbIsVwWOqwEwCsAtAOKgju5OgbrZxBBZlt/Z\n51otr1EPoATqiG4CgNOgZp3RAP5opa5xAB4GkB2oxwo16H7S8lhZll+HupBuNdR+vPMBbIcadi8I\nHDb24N8a6mxC4YY+FEGEEAZvmvTm9H88PjQhO5+tYDT27RtP1cX/cgpuu7AmpKeOfPUjsPiB2Oo5\n12YfdDHj7VP/qH5/5ulWo54/ftS5quvrMfWxx9Z+uXHjMVV1dasO/gyi0McRXYoYgW19n5hw6f8N\nYsgNDt1GTDS9/FnK0TRxDwouB1C7s4GjBhTUrEYjnpw+Pa1HQsKrJoMh8+DPIAp9DLoUMdz+rDtL\nzrp6TEavIfsuSCCNGIwmxOUMblqhZVv7duC0AXVVDVxxTUEvxmzGUzNmZBQkJLwhhPBrXQ9RR2PQ\npYjg8WddOnDamTPzR04Oky0KwsfgOedHnXtHfJXWdRwNWwzQWBfyA9MUIewWC56eMSOrZ2LiW0KI\nI95AgSgUMOhS2HMlpc/PGzl58cBpZzoOfjR1No8/C6UNabqamoMfG6yEAIw6zlyg0OGMisIzM2dK\n+T7fm0IIvjZS2GLQpbBmi0sck9F76LWjz7rGe/CjSSuFExbqrrjbFNKNaI2G/VeZEwUzT3Q0Hjnh\nhDzJ631dCBHSC0KJ2sKgS2HLanP0Su7a675Jl92TIASnTwaz7iUnmJ5d4Q/pPXSNBuZcCj1pLhfu\nnjChT5bb/bwQgpmAwg5/qCksCSEyErr0eGb6DY+n6NjGKegZTGa4svrj8++1ruTImdiRmUJUYVKS\n4ZqSkmGZLtfdWtdC1N4YdCnsCCHikrv1fn3Orc+lG82Wgz+BgkLRnPOtS5Z5Q3ZRmtmo8G0DCllj\nJSnqtP79p2a4XJdoXQtRe2LQpbAihIj1Zee/OfvW5yRLrF3rcugwxGXk4s9qv64uRHfSjbYCuytD\nevYFRbhT+vZ1jpGkxYk22wStayFqLwy6FDaEEAZves5rM258oofNm6B1OXQEehx3srj2fmNILkqL\nd0PZvilEUzpRwLUlJXE9EhLuiDWbu2tdC1F7YNClsOFJzb73+Ivu6OtN7cK3kENUz7EzzI+9G5qL\n0nweiLLNtVqXQXRUdDodlk+ZkpIbF/esECJe63qIjhaDLoUFd3LG6X0nnTwxs08xdz0LYUaLFfa0\n3vhmldaVHNyWUj16Tc/Ew6+oLUh9Huh3bFUHo797px73nrkbd560Cy/eUo1dZftvJrGmogL9n3gC\nKzZu7NS6iQ4mymTCIyeemJ0bF/e6EMKqdT1ER4NBl0JelN3dL7XHgL8NnrXEqXUtdPSK5l5gXXyr\np1rrOg5kd7XA2TckYne1Ds1vH/g8MO7YVtf05++NeOtftYh1C+QVG7Hh50a8tGz/3TDu/e47dHO7\nMTgpqXOLJzoEvthY3DdpUo8st/spwf6MFMIYdCmkCSHi49KlRyddfm8iX4vDgy87Hxt2+kVDkE5g\n2LjVgNmXpuC7X/fu6OF2QFSX1dX/8H4DrNHACZdZMWy2GSMXmLF5dRO2rGncc6xcVob3N2zAqd05\nDZKCV77Pp18yePDQNKfzQq1rITpSDLoUsoQQpvjMrq/MuPHJLIOJMxbCSf7YubjxQf0RRd3bHnMj\nZ0J2qx9Lb/btOe7pt+wYuiAdvWdk4ozrErGldP9+y6vXm5A7MRvvfxENAHjoZQeOW5SKX9aZ0T9/\n705oLjtQU1nfVLGtCc4EHQxG9Q8vr189b8W2vzaUuO+HH9AzLg59fT4QBbNZPXvaBqWlLXJarQO0\nroXoSLDFOYUsT2qXBydddm8BOyyEn8Lj5lgeXHh79SWn/H7Yr1E/rzHDZFSwcHLZfo91SVW7Ivzw\nmxmX3xWHgi41KMipxgvv2nH2DYl4+qb1ex1/+2Nu9MiuQXHv3QCAf7/iRHJ8A64+YwvWbDThk++j\n9hzrtAG1uxoUS5rAztK/Qm1dtfpvc+DQJhh0n27ejPtHjTrcL41IE7cee2zCr9u2/VsIMUBRlO1a\n10N0OBh0KSS5ktNPHTT97LH+7v2MWtdC7c9kjUZ0UgF+Wv07umYe3nN/WWdGVkodzpq2f9Bt9tx/\n7bDHNOHf126Ayaige3YNlt6SgB9Xm9EtU+2c8ONqM97+NAbLr9yw53nXnLkFA3tUQQjg9w2mvc7p\ntAH11Y0iIVuHVSsb8OvnDfB30+PL1+thNANxqerIbpPeauyXkIDuXu/hfWFEGjEbDFg+dWrW5Ecf\nfVEIMVRRlMaDP4soOHDqAoUcoyUqP7lr78sGTj+Li8/CWNHcCy2Lb3Yd1qK0XVU6bNpmgJR24DZf\nG7YakZ5YB5NRHW2V0tXjN2796++m2x7zoE+3agzo8VcJgwrUkNsakwnQQUH34UbEp+vw0rIa3Hny\nbvzyWQOGTDfDGiuw9rsGKDDoFubl7Xlek6K0fkKiIJJst+O60aN7Zbhcd2hdC9Hh4IguhRQhhC2p\na6+np1z5QBIXn4W3pNye4rVSv2hqKoPuEP8k/3mtOsp6sKBri27En9v/evnbXaWOtsZGqQNVX66y\n4MOvovD49etbfX5bTAZFMRgFZlxlxa+fN2D3DgUpuXrEpannX/F0HYRS3yi5XPoft2/HNZ9+irWV\nlUiNjcWFffqgMJ5tSyl4jcjKskzOyzsh2W5/b0NFxbNa10N0KDiiSyFDCCG8adKz0677t2SyRh38\nCRTyupXMxG2P6PZvQtsGea26KLGsQo/5lyehz4xM9J2ZiUX/SMCajX+N1vbMqcGvf5jw30+isatK\nh4dedsBqVpATGNm97VEPigqrUJi7f1uwAzHq1dFZvUEgZ4ARvcaa9oTc375owJbfm6Bvqq5vaGrC\nhStWwG2x4PbiYmQ6HLjgww+xM1T3P6aIcXFxsadrXNzNQogMrWshOhQMuhQyXMmZV486/YoBHn82\nh3IjRO8JCyz3vOI/5O3GmoPu8hddiI1uwomjK9C9Sw3e+jgGJ5zvx89r1BHfE0oq0C2zFmfdkIje\nMzLx5spYnDd3G5y2Jnz0TRQ+/9GKJTP/WnPTdIhR29jGe2SKomDF03XIGWiAQKPy8ebN2FZVhaW9\neqFfQgIu7tsXu+vr8cbatYf6pRJpQqfT4f7Jk1Pz4uNf4GYSFAo4dYFCgjkqtrDH6BNOyh81JUbr\nWqjzmKNjYfHl47d1a5GVevDjDXogKa4eNyzegj7d/ppb+8r/YnH+Mh8uudOH52/9Ayajgif/8Qfe\n/iQG28oN6NutGrkZzaO5bozqvwvdMmvx3S8WXHJnPFZvMCE9sQ5Xnr4VffPanjZsMiit/hH288oG\nlG1qwvFLLfjtf8CGnTsBAP7YWABArMkEh9mMjbt2HfL3hkgrDqsVd0+c2G3Bs8/+G8BUreshOhCO\n6FLQE0JYPP6sf487/1b2EYtAg+dcYF50s/OQFqVdfupWvHPf2r1CLgAcN3Qn+nSrxk+/m/dMYTAa\ngGMG78Lc43bsCbn//TQaP6y2YNGMUtQ3AGffkACPowEPXLERXVLrcNb1iajc1fbLptkINNTvPfzb\n1Khg5bN16DbEAKdPfW5jYAFaY4uFaHWNjeBbFRQq8n0+/ayCghFJdvuJWtdCdCAMuhT0PP7suyZc\n+n+clxuhUvL66lb9mSwOdfpAW3Jb6azQkqIAtz/mwTGDdyLbX4cPvozGljIDLj1lGwYVVOHqM7Zg\nd7UOr/wvts1reJ1QyrbsPc/2xw8aUFmqYMCkv9qRpdpsAIDvtqvTI9ZWVGBXfT38gfuJQsHZAwc6\nM12u64UQcVrXQtQWBl0KajGuuFFdhx1/bEpeH06ziVBCCOSOmIZ7nhIH7N3Z2Kj2vv3u19Z3yaup\nU8dLzabW23m99mEsft9gwqIZpQCAdZvVYJqaoAZXW0wTnLZGrN/SduvmeA9QuumvoNvYoGDl83XI\nH2aEzfPXy21fnw9xUVG48uOPsezLL7H0f/+Dw2zG6NRDmJ9BFCR0Oh3umTgxXfJ4nhJsg0NBikGX\ngpYQwuZKzvhnyRlXcbQgwvWZdIrljhdSD9iSoKFR4ITz/TjlquT9Fo8pCvD1z1YY9Mqekd29nwvc\n+YQbE4ZVIjWhfs99ANDY9Nfv79o6saePrhDYr6dugge68hYjut/+tx7VOxUMmLh3ODbr9bituBgJ\n0dF47tdfEWM04pahQxFl5P4nFFqS7HacU1TUK9XhWKx1LUSt4SgZBS1vmvTQlCvvz9bzl3/Es8Y6\nYHR3xbqNa5Ga1PoxZpOC4t678c5n0bjvORdOm/rXzmjLX3Ti1z9MmDisEjFR+8+BeOFdGzZtM+DM\naaV77stIVgPrNz9bMKBHNVZvMGJnlQ5piWoQnji8EhOHV+51Hp8Xxh2/1e/5vHCMCYVj9t5BrVmm\nw4F/lZQc2jeAKIhN69Ej9tVVq5YKIV5WFOV3reshaolBl4KSPT75xIHTziyOS8/h22EEABg4+zzT\n4ps/qn5xWUWbLY0uWrANX/9swW2PufHZD1ZIabX44TcLPv/Rimx/HS46adt+z6mrB+56yo2poyqR\n6G3Yc39Rzyr43A04f1kCji3aiXc/i4bT1ohxQyr3O0czrxP6naV19QD41xlFlDvHj0/5vazsKSFE\nf24RTMGEUxco6AghvN7U7OuK5izlFr+0R1rBIP2361MOuCgtxVeP5275AxOHVeKXdWY88qoDm7cb\nsGBCOZ644Q/YY/Z/8lNvOlBeqcfpJ5Tudb/ZpOC+yzciOa4ej//HjtjoJtxz6SZEW9vestdlB2p2\n1DW0eQBRmHJGReGakpK8dJfrOq1rIWqJI7oUVIQQwpue8+SUq/6VoTvUfV8pIgghkD10Mh568Yem\nBZPa/iM9wduA6xdvOeTzzh63A7PH7Wj1sS6pdXjyxkPfBthlB2p31redhInC2IisLEtRWtrsGLP5\nqV21tV9pXQ8RwBFdCjLOpPSzhsxZ2tvhS9G6FApC/aacZrnpyUPfKa2zOe1AXVUjp9tQxLphzJiE\nLLf7YSGERetaiAAGXQoiQohUX1beeb3Gz2UzUWpVtMMNnSMHm7ZqXUnrYqKApvqjbPhLFMKiTCYs\nGzcuN9Ptvl3rWogABl0KEkIIXXxm16cmX36vn+0Y6UAGzjzXvOQm2yHtlNbZhACMes5coMhWmJSk\nH5yaepzVaMzTuhYiBl0KCs7EtLOHn3JpfrTTo3UpFOQy+hTrPv89OWj/GjLqwaRLEe+akpKETLd7\nOTeSIK0x6JLmhBBOV3L6ovxRU7jHLx2UEAIZA8crj70anIHSZAjKsog6lc1iwdKioq5+h+MsrWuh\nyMagS5rzpnW56/iL7kjnH/50qAZMO8t6/aMpQTl9wWQAf5CJAEzq1i060+U6Rwjh1roWilwMuqQp\no9lakN1/VLE3TWI4oEMW44pDY3S22Fp68GM7m8mooOlAzX6JIoQQAsvGjUvv4vE8oHUtFLkYdEkz\nQgjhSc2+v+Ssa3xa10KhZ8CMc4xLb42p0bqOfdmigd0VDLpEAJDmcmF8166DPdHRw7SuhSITgy5p\nxuHznzZ03nm55qgYrUuhEJTVf5Tho5+Dr99yvBvK9k1B2+qXqNNdMGSIJ9luv1MIwa2xqdMx6JIm\nhBB2V3L6uT3GTIvWuhYKTTqdDql9xiovvK11JXvzeSDKNtdpXQZR0DAZDPh7SUl2usv1d61rocjD\noEua8KZJdx5/0R0ZXIBGR2PgjMXWK5YnB9WitAQv9OVbGXSJWipKTzcVJCTMEEKkaV0LRRYGXep0\nRou1e2afYaPiMnKZcumo2LwJqDVnoGyH1pX8Jd4NU+V2bo9GtK+bjjkmOTcu7l9a10GRhUGXOpUQ\nQnhSsh4Ys+haLkCjdtFv2hLj+bdFBc2kWJcdqC6vrde6DqJg44mOxoSuXXu4oqKGal0LRQ4GXepU\njgT/yUWzz8k1R8dqXQqFiZzBYw3vfu8Pml0aXHagpqK+Ues6iILR4kGD3Ml2+63cMY06C4MudRoh\nRKwzIfWCnuNmsc0CtRudXo/kgpHK6//TuhKV0w7U7W7QugyioGQxGrGwb99sX2zsDK1rocjAoEud\nxpOafcv4C2/jAjRqd4NmnWO99L7EoFiU5rQB9TWN/CEnasOMgoLYRJvtUrYbo87AoEudQggR58vK\nK/Fl5/Nnjtqdw5eC3YYMVO7SuhLAYAD0CJqZFERBR6/T4ZLi4oxUh+NcrWuh8MfQQZ3CmybdPPrs\na1O1roPCV98pZxkuus0aFH29TAYmXaIDGZmdbU51Ok8WQnDBBnUoBl3qcEKIxASpR7HHn6V1KRTG\ncovHG9/4OiUo2noZDMy5RAdzTUlJWqbbfYPWdVB4Y9ClDudNz7llzNl/D769Wims6A0G+PKKlXc+\n1roSwKTXugKi4NcjIUGfFx8/TgiRoHUtFL4YdKlDCSH8yV17DXYmpmldCkWAwbPPtV54V0KV1nWY\njAoXoxEdgmtKSvxdPJ47ta6DwheDLnUob3rOstFnXZOsdR0UGVxJ6ahQ0nRVGkddiwmivjYoZlEQ\nBTW/w4Gi9PSBQogcrWuh8MSgSx1GCJHh796/vz2eOZc6T+/JZ+gu/T+LpovSvE40lf4ZFOviiILe\nJcXFCble761a10HhiUGXOkxcRu5tJWdclah1HRRZuo2YaHrpE20Xpfk8QOmmoNmVmCiouaOj0dfv\n7yGEyNC6Fgo/DLrUIYQQXdILi3rbvFxjQJ3LYDTBKw1SVnypXQ0JXujLt9RrVwBRiLlw6NBEyeu9\nUes6KPww6FKHiMvoevvI069gyiVNFM0933renfGa7ZTm88BYvo1TF4gOVaLNhu4+X192YKD2xqBL\n7U5vMHbL6je8Z4zTq3UpFKE8/mxsr08VNTUaXd8JXVVpHYd0iQ7DJcOGpWR7POyrS+3KoHUBFH68\nadLtIxb+LV7rOiiyFU44VXflPd/U37CkztjZ13bZgZoddQ0AOv3akWZ7dTXu//57fLRpE8pramAz\nmdDX58PC7t2RFBMDADj+pZfw5+7dBzzPZf37Y1yGOkX0gw0bcPtXX6G0pgZd3W5c0Ls30uz2vY7f\nUVODiS+/jDMKCjC1S5eO+eIiTLrLhS4ezxAhhFtRlFKt66HwwBFdaldCiJysfiPyouwurUuhCNe9\n5ATTMx/6G7S4tssO1O5q4PZoHWx7dTXmv/kmXvztN2TY7ZgmSejmduPNdesw7803sX7nTgDAdEnC\nKfn5+33M7toVOiFgNRjQze0GAGyrqsIlK1bApNdjUnY21u/cicXvv4/axsa9rv3wTz/BZjZjYhZ3\nfGxPlw4blprpcl2jdR0UPjiiS+0qLj3370PmnsvRXNKcwWSGK6uf8vn3v6FPfude22kD6qo0ydgR\n5f7vv8fWqiosKSzE9Jy/2rC+sWYNrvj4Y9z+1Ve4eehQTMtpvUXrzV98gSZFwbm9eiE9MGL7xtq1\naFQU3Dl8ODxWK4anpOCkt97Cyk2bMCxF3eBxW1UVnvv1V5zbuzcMOo4Xtaeu8fEi3eUaLYSIVRRl\np9b1UOjj/6HUboQQDm9al17stEDBomjOBVHnLPN2+vYRUVZAaeDuaB3tfxs2wGmx7BVyAWBMejqS\nYmLw6Z9/tvncb7ZuxTO//IJ+CQk4LjNzz/2bdu+Gw2yGx2oFAGQ7ner9u3btOWb5jz/CGxWF4zLY\nDasjXFRcnJrudF6qdR0UHhh0qd24kjPOHzr/Ar/WdRA1i8vIxeZqv66ukxsgCAEYDZy50JGaFAXz\nu3XDKfmtD9cb9XrUNzWhoan1lsq3f/01DDodzuvVa6/7Y00m1DT8NRq/u15dUxhjMgEANu7ahZdX\nr8Yp+fnQCf4t0xF6JyfrUxyOSUIIq9a1UOhj0KV2IYTQx7p9k1Ly+vBnioJKj+NOFtc9YOj0eQRG\nPZh0O5BOCJwoSZicnb3fY2srKrCushJJMTGtTi14b/16/FRainEZGfDbbHs9lu/xoKqhAU/JMqrq\n6/HoqlUQQiDf4wEAPPD99/DHxmJMWlqHfF2kOreoKC3Fbj9b6zoo9DGUULuIdnqn9Ju6kKO5FHR6\njp1hfvSd1E5v9WXiiK4mmhQFN33xBRRFaXOh2OM//wydEJidm7vfY0VJSRiclIRbv/wSw555Bo+t\nWoVZubnIsNuxtqICb6xdi4Xdu+91PWp/Q9LTjd6YmFlCcNicjg4Xo1G7sMUlLu0+amqU1nUQ7cto\nscKW1gvfrFqNgv1zTcdd1wD+gu5kiqLg+s8+wxdbtqCr241pkrTfMXJZGb7btg3DUlKQHBvb6nlu\nGToUH27YgPW7dqGb240eXrUn+L3ff49spxPDUlKwfudOXLFyJX4qK0Oc1YrFhYUY4eff+u1FCIHJ\neXkp8rZtQwD8T+t6KHRxRJeOmt5g7C4NGp2uN7JlKAWnIXMvtC5Z5u7UndLMRgVNbcwPpfbX0NSE\naz79FC+vXo3kmBjcNGRIq9MWXluzBgAO2hasKDkZM3Jy9oTcX8rL8d769TgtMJp7+cqVqG1sxLKh\nQ1GUnIzLVq7EH5WV7fxVRbbZPXs6Uh0OLkqjo8KgS0fNm9blmkEzF3MbNApavux8rK9MFQ2dOFPX\nHgPsLGOLsc5Q09CA8z/4AK/9/jv8NhvuGjlyT9eEfa3YuBF2sxl9fL7Dusbd336LfI8HAxMT8Wt5\nOX4qLcXpPXpgQGIilhYWwmYy4cXVq9vjy6EAm8WCfJ8vRwjBlpV0xBh06agIITxxGV0LuN0vBbu8\nMXNx00P6Tkue8W4opZs6ud1DBKqsq8MZ77yDlZs2QXK5cP/IkYiPan0W1brKSmzctQtDkpMPq2PC\nd9u2YeWmTTg1MJrbvBFFSmDqg16nQ2J0NDa2aEFG7WPRoEEpGS7XRVrXQaGLQZeOijsl66LiBRek\naF0H0cH0Gj/HsvyNzluU5vNAlP7JoNuRahsbsfT99/FjaSkK4+Nxz4gRcFgsbR7//fbtAIAC7+H9\nYX73d9+hd3w8eserA4uNgQVojS0WotU2NnJSdgfoFh+P+JiYsUIIrimiI8KgS0dMCGGM9frGJUoF\nfH2noGeyRiM6qQCrOund5cQ4GMq3dnqzh4hy1zff4Pvt29Hd68XtxcWIOsg6gV/KywEAUmATiEPx\n2Z9/4qstW3Bajx577ksNtCRrDs6VdXX4Y+dO+NtY3EZHZ17v3inuqKiJWtdBoYl/IdERi/X4pg84\n4YxUresgOlSD51xgWXTz+9Vv313W4Y3o41ww7vyxrhGAvqOvFYm2V1fj2V9/BaAGz4d+/HG/Y4QQ\nmNu1K0x69T/BhsCUA28b83dbc8+332JgYuKeProA0MXpRI7Lhdu++gq/lZfj623bAAATW+npS0dv\nYteuUXd+9NFSAM9oXQuFHgZdOmKxnoSzug2f0Pb7hERBJrlroXi9NEU0NZWhlQX57cplB6rL6+rB\noNshfti+HQ1NTRBC4JU2FoEJITA9J2dP0K2sq4MQYs8uZwfz4YYN+KmsDP8eM2a/x24cMgTXfvop\nXvjtNyTExOD6wYOREB195F8QtclkMGBwWlqaECJbUZRfta6HQguDLh0RIUTagGlnpugN/BGi0NK1\nZBZuf+T7xnPmNnVoAHXZgdqKOvYX6yDFKSn4dMaMw3rOAyUlh3V8UXIyPpk+vdXH4qOicMewYYd1\nPjpyZw4Y4Ht39eorAMzSuhYKLZyjS0fE489e3G/KKYfXn4coCPQ+fr7l7ldTO3yVmNMG1O5u7OjL\nEEWEFIcDyXZ7fyGEWetaKLQw6NJhE0IIq901Mj6jq9alEB02S4wNlvg8rF7Xsddx2ICG2ka+xhK1\nk1k9eyY5rdbxWtdBoYUvwnQkCnKKxiZqXQTRkRo0+wLzopudHbpTml4PGHQt+k8R0VEZK0mWRJvt\nNK3roNDCoEuHLS49d2mv8XNdWtdBdKT8+X11P/2ZLDp6h16jAQy6RO3EajQi0+XKFELYtK6FQgeD\nLqXLEoEAACAASURBVB0WIYQ+xh3Xxx6XpHUpREdMCIGcYSfi3qdFh06iNeqZc4na07xevZISYmNn\nal0HhQ4GXTos5qjY4d1LTmDKpZDXd8pCy+3P+Tt0UZqJTUmI2tWQ9HSDNzp6ttZ1UOhg0KXD4kjw\nL+kx+oQYresgOlrWWAcM7q5Yt6njrmEyKNw1kKgdGfR65Pt8KUKIeK1rodDAoEuHTAhhcSSk5Fhi\n7VqXQtQuBs0+z7TkZltNR50/ygLUVrOVLlF7WtC7d5Lf4ViodR0UGhh06ZBFu7xTeh8/n9MWKGyk\n9Rys/+aPFHTUojSPC0rp5tqOOTlRhCpMShLuqKgJWtdBoYFBlw6ZzZNwck7RMWzWTWFDCIHsoknK\nwy+hQ6JuggeidFOH701BFFGEEBjg9ycJITK1roWCH4MuHRIhhCsuPTfdYGLOpfDS74QzrDc94e+Q\nYVefB7qyLQy6RO1tXq9e8Zlu9yKt66Dgx6BLh8ThS5nfd/LJnLZAYSfa4YZw5GDT1vY/t88NY8W2\nOvYYI2pnXbxeuKOihmldBwU/Bl06JFEOz5S0wiK91nUQdYSBM881nXNLbLsvSvM4odtdVl/f3ucl\nImCg3x8vhOiidR0U3Bh06aCEENGu5PQEnY4/LhSeMvoU6z9bndLu53XZgZoddR26KQVRpJqcnx+X\n7nTO1boOCm5MLnRQlhj7qLwRk9izkMKWEALp/Y9THn+tfbfsdTmAul0NnLpA1AHy4uNht1pHaV0H\nBTcGXTooW1zSnC4DSyxa10HUkQZOP9t63SMp1e15TqcNqKtq4KYRRB1ACIEcr9cnhHBpXQsFLwZd\nOiAhhIh2uCVrrEPrUog6VIwrDo3R2bqtpe13TqsFQBMHdIk6ytT8/ERfbOxkreug4MWgSweTndpj\ngFfrIog6w4AZ5xjOvTWmXRelmQwMukQdZXBamt4dFTVN6zooeDHo0gE5k9Kn542czKBLESGr/yjD\nip/bd1GaUa8w6RJ1ELPBAF9MjF8IYdC6FgpODLp0QNZYx+gEqYfWZRB1Cp1Oh9Q+Y5QX/9t+5zTy\n1y9RhxqVnR2nF6KP1nVQcGLQpTYJIazOxFS2FaOIMnDGEusVy5Oq2ut8JgO4GI2oAx2Tk2PLcLvn\na10HBScmGGqTyRo9Im/EJJ/WdRB1Jps3ATXmTFFe0T7nMxsVNDU1tc/JiGg/focDDoull9Z1UHBi\n0KU2ORL8c9lWjCJRvxOXGM9fFlXbHudy2ICKbQ3tcSoiaoPk9cYLIbiehPbDoEttinJ4cqLsbE9I\nkSdn8FjDO9+ntMsiMp8bSunmdsnMRNSG8bm5PldU1Bit66Dgw6BLrRJCZPnz+/KvY4pIOr0eyQUj\nlf98cPTnSvBCV7q5/uhPRERt+v/27juw6TJxA/jzzepIZ7qBsgoECgg4QMWFinfnOM95jlPPcXfq\n3amnnvOc5/x57q3nZqiIoKKAC9kUaAtlpoO26W7SNk2a0az390cLV/Zq+2Y8n39ik3yTJ5WmT9+8\n3/ednJurztTrfyc7B4UeFl3ap9QBQ68YP/1SbvtLUWvqH+6Me+CdAUe9U1pOBjQ2i7c3IhHRfqTE\nxSElLm6Y7BwUelh0aZ9iE5PPHjB6kuwYRNKkZOfCqR4Oe8fRPU6GARqH1cdJukR9bHBKSrqiKEmy\nc1BoYdGlfdKnpGdyWTGKdpMv/Zvm/pfjjmo41pAMuNs6WXSJ+tj0ESMy4rTaqbJzUGhhk6G9KIoy\nKHP4mBTZOYhkG3PGb7ULi3KPam0wQzLgsfu4vhhRH5s6dGhsbnLypbJzUGhh0aW9xCYknzJiylk8\nEY2inlqjQfbYM8TPa478MVKTAJ8r0HuhiGifBiQlITEmZpzsHBRaWHRpL0kZA87LHT+ZG5cSATjl\n2rvi7n0j54hPSktOBPydAe6ORtQPshMTMxRFiZGdg0IHiy7tJUafODIhlQO6RABgGDgMtuBQxXWE\nmwKrVIBG1StL8hLRQZw+fHg6gGNl56DQwaJLu1EURaM3ZKTLzkEUSo675BbVv16PPeKT0nT8fISo\nX5w+bFhinsFwkewcFDpYdGlP+UMmnJwqOwRRKBl31sW6+WuO/KQ0rVpwSJeoH4xKT4depztJdg4K\nHSy6tJvk7Nyzhx17Kvf9JepBo9UhwzhVrCg8suO1HNEl6heKoiBdr89SFIXz4gkAiy7tITYh6awB\noyfKjkEUck697p9xd7+adUQnpem0gr90ifpJfmZmIoBc2TkoNLDo0m7ikwwDtTGxsmMQhZz0wSNh\n9Q5WPJ7DP1YfC7id3DOCqD9Mzs1Ni9VoOGJDAFh0qQdFURKTswZx2gLRfkz63V9Uj72tO+zGmpkG\n0VJ/VBusEdEhGp+drc1NTj5Tdg4KDSy6tIuiUh2fN+XMNNk5iELVMedcrvt82WDf4R6XnQ60NLDo\nEvWHISkpiNFojpGdg0IDiy7tkpY74vwhE06Kl52DKFRpY2JhyJsi1m06vONyMqBpazrsfkxER0Cl\nUiElLo6fThIAFl3qQRsbd2z64JGyYxCFtFOvuyf+Hy9mHNb2EVlp0LQ3dx7x8mREdHjS9foURVF4\nwgmx6NL/xCWmpHJFFqIDyxw+Bg3uISrvYcxESE+Bytnm5ZAuUT85YdCgZAD5snOQfCy6BKBrR7S4\nxJRk2TmIwsGE829Qnvqv5pBPSjMkA502X6AvMxHR/xw7YEBKTmIiN44gFl3aZXj2yPF62SGIwsGk\nc6+OmfHToZ+UZkgGOrm8GFG/yc/KQnJs7Omyc5B8LLoEANDGxo8dMHoiJ+8THQJtbByShh6Pku2H\ndv/UZMDnDnBeEFE/SY6NhV6nGyQ7B8nHoksAgJTsQVMzh49Ry85BFC5Ou+7euNtfTD+kndJidIAi\nRF9HIqIekmNj07gVMLHoEgBArdGNSx0wTHYMorCRPXI8zLZcxX+IMxK0ajZdov40NDU1AUC67Bwk\nF4suAQBi9IkGtUYjOwZRWBn3mz/iuQ/Vh1R1dRr2XKL+NDojIwHAUNk5SC4WXQIAxOiTkmRnIAo3\nx/322tj3Fx3aSWla/h1J1K9GpqcnJcfGjpKdg+Ri0SUoipKUkJbJFReIDpMuTg/9wInYVnEI99WA\ncwWJ+tHglBRk6vUTZecguVh0CQCGZ+WNZdElOgKnXHtv7O3PGw56UlqMVih+PzdHI+ovg5KToVap\nRsvOQXKx6BJ0cXpj1vB8bhZBdAQG5R+rlFtzleBBOqwhBUFbMzdHI+ovMRoN4rXaVNk5SC4WXUJS\n5oDjDbnD+W+B6Ajln/MHvDpTdcCdz7INQGvDYewbTERHTa/TJcjOQHKx3BBUau2Y1JwhsmMQha3j\nL7w+9vWvBx+wxWZnQNXSyKJL1J/idbpErqUb3Vh0Cbo4vUGji5EdgyhsxSYkITZrPCqq93+fnAxo\nOHWBqH/lJifHAsiQnYPkYdEl6OLi42RnIAp3U6+5J+a251P3e1JaRio0DmvnIW4vQUS9wci1dKMe\niy5Bo4uJlZ2BKNwNHj9ZtbVh/yelGZIBj83LIV2ifjQiLS0pOTZ2pOwcJA+LbpRTFEXFokt09BRF\nwehpl+OdOco+T0ozpACdDh+3RyPqR0NSUpCZkDBJdg6Sh0WXkvWpGdyziagXTL7kT7EvzR2yzzPO\nUpMAryvAk2KI+lF2YiLUijJcdg6Sh0WXMpIyB/JMNKJeEJeUCo1hDMwNe9+WlAAEvNwwgqg/xet0\nUKtUibJzkDwsupSelJ7Dk9GIesnUa+7W3fFckmfP6xUF0Ko5c4Gov8Wo1ZyeF8VYdKOcWqvLSszI\niZedgyhSDJ10irrYnIt9nZSmVYNNl6ifxWg0/NQyirHoRrkEQ+ZgfUo6/x0Q9RJFUTDi1IvFx1/v\nXWp1GvZcov6m44huVGPBiXK6OP0QfWq67BhEEeXEy26J+7/Zg/eavqDlaZ9E/U6rVnNEN4qx6EY5\nRaUaoE9h0SXqTfrUdCjJRjQ07369TiO46gJRP4vVaGIURWHfiVL8Hx/lFJU6LSYhSXYMoohz0lV3\n6f7xfOJuo7qJeiAYYNkl6k+G+Hg1AP6ii1IsulFOq4uJVRT+3iXqbXmTp6kLynN3uy7TABHwC77v\nEvWjzIQEHYA02TlIDr7hRjlNTCwn6RP1AUVRMOykC8Ts7/53Ulp2OpRAgEWXqD9l6vUxAAyyc5Ac\nfMONchodiy5RXznpir/FPflxrnvn1zkZUAf8nLpA1J8yExL0YNGNWiy6UU6t0fJsVKI+kpiWhYB+\nhGJp7fo6Kw3aYJBFl6g/GeLjNWnx8QNk5yA5WHSjnKJS898AUR868cp/aO96IcEDAGkpUBSw6BL1\np3itFnFabYrsHCQHS06UU3gmGlGfGnnSOZrl27pOSjMkAyoVf+aI+pNGpYJGpeKnl1GKRTfasegS\n9SmVSoXBx/9KfPVTd9HlTxxRv9KoVFAAnewcJAeLLvHfAFEfm3r1HXEPvzfQnZoE8E9Lov6lVauh\n5ohu1GLJiWKKoqhUajV/7RL1saSMAfDo8uB0AypFiIMfQUS9hSO60Y1FN7qpVRqt7AxEUWHK72/X\n/vPF+E61Giy6RP1Io1JBURQW3SjFohvdNBzRJeofo089V/NjSa5Q8V2XqF9xRDe68S03umnUGi2L\nLlE/UKnVGDTxbOH1yk5CFF20ajVHdKMYi250U6s5dYGo30y9+h9xPhHL912ifqRRqSA4ohu1+IYb\n3TQqtYYjukT9JCVnMBATF/zP8uXu77ZvD5ZZrfD6/bJjEUU0Tl2IbhrZAUgqtYpTF4j61TUz1qmq\nipbFFWxYFbCvL3Z5mmuVmGBQ6NUqpMXGKOOyMnFMVlbMiPR01eCUFGjVatmRicJa988QP76MUiy6\n0Y1zdIn6WXxSCvLP+C3yz/itGkB8z9s62qyoKlqGFRtX++0FGzweS70Sh6DQq1VKWmwsxmdl4Zjs\nrNgRaWlKbnIyNCzBRAfVPXWBRTdKsehGN05dIAohCanpGHfWxRh31sUa7PH+3NHajPLCpfhlw2q/\nY1Wx19vSoMRCCL1KpWTq4zA+Kwvjs7Ji87pLsJrLOxABALyBABSgU3YOkoNFN7qpWXSJwkOCIRPH\nTL8Mx0y/bK8SbLfUY3vhMvy0cbXPvny9z9faqMQpEHq1SsmMj8eE7CxlfFZW7PC0NAxKSoKKJZii\niMfngwCcsnOQHCy60a3T3+kOyg5BREcnKWMAJvz6Ckz49RVa7PERra2hBpuKlmNxyWqfY+kan6/N\nosQrEPFqlZKj1+OY7CxlXNdIMAYkJrIEU8Rx+/0QQrDoRikW3ejmcDvaArJDEFHfScnJxaTzrsKk\n867aqwS31lWiqGg5vtu42uvYuMrnb7eq4rtHggcmJopjsjLVY7OyYvLS0pCTmAhF4QdAFH48Ph+E\nEA7ZOUgOFt3o1uFx2DmiSxSlDAOHwTBwGHDBtTrssfySpboMBUXL8fWmgs6O4uWBgL1110hwblKS\nmJCTrc7PyIjJS0tDVkICSzCFLLffD18waJedg+Rg0Y1iQojgmNPO98nOQUShJ2PISGQMGQlcdENM\nz+uDwSAs1aVYWbgcX24p8HYU/uIP2tsUvVoRerVaGZycJCZmZ2tGZ2To8tLSkKHXswSTVB6fD95A\ngCO6UYpFN8oFAz5uSEpEh0ylUiFr2GhkDRsNXPqn3UaCg8EgmnZsxZKiFfhs05pO59qfAsGOdkWv\nVokEtUoZkpIsJmZna8ZkZuqGGwxIi49nCaY+5/L54Pb5OKIbpVh0o1zA7+eILhH1CpVKhZwR45Az\nYhxw+c17jQQ3lpbgx+KVmLV5Tadz1fqAcDmUBLVK6NUqZXhqqpiQk60dk5GhzTMYkBofv7+nITos\nHr8fjs7Odtk5SA4W3SgX9PtYdImoz6lUKgwYPREDRk8E8Ne9SnD99mIsKl6JGZvWeDpWrA0qHqei\nVykiQaNW8lJTxcScbO3ojAztcIMBKXFxkl4FhSOXz+f3B4Mu2TlIDhbdKBdg0SUiyVQqFQblH4dB\n+ccBuC22521+vx/12wqxoHgVPty8xuNctlqoPC7Eq1UiUaNWRqYZxMTsbJ0xI0Mz3GBAUmzsfp6F\nopXT6/UCYNGNUiy6UY5Fl4hCmUajweDxUzB4/BQA2L0Ee72o3boeXxavFK1b1rqdtSuh8XoQr1JE\nolajjEoziEk5OTEj09PVww0GJMbE7PtJKKK5vF4fALfsHCQHi26U49QFIgpXGp0OQyeejKETT1YA\n7Dafwe/1wrxpDT7fsEq0bVnrdtYtFxpfJ/QqBYlajTI6PU1MyM7WjeoeCdbrdPt5Fgp3jq6iyxHd\nKMWiG+UCfr8n4PdDreE/BSKKHBqdDsOPOw3DjzttrxLs9XhQU7Ian25YFWzbus7tql8Krd8LvVoR\nSVqtMjo9XUzMyd41Ehyn1e7nWSgcNNjtnQCssnOQHGw3ZPe6OxCXmCI7BxFRv9DFxiJv8jTkTZ6m\nwp4l2O1CVfEKzChZHWzbtN7lblii6AI+oVcrSNbplPyMDHFMdlbMyPR09bDUVMSyBIe8po4OH4AW\n2TlIDhbdaCdEW6fTwaJLRARAFxePUSefg1Enn6MCsNsaZ54OO6qKV+KjkjVB28b1LleDWYkJ+oVe\nrVJSY3QYm5mBY7KzY0akpamGpqYihp+UhQS3z+cRQnC7+yjFn8IoFwz62zpd3DCGiOhgYhOSMPrU\n32D0qb/ZqwS77DZUFS9HwcbVAXtRsdvdVINYERB6lUoxxMZgXGYmxmdnxYxIS1MNSUmBjiW433j8\nfp6IFsX4kxblvC5nvbu9TXYMIqKwFp+UgvzTL0D+6Reoscd0CFd7K3YULsPKjWsC9nVFLo+lTokV\nQZGg7irB4zMzd40ED05JgUatlvQqIpPH52PRjWIsulHObqnfZms0dwJTue4OEVEfiE82YNyZv8O4\nM3+nxh4jwR2tzdhRuBzLS1b729ds8HRa65U4ERR6jVrJiIvD+KxMjM/Kih2RlqYMSk5mCT5MQgi4\n/X6uuBDFWHTJbK0uswPIkB2EiCjaJBgyMX76JRg//RIN9vidbLc2wlS4DD9vXO13rCzydrY0KnEQ\nIkGjUjLi43BMVpYyPisrJq+7BKtVKkmvInTZOzsRCAZ5IloUY9GlWqu5zCM7BBER7S4pPRsTfnU5\nJvzq8r1KcHtTHbYULsMPJat9juXrfL7WJiVOgdCrVUqWPh4TsrKV8dlZscMNBgxMSoIqSkuw1emE\nPxg0y85B8rDoRjkhhMd4yq/5sQ4RURhJzhqISedeiUnnXqkFsNsaZ2311dhYtBwLN672dZSs8fls\nll0leEBCgpiQnaUam5UVm2cwICcxMaJLcHNHBxydnTtk5yB5WHQJXrezQ3YGIiLqHakDhiB1wBAc\ne/4f9irBLTUVWF+0AgtKVnsdxSv9AXvLrhI8MDFRTMjOVo/NyozJMxiQnZgIRVEkvYre0dTR4Wvq\n6KiUnYPkYdEleF0sukRE0SAtNw9puXnAhdfpAOza9zgYDKLFXIY1RSsxf/PqTmfhskDA3qrEq1VC\nr1YpuUldJTg/s6sEZyYkhEUJrrPbnQCaZOcgeVh0CX5vp7XT6UCMPlF2FCIikkClUiFjqBEZQ43A\nxTfstgpPMBiEpXI7lhevwBeb1nQ61/0SCDraFH13CR6SnCwmZmdpRmdm6vIMBqTr9SFTgs1tbS4A\njbJzkDwsuoRgwGeyNdUga3i+7ChERBRiVCoVsvLykZWXD1z6571KcFP5ZvxUtAKfbi7wdBT8FISz\nXdGrukrw0NQUMTE7WzM6I0OXZzDAEB/fryW42mZzgyO6UY1Fl+CwNpbYGmpE1vD80PgTnIiIwoJK\npULOqGOQM+oYALfG9rwtGAyioXQjvi9agZmbCzqdq9YH4HLsGgnOM6SKCdnZ2jEZGdrhBgNS4+P3\n8yxHzun1OoQQotcfmMIGiy6h0+nY0VJb4QCQJDsLERFFBpVKhYGjJ2Hg6EkA/r7bSLDf70fD9mJ8\nt2ElPt5U4HGuKAgqbqeSoFaJBI1aGWFIFROyc7SjM9K1w9PSkBwbu59nOTBHZyf3uI9yLLoEAGZr\nVWkHWHSJiKgfaDQa5I47AbnjTgCA3Vqs3+9H/Zb1+HrDSry/ucDjWrpKqDrd0KtVIlGrUUYaDGJi\nTrZuVHq6Ji8tDYkx+97Y09HZCY/fz2kLUY5FlwDA0rUNMBERkVwajQaDJ5yIwRNOBPYswV4vares\nxdzilaJ1yzq3c8kKqL0e6NUqkaTVKMa0NDExJ1s3Mj1d4/b50On3l8h5FRQqWHQJQojgiClnOWXn\nICIiOhCNToehk07B0EmnKADiet7m9XhQu7kAnxavEm1b17ptFZtVZputVVJUChEsugQA8LocnMdE\nRERhSxcbi+HHn47hx5+uAIhb9OqDTeUfPPed7FwkV+Tu+0eHpdPlbOh0susSEVFkaN6xzQXALDsH\nycWiSwAAT0f7ssbyzbJjEBER9QpPh71dCBGQnYPkYtElAIC9uW5lzeZ17bJzEBERHS0hBDwOm012\nDpKPRZd22lKzqYBvCkREFPYc1gYEfL5y2TlIPhZdAgAIIdwdrc0sukREFPYsVaVw2qyrZecg+Vh0\naRdXe2tLMBiUHYOIiOio1G/fYHO2WTbIzkHysejSLn5v54bW2grZMYiIiI5K9cZVbQB4hjWx6NL/\ntNVX/1y/fYNXdg4iIqKj4WhpahFC8PcZsejS/wR8ncXVG1a1yM5BRER0pDyOdnR2tFfLzkGhgUWX\nemqwVJfZZYcgIiI6UrXbioSzzbpIdg4KDSy6tIsQQrgdbdwXnIiIwtaO9b9YOlqbl8rOQaGBRZd2\n43V1mN32NtkxiIiIjkjd1qJ2AFxDlwCw6NIe3Pa2XxrKNsmOQUREdEScNqtFCCFk56DQwKJLu3FY\nG1ebSwo4fYGIiMJOe3MdvG7nVtk5KHSw6NKetlRvWMm5C0REFHZqNhV4bY01C2TnoNDBoku7EUL4\n7ZaGJn7qQ0RE4aa84Gerz+3k1r+0C4su7aXT5VjRXLlNdgwiIqLDYjWXtwkhmmXnoNDBokt7aTGX\nf1G26geb7BxERESHKhgIwGWzsuTSblh0aV+Ky9f+bJUdgoiI6FBZqk3welxrZOeg0MKiS3sRQvgd\n1oamYDAoOwoREdEhqS5e1dFiLl8oOweFFhZd2qdOp2NZ8w6u0EJEROFh+4qFzQDWys5BoYVFl/ap\npabii9KVi7nMGBERhbxgIABbo7lBCNEpOwuFFhZd2p+NFeuXtsgOQUREdDB124vhcbR/JzsHhR4W\nXdonIUTAYW1o5DxdIiIKdZt/nNvUVl/1mewcFHpYdGm/Op0dvzSVb5Idg4iI6IBqNq21CiEqZOeg\n0MOiS/vVWlsxt3Tl962ycxAREe2P294Gp63FJDsHhSYWXTqQkh2FnKdLREShq3T1D572ppqPZeeg\n0MSiS/slhAg6WpoagoGA7ChERET7tOWneY2dTsePsnNQaGLRpQNyt7d+bS5ZzTPSiIgo5Agh0FpX\n2SiEcMrOQqGJRZcOyNZYM6Pw64/rZecgIiLaU2PZJng62n+SnYNCF4suHZAQoqmhtKRBCCE7ChER\n0W62LPmqpaWmYpbsHBS6WHTpoFztrQvqthXJjkFERLSbysLlFgDbZOeg0MWiSwfVVl/1fuHXHzXI\nzkFERLRTp9MBp81aLviRIx0Aiy4dlBCitnZLIacvEBFRyCgr+NHrsNRz2gIdEIsuHRK3vfXHpvLN\nsmMQEREBAIq++aTO1d76tewcFNpYdOmQtNRUvFP4zSfNsnMQERF1ujrQVl9VxmXF6GBYdOmQCCEq\nzCWrOU+XiIik2/zTPJetsfZ12Tko9LHo0iFztbctt5rLZMcgIqIot+G7WXUeh22h7BwU+lh06ZBZ\nq0vfLF4w0yo7BxERRa+ONgvsloYSIYRPdhYKfSy6dMiEEFsr1v3SKDsHERFFrw3fzWpvqSl/QXYO\nCg8sunRYXPbWtW311bJjEBFRlNry8/z6gM+7WnYOCg8sunRYLJXbX9+wcHab7BxERBR9Wusq4Wyz\nruQmEXSoWHTpcBVvW/ZtHd9jiIiov62b977FUmV6UXYOCh8sunRYhBDC2Wb9un57MZsuERH1GyEE\nKtb+XC+E2Co7C4UPFl06bK21FS+vnP1arewcREQUPRrLSuC0tSySnYPCC4suHTYhRHNDackOn8ct\nOwoREUWJgi/erW+t3fGa7BwUXlh06YjYm+te3vTjXJfsHEREFPmCwSBqNq+rE0Lw00Q6LCy6dERc\ntpZvCr/+iG84RETU5yrWLvG5bNZZsnNQ+GHRpSMihPA7rI3LW2t3yI5CREQRbsWMF822xpp3ZOeg\n8MOiS0fMUmV6cvmMlxpk5yAiosjVWleJtvrqFUIITpejw8aiS0dMCFFZvWHVDr+3U3YUIiKKUEs/\n/E+9pcr0iOwcFJ5YdOmo2C31L5Z8P8cpOwcREUUer9sJc8maUiEE956nI8KiS0fF2Wadv27eezWy\ncxARUeRZP/+D9rb6qidk56DwxaJLR0UIEXC0NC1sLNskOwoREUUQIQSKv5tl7nQ6fpadhcIXiy4d\ntRZz+bPLPnmRo7pERNRrygt+9nW0NH8ohOCW83TEWHTpqAkhmhq2b9jqtrfJjkJERBGia0kx81uy\nc1B4Y9GlXtFcuf3epR/+xyI7BxERhT8uKUa9hUWXekXA79toWrnI5HVzAQYiIjo6XFKMeguLLvWa\n1rrKB1Z9+jrnLxAR0RHjkmLUm1h0qdd0Oh3LS76fU+73eWVHISKiMMUlxag3sehSr2pvrHl8/fwP\n7LJzEBFR+AkGgyj+diaXFKNew6JLvcrV3vrtuvnvVwYDAdlRiIgozJQs/txptza+wCXFqLdoZAeg\nyCKEEEmZA/5Tsvjztyaee6Vedh76n4oVC7B14QzYaiugjU9A5siJOO6K25A8YBgAYM7fzkFHb8P2\nIgAAIABJREFUS8MBH+OUm5/AyNMvBACY1y/B2hnPwW2zIiNvPE684UGkDBy+2/099jZ8cfuvcdwV\nd2DMr67smxdGRBEhGAhgxcyXK+3NdR/LzkKRg0WXep3D0jBr5ezX7p/wmyvyFUWRHYcAFH72Ckrm\nv4vknKEYfc4VcLU0obLgezRsKcBvn/4ciZmDkH/uNfC5O/Y61t/pweYFH0Kti0HGiHEAAFdrM355\n+S4kZg/B6LMvR+XqRfjh6Ztx0QvfQKOL2XVsyVf/RYw+GcazLu2310pE4an4u5kdDkv9M0KIoOws\nFDlYdKnXCSGCKTmD39i+7Nvnx5x+fszBj6C+ZCnfhJL57yI7/wScc99bUGt1AIAhBdOx5KU7sWHu\nWzj1licw9txr9nn8mg+fghBBnPjH+5EyMA8AULFyAYKBAH714LuIT0nH0CnnYMHDV6N2wzIMnTwd\nQFcZNv34OaZcdx9UGm3/vFgiCksBvx+rZr9eabc0zJadhSIL5+hSn2hvrHln6Uf/qeI0K/m2fT8b\niqLC1D89uqvkAsDQKdNhPOtSJOcM2e+xTduLsG3xbAw85mSMPOOiXdd3NNchNikV8SnpAIDUIcau\n6y31u+6zcd7biE/N3O04IqJ9KfrmY4fD0vAER3Opt7HoUp8QQvjszfUzKouW+2RniXZ1G1YgdfBI\nJGUP3uu2k296BMf87k/7PXbtJ89BpdZgyh8f2O16nT4J/k73rq93TnnQxSUAABzNtShd8iUmXnor\nFBXfZoho/wI+H9bMeXuHo6VxjuwsFHn4G4j6TFt91fM/v/tUlewc0czd3gKPow0pg/Jgq9uBn56/\nHTNvOAkzbzgJS166C47muv0eW7X2B1h3bMbIM36316hv5sgJ8Hlc2LpoJnweFzZ/8wEUKMgcNREA\nsOGLN5GUMxR5p5zXp6+PiMLfuvnv2+3NdY9ypQXqCyy61GeEEG5bQ/U8c8karjUmiavN0nXZ0oQF\n/7oKzpYGjJx2MTKNE1FV8D2+fegqdFj3vdLClm8/hqJSY9z51+91W+5xZyD32NNR8NEzmHH9FGz+\n9iOMu+D6XYW6YuW3mHTZrbvuL4L8NJKI9ub3ebH2y/9WdLQ2fyU7C0UmnoxGfaqlpuKx716676K/\nvLdkJFdg6H/+ThcAoHF7IUacdiFOufnf2Pn/YeuiWSj46GkUfPQMzrrr5d2Oa6nchubSDRg6efo+\npzwAwNn/fA3mwl/gaDQjfcQ4ZBmPBQAUz3kdhiGjMHTydNgbzVj2+v2wlm9CfFoWJv/hnxh64jl9\n+IqJKJwUzHm7vb2p9l8czaW+whFd6lNCCFd7Y+2LJd/PccnOEo0UpetHXKVSY/K196DnHxtjfnUl\nEjMGom7Dcvi9nbsdV77sawDAqIMsCzb4uDMw9rxrd5XclqrtqC74Acde9ncAwNJX70XA68H0+97E\n4OOmYemr96C9gdvXExHg6/Sg8OuPyl22loWys1DkYtGlPmdrNL+99MPnyvYsU9T3dPGJAICEjIGI\n0SftdpuiKEgdMgoBvw/OPaYv1BT9gpiEFAwYd+JhPV/RZ68gY9QEDJp0KlqrTbDu2Ixjf38bBk6Y\nisnX3QtdQhJKf+L5JkQErPn8TVt7U+0DHM2lvsSiS31OCBFsb6y9e/nHL7TKzhJtEjIHQVFUCPr3\nvfhF0O8HAGhiYndd115fCUdzLQYfP+2wVkxoLt2A2g3LcezlXaO59kYzACApu+tENpVKjcTMXDia\na4/otRBR5PB53Cj6dmaZ09byvewsFNlYdKlfOG3WHzcu/nyj09YiO0pU0ehikJ43Dh0tDbuK507B\ngB+t1SbEJqYi3pC16/rmshIA2DUd4VAVfvYKcsZOQc7YybseHwBE8H/nIga8HgCcq00U7ZZ/8kJr\ne2PNvbJzUORj0aV+01Sx5ebFrz5Yf/B7Um/auf1uwUdP7yqfALB5wUdwtTUj79QLdpu721q1DQCQ\nNmzMIT9H/aY1aNy6Dsf+/u+7rkseMAxA10gvAHR2tKO9ofqAG1QQUeSzWxqwcfGcIld76xLZWSjy\ncdUF6jdCiNLMYcZfLNWlV2UMGSU7TtQYecZFMBf+AvP6n/HVvZdi4ISpaK/fgdoNK5CcMxSTLr11\nt/s7mmoAAPGpmYf8HEWfvYJBE09F5sgJu65LGzoa6cPysfaT59BqLkXT9kJACBjPvrx3XhgRhaVv\nnruzpnnH1htk56DowBFd6leWqtLbv/3P3ZWyc0Sbaf94AZOvuQcAsP37T9FaXYox51yJ8x6fAW2c\nfrf7dna0Q4ECXXzCIT22ufAXWCu6Tjrb05l3vYzMURNQ+tMcBLydmPaPF5CQMeDoXxARhaXKohW+\nBtPGz4UQNbKzUHRQeLIj9be03LwnL374rbvzjj9DJzsLERH1j2AggDf/eJqpdsu6CUIILsND/YIj\nutTvWmt3PL741Yd2BLlbFhFR1Fg1+zWbvbnufpZc6k8sutTvhBCd9ubaZ4oWfNIhOwsREfW9jjYL\n1n/1YUl7c9082VkourDokhTtTXUfr5z1qsnnccuOQkREfey7F+6pa6rYcr3sHBR9WHRJCiGEaKur\n+uvi1x5qkp2FiIj6Tu2W9YGazWu/EkLskJ2Fog+LLknj6WgvMK1YuLB++4bAwe9NREThJhgM4pv/\n3FVhrS77p+wsFJ1YdEkqq7nsr/Oe+mt5wO8/+J2JiCisrJv3nr29seZxIYRLdhaKTiy6JJUQwmVr\nMN++5L2nuTcwEVEEcTtsWP3Zm1vam2pnyc5C0YtFl6RzWBsXlyz+fKmlupSLOhMRRYiFL91Xb6ky\nXS+4YD9JxKJLIcFSZbrxy8dv5tq6REQRoKp4hbeycPlXAZ/XJDsLRTcWXQoJQghba13lg6tmv2aT\nnYWIiI5cp6sDXz/7j+1Wc9kdsrMQsehSyGhvqv1s7Zf/LbA1mGVHISKiI/T1s3fUW6tLrxZCeGVn\nIWLRpZBiqdx+7Zf/vrmKU7qIiMKPacUid1Xxyg+9Htdm2VmIABZdCjFCiGaruez/ir75xCE7CxER\nHTpXeysWvvLA5tbaiodlZyHaiUWXQk5bffVbyz95YWNHa7PsKEREdIi+/PctNU3lm68QQnATIAoZ\nLLoUcoQQoqli69XznvwrJ+sSEYWBDQtndzSYNrzCbX4p1LDoUkgSQpgbS0veKP5uFqcwEBGFMLul\nAUvee7aota7yedlZiPbEokshq7Wu8v+W/Pfp1VZzOc9MIyIKQUIIfPHoTVXNO7ZewY0hKBSx6FLI\nEkIIS5Xp8s/+dV2Z39spOw4REe1hzZy32i1VpseFEA2ysxDtC4suhTQhRHuLufzar5+5nW+iREQh\npLWuEqs/fWN1W331B7KzEO0Piy6FPFd7a0HF+qXvlnw/xyk7CxERAcFAAF88clOFpcp0tewsRAfC\nokthobW24tEf3/p3QWtdpewoRERRb/Fr/7K01FTcIYRolZ2F6EBYdCksdM3X3X7pZw9cW+b3cVdJ\nIiJZti9f6N6y5KsZ7c11C2RnIToYFl0KG0KINqu5/IYFz/2jUXYWIqJo1FZfjYUv3buuxVx+t+ws\nRIeCRZfCitNmXVG25uf3N/88n/N1iYj6kd/biU8fuKasuXL7RUKIoOw8RIeCRZfCTmttxUM/vP7w\nelsDN04jIuov85/6W72lynQN5+VSOGHRpbAjhAg2V26/5NMHry0P+Hyy4xARRbzCrz+yVxYue9XV\n3logOwvR4WDRpbAkhGixVpX++dsX/tkkOwsRUSSr21bs/+WD535qrat8VnYWosPFokthq6PNssS0\ncvEb6+a93y47CxFRJOpos2DOIzduslaXXsUtfikcaWQHIDoaLTXlj2cMGWVMy827ZPjxp8fIzkNE\nFCkCPh9m3XPljqbyzRcIITyy8xAdCY7oUtizmsuum//U39ZYzWUcbSAi6iVfPXNbQ/OObTcIIepk\nZyE6Uiy6FPaEEH5LlemCWfdetcVpa5Edh4go7BV88W57xbolL3a0Ni+VnYXoaLDoUkQQQjgaTBvP\nm3n35Tu4cxoR0ZGrLFruXTHz5YUtNRXPyc5CdLRYdCliCCHMjeWbr57z8A31PGeCiOjwNZVvCcx7\n4taV1urSa2RnIeoNLLoUUVztrWvMJWvu++GNRy2ysxARhRNbYw1m33/1BkuV6XwhhF92HqLewKJL\nEaetvvqTjYs/e3/DwtkO2VmIiMKBq70VM+66bGtTxdZfCSFcsvMQ9RYWXYpIrbU77v/p7ScWV29c\nw63TiIgOwOdx45M7L6mo21b0GyEEz+iliMKiSxFJCCGs5rKrvnz8z+tb6yplxyEiCkkBvx8z77nC\n3GAquUQIYZadh6i3sehSxBJC+Jort583859XbHPb22THISIKKUIIzH38z/X124tv9DjtG2XnIeoL\nLLoU0YQQbfXbi8/98LYLyz0ddtlxiIhCxqJXHmyuLFx2n93S8KPsLER9hUWXIp4Qoqpua+FvPrr9\nwopOV4fsOERE0q2a/Xrb5p/mvdxWX/2J7CxEfYlFl6KC3+ctr9tWfMHHd1xU6XXzhGIiil6bfvzS\nueqz12e11JQ/JTsLUV9j0aWo4XU7t9VtK/rdJ3deUu3r9MiOQ0TU7yqLlnu/f/3hRS3m8r/LzkLU\nH1h0Kap4OuwltVvXXzLjrkvN3CqYiKKJuaTAN++JW5dZq0uvENw+kqIEiy5FHbfdVlized3vZ959\neU3Ax2V2iSjyVRWv8M555MafLVWmc7nrGUUTFl2KSq721jU1m9ddM/u+q2sDfr7nE1Hkqli7xDP3\n8b/8YK0uPV8Iwb/uKaqw6FLU6mhtXmretOb6z/51XV0wEJAdh4io15Wu/sE9/6m/LbJWl13IkVyK\nRiy6FNXsloYfq4tX3jznkRsbgsGg7DhERL1m29IFrm+evWOB1Vx2iRCCf81TVGLRpajX3ly3oLJw\n2d/m/fvmRp6fQUSRYPNP85zfvnjPfKu57AohBP+Kp6jFoksEwNZY82X52iV3znn4hgbO2SWicLZh\n0acdi159cE6LufwPLLkU7RSOYBH9T1JGzvSBY45976pnZ+dqY+NkxyEiOiyF33xsX/LeM5+2mMtv\n5hJiRCy6RHuJ0SdOHGCcOPfaF+cOj0tKlR2HiOiQrP3yvfZlHz3/SUtN+W0suURdWHSJ9kFRlCED\nxkxaeO0Lc8ckZw2SHYeI6IBWffqGbdXs196zmsvulp2FKJSw6BLth6Io6Vl5Yxdf9eysiZnDx3A+\nOxGFHCEEfn7nSWvxd7PesZrLHpSdhyjUsOgSHYCiKPqMocbvLn74rROHTpyqk52HiGingM+HLx77\nU0N18cqHWuur3pOdhygUsegSHYSiKNr0wSPn/Pr2p84eO+1Cvew8RERuhw0z7rq8qrGs5FqnrWW5\n7DxEoYpFl+gQKIqipA0e8dYpV9/2+xMvuzlZdh4iil6ttTsw854rttZv33C+EKJSdh6iUKaRHYAo\nHHSfwfyX9MEjGuyWhr9Ov+XRdEVRZMcioihTvWGVd+6/b15nqdx+vhDCJjsPUajjiC7RYTIMGnZT\n3gnT/n3Rg29kq9Rq2XGIKEoUfzvL8fO7Ty62msuuEkL4ZOchCgcsukRHIClzwK+zhue/ddWzs4Zw\nrV0i6ktCCPzwxqOWku8/f6+lpuIBrpFLdOhYdImOkKIow3JGTfjq909+ODYrbyyXHyOiXuf3eTHn\n4RvqzSVr7murr/5Edh6icMOiS3QUupcfmzv91sdOHX/2xfGy8xBR5HC1t2LG3ZdXNpZtusrV3rpG\ndh6icMSiS3SUuldkeHbcmb/74zl/eyJDpeLgLhEdHUuVScy+/w9bGkwbzxNCmGXnIQpXLLpEvSQ1\nZ/AlWSPGvXjFkx/nxiZyBTIiOjLrv/rQvvTD/6y0Vpf+XgjhkJ2HKJyx6BL1Io1WNyJzeP6Xlz3+\nXn7OqGO4JAMRHTKv24V5T9xSZy5Z81JrXeXzPOmM6Oix6BL1MkVR4tOHjJp12rV3TjvhohuSZOch\notDXvGNb8POHrjdZqkuv6nQ6NsjOQxQpWHSJ+oCiKIph0PA7h0w8+c6LHnxjgDYmVnYkIgpRa798\nr33FjBeXWapKrxJCdMjOQxRJWHSJ+lBcYsoJGcOMM6546pNRhoHDZMchohDidTsx97G/1NVsWfuf\ntrqqlzlVgaj3segS9TFFUVIzho3+7JSrbptywsU3JnHrYCJqqtgS/PyhG7ZbKrf/3utxbZadhyhS\nsegS9QNFUZTUnCF/yswb88Clj/53SIIhU3YkIpJACIG1c99tXzHzlV+s1aVXCyGcsjMRRTIWXaJ+\npCjKoKy8/M/O+svDE7nBBFF06XR1YO5jf66t3Vr4f621O16VnYcoGrDoEvUzRVFUhkHD7x+Uf9zN\nFz34xiCuuUsU+aqKV/q+/r87TNaq0iu8HtcW2XmIogWLbhgyGo0fArgWwESTyVRymMeqAdwC4H2T\nyeTqg3i9zmg0ngAgxWQy/dD99VAAOwB8ZTKZLpKZ7WiotbpRmcNGf3reXc/lj5h8ZozsPETU+zpd\nHVjwn7saKouWf9FiLr9bCOGVnYkomnCv0vA0D8CjAJqO4NhZAF4BoOnNQH3FaDSeB2ANgDE9rm4D\n8BiA2VJC9ZKAz1vaWLZp8vwn//rGvCdubfR53LIjEVEvMq1c3PnW9adv2LDw019Zq8tuY8kl6n9h\nUXZodyaT6SsAXx3h4VkAwmkYPwPAbssUmEymdnQV3bAnhPADuDMuMWV27dbCDy568PXRg8Yezx3V\niMKY296G+U//vbZ2y/qPWmt3PCKECMjORBStWHSjV7itcRVueQ+L22FbpyjK8Z89eN2bY8+88Lzp\ntzyWodZqZcciosNU8sMXrp/feXJTU8WWPwghymXnIYp2nKMbhnrO0QVgR9d81ccAFAP4F4BxABzo\nGvW932QytXQfF9zjoZaaTKZp3bfpANwF4BoAw7of9wcAD5lMpsoez/0ogIcBnA3g6e4MlQCONZlM\nLqPROLX7cU4CYADgBLAewFMmk+mXPV5HOoAHAfwOXSPNNQDmAHjaZDI5e7zOnoaia8rNXnN0jUZj\nDoBHAJwHIBNdUzu+BfCYyWRq3MdrGAPgOgB/6L5/OYBXTSbT25AowZA5PXXAkFfOu/O5vKGTprLt\nEoUBh7UR8568tbqxfMvrbXWVzwsh9ny/JSIJOEc3clwA4EsAdQBe7r68CbtPcXgMQHX3fz8D4AMA\nMBqNWgALATwJoB3AqwAWAbgEwDqj0Th2H883E10l9mUAS7pL7oUAlgKYDGAugBcArAJwFoDvjUbj\nhJ0HG43GbHQV4NsBVAB4DV1F9wEA87tPmpvXI/8idM1Lbu+RQfR4vDx0Ff0/A9iKrnnI2wH8BUCh\n0Wjc17ZkM7q/RwsAvAtgIIA3jUbjTfu4b7/paG3+oXbL+olfPPqnV2bec2WNw9p48IOISAohBNbN\ne9/+35t/vXTb0gWntNbueI4llyh0cOpC5DgWwGUmk2kuABiNxofQVfxONhqNRlOXx4xG4zQAgwE8\nYzKZ7N3H3gFgGoBnTSbT/Tsf0Gg0voKuovo+gCl7PJ/ZZDKducd1z6LrRLFJJpPJ0uNx/tl92+UA\nNnZf/X/dOe4wmUyv9LjvW+gqqxeYTKb5RqMxFcCFABbtvF/3dXt6B12jsjeZTKb3ezzezQDeQFeR\nPXuPYwwAxvQY8Z4FYCWAGwH8dx/P0W+EEJ0A7lYU5ZXmHVvfm3juVceeds2dBk5nIAodtgYz5j7+\nl0qruexpW4P5v9zClyj0cEQ3clTsLLkAYDKZ/AB+6v5yyEGOvRFdBfXBnleaTKZCdE0lOMFoNI7Z\n45h5Pb8wGo0KgPsAXNOz5HZb2n2Z0X3fGAAXAyjtWXK7PYWukeWGg2Tu+dy56Crqy3qW3O7X8BaA\ndQDONBqNg/c49P2dJbf7vqvRNWJ8sO9XvxFCmJsqtk5fNevVP7z5x1O3VqxbwrO2iSTzedz4/vVH\nLB/8/YJvywt+mtJWX/0uSy5RaOKIbuQo3cd1Oz/m3+8arUajMQHAKACNAB42Go173iW7+3IigG09\nrq/seSeTySQAzO9+zCHomiecByAfXSUUAHauJpAHIB7A6j2fzGQymQE8tL+8+zGx+3LZfm5fBeAE\nABMAmHtcv6/vmR1AwmE+f59ztDQtVBTlpy//fcsj2SPGXf3be18akpw1SHYsoqgSDAZRtOCTjpUz\nXzG11Vfd5umwr5KdiYgOjEU3cnTu47qdIwwHWrFg57Zc2eg6QWtfBLo+5u9pr0VfjUbjeHTNjT29\n+yofgC3oGlEd2SPHzqkHdvSOpO7L9v3cXt99ueeWu/v7noXkCg/da3A+qCjK61Zz2bvjp1865Ywb\n7k3TaHWyoxFFvMrC5b5Frzywo7257qn2ptpPOIJLFB5YdKmj+3KZyWQ640gfxGg0JqJrlYZEdK26\n8AOA7SaTyW80GicDuGofz5m4n8fSm0wm52E8vaP7cuB+bt9ZrFv2c3tYEULUAzhPn5px5ralC146\n56+PjzJO/RV3ViPqA1ZzOb59/u6q5srts1trdzzWPX+eiMIEi2702W0UwmQytRuNRjOAcUajMdZk\nMnl63m40Gq9E12jsRyaTqRr7dya6TgZ7zmQyvbjHbfndlztHSk0AvOhanWE3RqNxIIAao9H4jslk\nunnPvPuxofvylP3cfhqAILpWY4gYzjbLz4qiHPvVM7fdnzls9PXn/uPZYZnDRsuORRQR3PY2LHr1\nX/WVhcuWWapMfxdCWGVnIqLDx5PRoo8PXYWz5wjgh+iamvBM90llAACj0ZgP4C0Ad+Lgo6E7pzJk\n97yy+wSwR7q/1AJAd5meCyB/H0t5PdB9+WOPvMAB5hl3z+tdAuD47lUWej7/TQBORtcSaPX7Oj6c\nCSH8rbU7/m1asfCET+689P2Z/7yi2mrmGvVER8rv82LpB8+1vn3jmcvWf/XhWc2V269kySUKXxzR\njT613ZfvG43G700m06voWlP3VwBuA3Cq0WhcCiAFwGUA4gBcbTKZOvb5aP+zAkAVgGu6N4IoAZAL\n4LfoWrJrEID0Hve/G10jsO8YjcaL0TWXdzKAUwHMM5lMX+yR9xaj0ZgG4KX9PP9fACwH8Eb3420C\nMB5dS4rVoWvJsoglhGgBcKOiKJlNFVueyzFOmPbrvz+RmzpgqOxoRGFBCIEtP81z/fze0xXtTXX3\nONssi2RnIqKjxxHd8CRwaB/p7+t+TwIoADAdwK3ArhHWaegaeY0FcAuA36CrOE4zmUyfHey5TSaT\nq/sxvwRwHLo2gjCia77uuegqsqcYjcb47vs3oKvYvg3gmO775wL4N4ArejzuMgCvo2vE+VZ07Wa2\nF5PJVA7geHStl5sP4K/oWt3hZXSt69tzlYgDff/C+gQTIURzc+X260oWf37iB3//7aefP3R9ra2x\nRnYsopAlhIBp5aLOt2+cVvrti/fc32DaOIkllyhycAtgogimKEpu5rAxL+aOn3ziOX99fGBSRo7s\nSEQhIRgMYvOPc10rZrxc7bA2vGNrrHmTJ5oRRR4WXaIooCjK0Mzh+S8NnTT1hOm3PDIgwZApOxKR\nFAGfD0ULPnGs+eLt6g5r04t2S/3HQgi/7FxE1DdYdImiiKIoI7Ly8l8efvwZx571l4ey9SlpsiMR\n9Qufx42CL96xFX87c4fd2vBUR0vTl1wLlyjysegSRSFFUUZnjRj7wsDRk8afccO9gzKG7rUjHlFE\n6HQ6sHzGSy1bfp5fZmusedjjsP3IgksUPVh0iaKYoigD0gaPeCglO/esU66+ffCoqb+OUal4jiqF\nP6etBb+8/2xz+Zoft7XU7rjP63aukZ2JiPofiy4RQVGUmKTMgTckpmX9+ZhzLhsy+ZI/pcYmJB38\nQKIQYzWXY9lHz9eZN60psVSa/hnw+7bIzkRE8rDoEtEuiqIo2pi4qWm5eY8MzD9uzBk33DMwffBI\n2bGIDijg92PbsgWe1Z+9UWNvrvvFWl32hBDCLDsXEcnHoktE+6QoyqD0ISMfTsnKnXbKNXcMHnnS\nOTpOa6BQYrc0YNXs1yxla36oslsa3+hoaZwlhPDKzkVEoYNFl4gOSFGUuOSsQTclGDJvmvCbKwYf\nf+EfU+ISU2THoigVDARQtvoH3+rP36xpq6sqaq7c9qgQgtMTiGifWHSJ6JAoiqLo4vSnpw4YemdK\ndm7+cb+9buDo086L1cbEyo5GUaC1rhKrZr/WWFm0wtzR0vS+3VL/iRDCJTsXEYU2Fl0iOmyKosTq\nU9MvSkzP+XPGsNHDplx806Bhx5+u5tQG6k0+jxubfpzrLPz6o1q7pWGFtbr0GSFEuexcRBQ+WHSJ\n6KgoipKSnJ17rT4l/apB+cflTrn0TwNyjBOhKIrsaBSGOp0ObF36jbtk8ef1tsaaCrul/i2XreUb\n7l5GREeCRZeIeo2iKAMNg/Ju1aeknZd3whkDTrj4xgzDwGGyY1GIc9vbsOXn+c6SH+bWOyz1pbbG\nmnc8He2LhRCdsrMRUXhj0SWiPqEoypiMocZ/6FPSp+ZP++2AsWf+LoWll3bqaLNg0w9fOLYu+are\nYW3a1lpX+YbP41rCkVsi6k0sukTUp5SuOQwnZAw1/jE2MXlK5rDRmeOnX5Y9/LjTNNrYONnxqB/Z\nLfXYuPhz27alCxqdbdaNLTXlbwR83pVCiIDsbEQUmVh0iahfKYqSlpCWdWGCIfPyBEPG0BFTzs4c\ne+aFqWm5IzivN8L4fV7UbForti//ttm8qcDqamspslSXviaCgXWCv3yIqB+w6BKRNIqiqABMTB8y\n6prYhKSp6YNHZh9zzmXZeZOnaXVxetnx6DAFfD7Ubi3E9uXfNtdsWtvqtFlrXO2t39qb674BUMly\nS0T9jUWXiEKGoigp8cmG3ySm51ytN2TkDTv21LS8E6ZlDBxzLHRx8bLj0R4Cfj/qTRt9F+QIAAAD\niUlEQVRgWv6dpWrDqlZnm6XO7bAtbG+s+RpAGYstEcnGoktEIal7bu/oxIycs+KTDNNjE1MGJWUO\nSMs7/oy0IRNPSsgcNgYqtVp2zKgS8PvRWLYJphULWyqLlu8stj/YGszzAWxjsSWiUMOiS0RhQ1GU\neACT0nLzfqON05+sT07LTBs8ImXkiWdl5I6fokvOHCg7YsTwONrRULYJNZsK2syb19qdrZY2t8Nm\n9bqdK9rqKr8EsEUIEZSdk4joQFh0iSisKYqSpY3Tn5iSNeh8XXxCvj4lLWNg/nHJuWNPSMsYNlqd\nOmAo1BqN7JghSwiBtvoqNJg2+quKV1qbdmxzuB1tLV5XR43b3vaLw9q4GsBWrmlLROGIRZeIIkr3\nCW4jtLHx45OzBp2s1urGxSYkGWL0SSlpg4bHD8w/NjVj6Og4w6Bh0KekR81KD54OO2wN1bA11gSt\n1WX2um1FDoe1weayt7X4Oz2bbY01P/k8rkIAtZyCQESRgkWXiKJC95zfLABjkrNzJ8bEJ05Ua7S5\nMfrEpJj4hKTEjJy4rLyxCWm5eUn6lHRVfEoa9KnpiE1IDvkyHAwE4LA2wNZYg9a6yk5rdWm7tbrc\n43bYXD6Ps6PT1dER8HmbAz7vNrulYZPP46pE18liDtnZiYj6EosuERG6VnwAMEyt1Q1JSMsarIuN\nH6yo1AMUlZKmjYmL1ehi4jS62Fi1Vheri42PSUjP1iRl5MQkZQzQJxgyY+NT0qDWaKHSaKBSqaGo\n1FBrNFBUaqjUGqjU3depNVDUaqi6r1fUagQDfnQ6HfC6HPB02NHp6kCn0w6Po93rdtg8bnur1223\n+VztrQFPR3sg4PP5gwGf1+/z+nweV4fX5XQEg4HKTqdji725bisAM7pGZj2Sv61ERFKx6BIRHSZF\nUbQADAAyAKTHxCdmx6cYchWVJk6lUmmhQKsoigZQtIqi0gBCi66vNYDQANjzv32AsIlg0Bbw+6w+\nj9vqtrc1+72edgAOAPYel07uJEZEdGhYdImIiIgoIqlkByAiIiIi6gssukREREQUkVh0iYiIiCgi\nsegSERERUURi0SUiIiKiiMSiS0REREQRiUWXiIiIiCISiy4RERERRSQWXSIiIiKKSCy6RERERBSR\nWHSJiIiIKCKx6BIRERFRRGLRJSIiIqKIxKJLRERERBGJRZeIiIiIIhKLLhERERFFpP8HP6nrdsRE\njLMAAAAASUVORK5CYII=\n",
"text": [
""
]
}
],
"prompt_number": 48
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice how the residual error variance is 27%. This is $100 - R^2$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Aside: Types of Sums of Squares\n",
"\n",
"In unbalanced designs, the \"Types of Sums of Squares\" become relevant. The different Types of Sums of Squares (Type I, II, III...) basically correspond to different ways to compare models up and down the graph diagram shown above. That is, they correspond to different nested model comparisons -- and with unbalanced data, it matters which comparisons you do.\n",
"\n",
"A complete discussion of the different sums of squares is beyond what we want to cover in this course. For some starter reading, see [here](http://goanna.cs.rmit.edu.au/~fscholer/anova.php), [here](http://www.uni-kiel.de/psychologie/dwoll/r/ssTypes.php). \n",
"\n",
"The main thing we want to emphasise is: **ANOVA is about model comparison**. Think about what hypotheses make sense to test for your dataset."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ANOVA as model comparison: summary\n",
"\n",
"We've provided here a brief overview of one way to think about ANOVA: it's comparing the incremental change in the sum of squared residuals as you add terms to the model. In this sense, ANOVA is about **model comparison**, using the change in residuals to decide whether a term is important. \n",
"\n",
"F-ratios are not the only way we could compare the models. ANOVAs are just a historically popular way to do model comparison -- mainly because they are computable by hand. One advantage of thinking about ANOVAs as model comparisons is that once you understand the idea, it's easier to think how to use other types of comparisons.\n",
"\n",
"For example, if we wanted to emphasise parsimonious models (i.e. keeping the number of parameters low), we could consider using an information criterion like the [AIC](https://en.wikipedia.org/wiki/Akaike_information_criterion) or [BIC](https://en.wikipedia.org/wiki/Bayesian_information_criterion) to compare the models. The structure of the comparisons and their hypotheses are the same, however. For a comprehensive treatment of this type of model comparison, see [3].\n",
"\n",
"Another possibility is using *regularisation* -- we fit the model but penalise the fit for having many coefficients far from zero. For example, [*Lasso regression*](http://statweb.stanford.edu/~tibs/lasso.html) (also called L1 norm regression) penalises the sum of the coefficients from growing too large. This has the effect of setting many coefficients to zero -- i.e. removing those terms from the model. This is equivalent to using a model lower down the nesting hierarchy drawn above.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ANOVA in GLMs\n",
"\n",
"We've demonstrated ANOVA for the linear model. The concepts inherent in ANOVA can be extended to Generalised Linear Models with non-Gaussian errors and non-identity link functions. To do this, we need to consider a generalisation of the squared residuals, called the *deviance*. \n",
"\n",
"The *deviance* is the difference in log likelihood (i.e. the likelihood of the data under the model) between a *saturated* model (that contains a parameter for every data point, and thus predicts the data perfectly) and the model in question. Intuitively, think of *deviance* as the \"amount by which our model predictions *deviate* from a perfect prediction\". In the Gaussian + identity model, the deviance reduces to the sum of squared residuals:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# anova table from interaction model:\n",
"print(anova_table)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
" df sum_sq mean_sq F PR(>F)\n",
"C(design, Sum) 1 0.385641 0.385641 0.484016 0.499861\n",
"C(paper, Sum) 1 1.718721 1.718721 2.157158 0.167628\n",
"C(design, Sum):C(paper, Sum) 1 23.386896 23.386896 29.352777 0.000156\n",
"Residual 12 9.561029 0.796752 NaN NaN\n"
]
}
],
"prompt_number": 49
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# summary table from GLM fit:\n",
"fit_glm = smf.glm('distance ~ C(design, Sum) * C(paper, Sum)', planes).fit()\n",
"print(fit_glm.deviance)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"9.561029\n"
]
}
],
"prompt_number": 50
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For GLMs with non-Gaussian errors and non-identity link functions, ANOVA becomes \"Analysis of Deviance\", but the conceptual principles remain the same. We perform nested model comparison, comparing how the deviance is reduced by adding terms to the model.\n",
"\n",
"Unfortunately Statsmodels has not implemented analysis of deviance as a function yet (as of version `0.6.1`), so it's not as straightforward to do as in the linear model cases outlined above. However, analysis of deviance can be easily performed in [R](https://stat.ethz.ch/R-manual/R-patched/library/stats/html/anova.glm.html) or [MATLAB](http://de.mathworks.com/help/stats/generalizedlinearmodel.deviancetest.html?refresh=true), or by hand in Statsmodels using model comparisons as before. Remember to take degrees of freedom into account."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# References\n",
"\n",
"[1] Venables, W. N. (1998). Exegeses on linear models. In S-Plus User\u2019s Conference, Washington DC. Retrieved from [here](http://www.sefsc.noaa.gov/sedar/download/S10RD01.pdf?id=DOCUMENT)\n",
"\n",
"[2] Howell, D. C. (2002). Statistical Methods for Psychology (5th Edition).\n",
"\n",
"[3] Burnham, K. P., & Anderson, D. R. (2002). Model selection and multimodel inference: a practical information-theoretic approach. New York: Springer."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from IPython.core.display import HTML\n",
"\n",
"\n",
"def css_styling():\n",
" styles = open(\"../custom_style.css\", \"r\").read()\n",
" return HTML(styles)\n",
"css_styling()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
""
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 1,
"text": [
""
]
}
],
"prompt_number": 1
}
],
"metadata": {}
}
]
}