diff --git a/docs/_quarto.yml b/docs/_quarto.yml index dfdd1f2e7..6bb93e348 100644 --- a/docs/_quarto.yml +++ b/docs/_quarto.yml @@ -31,6 +31,10 @@ website: file: changelog.qmd - text: "Contributing" file: contributing.qmd + - text: Tutorials + menu: + - text: "Variance Reduction in AB Tests" + file: panel_variance_reduction.ipynb - text: Learn more menu: - text: "Regression Tables and Summary Statistics" diff --git a/docs/panel_variance_reduction.ipynb b/docs/panel_variance_reduction.ipynb new file mode 100644 index 000000000..604fde94f --- /dev/null +++ b/docs/panel_variance_reduction.ipynb @@ -0,0 +1,1838 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7671c77a", + "metadata": {}, + "source": [ + "## Panel Estimators for AB Tests with Repeated Observations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b206e77", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "
\n", + " \n", + " " + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "
\n", + " \n", + " " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "from IPython.display import display\n", + "from tqdm import tqdm\n", + "\n", + "import pyfixest as pf\n", + "from pyfixest.utils.dgps import get_sharkfin" + ] + }, + { + "cell_type": "markdown", + "id": "0e15e974", + "metadata": {}, + "source": [ + "In this tutorial, we show how to use pyfixest to reduce the variance of your estimators in AB tests with repeated observations.\n", + "\n", + "For example, we may be a streaming platform and we want to estimate the effect of a new feature on the amount of time users spend watching videos. To do so, \n", + "we randomly assign the treatment to half of our users. For 15 days prior to the experiment, we track the desired outcome (minutes watched) for each user. If\n", + "users are not seen on the platform on a given day, the number of minutes watched is 0. Our experiment runs for 15 days. All in all, we have 30 days of data for each \n", + "user. \n", + "\n", + "To get started, we simulate a panel data set of 100_000 users, with mentioned 30 days of data, with 15 days of pre and post data. " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "337fdce4", + "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", + "
userdaytreatminutes_watchedever_treated
00008.1144880
10107.6330580
20206.7123320
30306.2308200
40406.0044890
\n", + "
" + ], + "text/plain": [ + " user day treat minutes_watched ever_treated\n", + "0 0 0 0 8.114488 0\n", + "1 0 1 0 7.633058 0\n", + "2 0 2 0 6.712332 0\n", + "3 0 3 0 6.230820 0\n", + "4 0 4 0 6.004489 0" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def get_data(num_units: int, sigma_unit: int) -> pd.DataFrame:\n", + " \"\"\"\n", + " Random example data set.\n", + " num_units: int\n", + " The number of users\n", + " sigma_unit: int\n", + " The user-level idosyncratic error term.\n", + " \"\"\"\n", + " data = get_sharkfin(\n", + " num_units=100_000,\n", + " num_periods=30,\n", + " num_treated=500,\n", + " treatment_start=15,\n", + " seed=231,\n", + " sigma_unit=18,\n", + " )\n", + " data = data.rename(columns={\"Y\": \"minutes_watched\", \"unit\": \"user\", \"year\": \"day\"})\n", + "\n", + " return data\n", + "\n", + "\n", + "data = get_data(num_units=100_000, sigma_unit=18)\n", + "data.head()" + ] + }, + { + "cell_type": "markdown", + "id": "b125d2a3", + "metadata": {}, + "source": [ + "We can inspect the data generating process via the `panelview()` function: " + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "005e6a4d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAHbCAYAAABFidl8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJzFJREFUeJzt3Qm4ndO9P/AVGaghEaVJSCTmeR5Ke40hSM3ci15EEFPNNTzxlLSXXkNvQ5FWrxalbqul6FXVFkGLXkXDpTqYKhoSUyYhJNn/57f+/33++wxZOSdOss/w+TzPFvvd79577fW+Z7/fvd611tujUqlUEgDAQiyzsAcAAIKwAAAUCQsAQJGwAAAUCQsAQJGwAAAUCQsAQJGwAAAUCQsAQJGwACy2Y445Jg0bNqzexei2dt1117TpppvWuxh0A8ICS1SPHj1adXvooYeWSnm+/e1vp5tuuil1Bn/605/SV7/61fTqq6+2+bnnnXdertfDDjtsiZStu5szZ07eNm3db6dOnZrOOeectOGGG6bll18+rbDCCmmbbbZJl1xySZo+fXrqyKZMmZI/86RJk+pdFOqgVz3elO7jlltuaXT/5ptvTr/5zW+aLd9oo42WWlhYddVV8y/izhAWvva1r+Vfj2359R6Xe/nRj36Un/Pf//3fadasWWmllVZaImW8/vrr04IFC1J3DAuxbUJsn9b4wx/+kEaOHJlmz56djjzyyBwSwpNPPpkuu+yy9Mgjj6Rf//rXqSOHhfjMsV9tueWW9S4OS5mwwBIVX4q1fv/73+ew0HR5S1/G8cuLtotfu6+//np68MEH01577ZV+9rOfpVGjRi2R9+rdu/cSed2uJloNDjrooNSzZ8/0xz/+Mbcs1Pr617+eg1dHNG/evG4ZCGkirjoJS8uXvvSluMppo2W77LJLZZNNNqk8+eSTlZ122qnyqU99qnLGGWfkxz788MPKRRddVFlnnXUqffr0qQwePLhy7rnn5uW1brjhhspuu+1WWW211fJ6G220UeXb3/52o3WGDh2a37v2Fu8dbrzxxnz/t7/9beW0006rrLrqqpV+/fpVTjjhhMrcuXMr7733XuWoo46qrLzyyvkWZViwYEGj158/f37lyiuvrGy88caVZZddtvKZz3wmP//dd99tVo4vfOEL+b222267vO5aa61V+cEPftCwTrU8TW8TJ05cZB0fd9xxuQxhn332qey5554trnf11Vfn9aK+4zNts802lVtvvbXh8ZkzZ+btEOWNOo263WOPPSpPPfVUwzqjRo3Kj9d6++23K0ceeWRlpZVWynV49NFHVyZNmpTLH5+r9rkrrLBC5fXXX68ccMAB+f+j3r/85S9X5s2b17DeK6+8kp/7jW98o3Lttdfmuooyx+d67bXX8nb4t3/7t8oaa6xRWW655Sr7779/5Z133mn2ee+9997KP/3TP1WWX375yoorrlgZOXJk5bnnnmu0TmvKVC1P09u4ceMWuk0uu+yyvE5t/S7KhAkT8vaJuh80aFDllFNOyfthS387zz//fGXXXXfN9bL66qtXLr/88mavN3Xq1Mqxxx6b98vY5zbffPPKTTfd1Gid2rqOfXnttdeuLLPMMvn/W/rM1e3517/+tXLwwQdXBgwYkF87tsVhhx1WmT59eqs/Lx2bsECHCAsDBw7MB6M4UH/3u9+t3HXXXfngO2LEiPzlfuaZZ+blp556aqVXr175i7xWHHSPOeaY/KV2zTXX5OfF+8TBperOO+/MYWPDDTes3HLLLfn261//utHBecstt6zsvffe+Ys6wkEsO++88/JB5otf/GIOIPvuu29eXntwD8cff3wu25gxYyrXXXdd5fzzz88HmyjbRx991LBeHFw32GCD/MV6wQUX5DJuvfXWlR49ejQcvF566aXK6aefnt8n1qmW98033yzWb4SoOPBffPHF+f7NN99c6dmzZ+WNN95otN5//ud/5tc+9NBDc71+61vfyiEj3rMqPm8cqM4+++zK9773vXwA2m+//So//OEPFxoWYpvtuOOO+T1jW8Vni4P6Flts0WJYiIN7HOziIPad73yncsghh+T1aoNe9QAW2yYOnuPHj6985StfyWXbYYcdcv187nOfy+Enyh/1OHr06EafN+ohlse2jf0jPsuwYcNyXcXrt6VMs2fPzstj2UEHHdSwbZ555pmFbpcoXxzII3i2RgSPeP0IZ1HeqMuo06b7UvztRDgYMmRIDnZRxt133z0/N8JR1Zw5c3KA7t27d+Wss87KdRXBPNa76qqrmtV11HMEhQg58Tf16quv5kAWj0UArn7m2E/jM0WAi3JccskleV/52te+lssaz6NrEBboEGEhlsUBtlZ8GcWvmvgFXivWi/UfffTRRl+GTe211175C69WHASqrQm1qmEhnlPbYhAHvjjInHTSSQ3L4hdmhI7a14kytvTL8b777mu2vNrC8cgjjzQsmzZtWv5FFr9gq37605+2ujWh6vbbb8/P+dvf/tbQOhAHv/jCrxVhK+qiJFoFYnuVNA0Ld9xxR7MDUASI6gGsaViIZXEQqrXVVlvlVo6mB7AIk7W/VMeOHZuXRxD5+OOPG5YfccQROUhUW59mzZqVQ0GEuFoRvOIz1i5vbZneeuutRbYm1Orfv38uZ2vEvhDlj8AbdVcVwSveM1rRmv7tRBiqioN3hO8IOVWxPWK92qAXoSP272hlif2ktq779u2by1HrD3/4Q7NtGP74xz/m5bG/0nUZDUGHsOyyy6bRo0c3WvbTn/40d3yM87tvv/12w2333XfPj0+cOLFh3U996lMN/z9jxoy83i677JJefvnlfL+1jjvuuDyKoOqzn/1s7jAYy6vivPO2226bX7u2rP369Ut77rlno7JGJ7YVV1yxUVnDxhtvnHbaaaeG+6uttlraYIMNGr3m4rj11ltz2dZdd918Pzo2fuELX8jLa6288sq5X0N0uluYWOd//ud/cse21rrvvvtyP4YxY8Y0LFtmmWXSl770pYU+56STTmp0P+qlpXr453/+51zHtdsmRP+XXr16NVr+0UcfpX/84x/5fvSRiT4DRxxxRKNtE9sx1m26bdpSptaaOXNmqzuZ3n///bn8Z555Zq67qqjTvn37pl/84heN1o/9q7YPUJ8+fdL222/fqLz33ntvGjhwYK6DqthOp59+eu5w+fDDDzd6zUMOOSTvk61R3Sa/+tWvcl8juiZhgQ5hjTXWyF9ytf72t7+l559/Pn9p1d7WX3/9/Pi0adMa1n300UfTHnvskYeixUEu1rvgggvyY20JC2uuuWaLX4RDhgxptvy9995rVNZ4n8985jPNyhtfxrVlbel9Qv/+/Ru9ZlvFATEOChGSXnzxxYbb5z//+dzj/q9//WvDuueff34+yMRBZb311ssH86jDWldccUV67rnn8meP9WLY3KIOmH//+9/ToEGDmnVOrYaXppZbbrlmB6WF1UNbtk2ovkZsmxAhs+m2idEHTbdNW8rUWnGQj1EprRF1GCI81oq/j7XXXrvh8arBgwc3CrgtlTeeE9u5NnzUjkJq+pprrbVWq8paXffss89O3/ve9/JIo+hUO2HChDb93dHxGQ1Bh1DbMlAVPbA322yzNH78+BafUz1IvPTSS2n48OG5BSLWjeXxxRoHziuvvLJNPbnj12Zrl0eLQ21ZIyg0/QVf1fTgs7D3qX3NtorWjblz56ZvfvOb+dZUlK063C8OEn/5y1/SPffck1sD7rjjjjys9KKLLmpY51/+5V/yL+o777wzH1S/8Y1vpMsvvzyPrthnn31Se1hYPbRl3UXVZXX7x3Dd+HXdVG2rRFvL1Fqxb8b8BNFi0DQUf1JLYl9q6e+xJPa3GI589913530lWiwuvfTSPPopwgydn7BAh7XOOuukZ555JgeBpr+casVcAnGQ/PnPf97o12dLzcul1/mkZY3m4/gV39Yv2oVpa1kjDMRsfuPGjWv22He/+930X//1Xw1BIEQrTEzaFLc4iB188MF5CN/YsWPzr+sQrQSnnHJKvsUv8K233jqvs7CwMHTo0FzvTYe+RgtHvcS2CRHmovWpHttmv/32S48//ngOZbWnAhZWhyHCXLQkVMU2euWVVxbrM8RrPvvsszk41bYu/PnPf270np/kM0ewj9tXvvKV9Nhjj+W/heuuuy5POEXn5zQEHVb8so3zzi2NP//ggw/S+++/3+iXVe0vqWgCvfHGG5s9Lw6QS2KmvCjr/Pnz08UXX9ziOPXFec8oa2jNcydPnpwn9YlyHHrooc1u0R8kDtjRByG88847jZ4fv3ajH0XU4ccff5w/S9Nm5DjYrr766jmYLUw0Qcfza7dZHKCiWbpeokxxGuDf//3fc9maeuutt9r8mtUg1NrtGn0gInh9+ctfbnQ6qCqCWPWgGmEgtsfVV1/daJ/+/ve/n7dJ9EFpq5gM6s0330y33XZbo/3ymmuuyaej4tTV4u6P0R8jXqtWhIYIJaV9hc5FywId1lFHHZV+8pOf5C/a+LUav1TiIBa/hmJ5dKiKznwjRozIX67x6+3EE0/MfQTiYBUHtzfeeKPRa0aHw+985zv5iznOo8c61Q6Tn0R82cZ7R9NrNDdHmaIDWZwvj9MD3/rWt/JBuy1ilrwIQtH0HweJ6AQaZY0yNxWtBnFg2X///Rd6sIjm9mh9iE59Ub5oko86HTBgQHrhhRfStddemw9E0REvDgjRfBxl3mKLLfIBJVpOokNkS6c4qg488MDcvyEOihFOovk9WnzefffdJdqyUxJBIbZ57E/RMnL44Yfn00KvvfZa7iwYdRCfvS2i9SjCVRx8ow/NKqusklt1FnadhuhDEKdzYjvEdq2dwfHpp5/OM27uuOOO+X6ULVp3ohVo7733zts0WhniNNF22223yAnNWnLCCSfk1qU4VfDUU0/lWRhvv/323E/lqquualXny2ihif5A0VoQ60d4iH0pWv9OPfXU3AE16iKCQ5zyiX03OkrSRdR7OAbdS2lSppbE8K4YEx+Px9DCGIIWQ9hiHPeMGTMa1vv5z3+eJ5mJYYIxfj6eE0PM4r1qx9HHcLmYECkmDGppUqYYHtbSePcYKtfS5D1NxfwFUb4YUx/vsdlmm+V5GqZMmdJsUqamoixNh3Vef/31efhnjLEvDaOM91lzzTUrJTFpT0zIE8MMY26FnXfeufLpT38612tMehUTTVXrNIbfxf0Y7hefIz5r/H/Tia5ampQp6irmaKhOyhTzX8Qw1yj/j3/840XWYbXOW5ooqFbURUtD9ha2LWP9GBobZYr9JD5zlC0mA2trmcJjjz2Wt3UMc2ztMMrYD2Keg/XXXz+XIeYQidf4+te/3mh/rg6VjDlBYm6EmJPj5JNPXuikTE21tF1iUqaYfyImmYoyxz7TdBjkwuq66u67785zMMR8ItVhlC+//HKekyLqMz7TKquskidIu//++xdZH3QePeI/9Q4sQNd211135emOf/e73+Vf8kDnIiwA7Sr6k9R28oxTR3HaI4Zvxnnz9uoACiw9+iwA7eq0007LgSHOwUcHtxhqGb3jo4OhoACdk5YFoF1FZ8voBBkdHD/88MPckfTkk0/OneCAzklYAACKzLMAABQJCwBAkbAAABQJCwBA9wwLMRd9TGkaF8SJKUmfeOKJehepQ4vLD8dUvLW3mKqXxuL6CzGtdFwjIeooJhuqFf2F48qNcR2AGCYY8/xXL5HcXS2qzmIK4qb7Xkxz3J3FtOExtXNMqxzTe8c02jHlc60YaRKXFv/0pz+dp+OOqZWnTp2aurPW1Nuuu+7abH+LKeXphmEh5muP66vH1fdi3vWY2z4uJtP0uvU0tskmm+RrKVRvMdsejcXFq2J/WtiFka644op8AaCYPz8u2hTz58e+F1/s3dWi6ixEOKjd9+JaCd3Zww8/nINAXOL5N7/5Tb4AVkxsVb14WjjrrLPyFVfj2iOx/pQpU/KVQ7uz1tRbGDNmTKP9Lf5uWYRKF7T99tvnaxBUzZ8/v7L66qtXLr300rqWqyOLee1j7n9aL/587rzzzob7CxYsqAwcOLDRvPrTp0/P11740Y9+VKdSduw6q17H4IADDqhbmTqDadOm5bp7+OGHG/aruGZE7XUxXnjhhbzO448/XseSdux6q15P44wzzqhruTqjLteyENd8j6uq1V7zPS6VGvfjevIsXDSXR1Px2muvnf71X/81X5WP1nvllVfydMa1+16/fv3yaTD7XtlDDz2Um4032GCDPIFT00tod3fVy4XH1S1DfMfFr+bafS1OG6655pr2tUK9VcXVV1ddddV8ldC4wuecOXPqVMLOo8tN9/z222/nuejjsru14n5c2piWxQHtpptuyl/W0SwXl8fdaaed0nPPPdeqy9eSclAILe171cdo+RRENJ+vtdZa6aWXXkoXXHBB2mefffJBLy5z3N0tWLAgnXnmmfkCXNVLYMf+FJdlj0tG17KvlestfPGLX0xDhw7NP4yeffbZdP755+d+DTEtOd0oLLB44su5avPNN8/hIf6gfvKTn6TjjjuurmWjazv88MMb/n+zzTbL+98666yTWxuGDx+eurs4Bx+hXR+i9qm3E044odH+Fp2RYz+LoBr7HS3rcqchomkpfo007RUc9wcOHFi3cnU28Ytl/fXXz/P70zrV/cu+98nEabD4O7bvpXw9jXvuuSdNnDgxDR48uGF57E9xynX69OmN1revleutJfHDKNjfullYiKa5bbbZJj3wwAONmqPiflwFj9aZPXt2TtqRummdaEaPL+rafW/mzJl5VIR9r/Vef/313GehO+970Rc0Dnh33nlnevDBB/O+VSu+43r37t1oX4um9Ohn1J33tUXVW0smTZqU/+3O+1u3PQ0RwyZHjRqVtt1227T99tunq666Kg+dGT16dL2L1mGdc845eSx8nHqIIVgx7DRaaI444oh6F63DhajaXyDRqTG+bKIDVXQui3Okl1xySVpvvfXyF9WFF16Yz43GeO/uqlRncYv+MTFHQAStCKjnnXdevlJlDDntzk3ocfXOu+++O/cZqvZDiA6zMX9H/BunB+O7Luqwb9+++dLgERR22GGH1F0tqt5i/4rHR44cmeeniD4LMQR15513zqe/KKh0Uddcc01lzTXXrPTp0ycPpfz9739f7yJ1aIcddlhl0KBBub7WWGONfP/FF1+sd7E6nIkTJ+ahWE1vMfyvOnzywgsvrAwYMCAPmRw+fHjlL3/5S6U7K9XZnDlzKiNGjKisttpqeSjg0KFDK2PGjKm8+eable6spfqK24033tiwzgcffFA55ZRTKv37968sv/zylYMOOqjyxhtvVLqzRdXba6+9Vtl5550rq6yySv77XHfddSvnnntuZcaMGfUueofnEtUAQPfqswAAtC9hAQAoEhYAgCJhAQAoEhYAgCJhAQAoEhYAgO4bFubOnZu++tWv5n9pHXW2eNRb26mzxaPe2k6dfXJdelKmmJc/pvmMa5rHdKgsmjpbPOqt7dTZ4lFvbafOPrku3bIAAHxywgIA0HWvOhmXno4rJMbVxXr06NFi01PtvyyaOls86q3t1NniUW9tp85aFr0QZs2ala+Mu8wyy3TdPgtx3fshQ4bUuxgA0GlNnjw5DR48uOu2LESLQvinNDL1Sr3rXRwA6DTmpY/T79K9DcfSLhsWqqceIij06iEsAECr/b/zCi2dxm9KB0cAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoEhYAACKhAUAoOOHhQkTJqRhw4al5ZZbLn32s59NTzzxRL2LBAB0lLBw2223pbPPPjuNGzcuPf3002mLLbZIe+21V5o2bVq9iwYAdISwMH78+DRmzJg0evTotPHGG6frrrsuLb/88umGG26od9EAgHqHhY8++ig99dRTaY899vj/BVpmmXz/8ccfb7b+3Llz08yZMxvdAIAuHBbefvvtNH/+/DRgwIBGy+P+m2++2Wz9Sy+9NPXr16/hNmTIkKVYWgDonup+GqItxo4dm2bMmNFwmzx5cr2LBABdXq96vvmqq66aevbsmaZOndpoedwfOHBgs/WXXXbZfAMAuknLQp8+fdI222yTHnjggYZlCxYsyPd33HHHehYNAOgILQshhk2OGjUqbbvttmn77bdPV111VXr//ffz6AgAoP7qHhYOO+yw9NZbb6WLLrood2rccsst03333des0yMAUB89KpVKJXVSMXQyRkXsmg5IvXr0rndxAKDTmFf5OD2U7s4DBvr27dt1RkMAAEufsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAFAkLAECRsAAAtH9YuPnmm9PcuXObLf/oo4/yYwBANw8Lo0ePTjNmzGi2fNasWfkxAKCbh4VKpZJ69OjRbPnrr7+e+vXr1x7lAgA6iF5tWXmrrbbKISFuw4cPT716/f+nz58/P73yyitp7733XhLlBAA6Q1g48MAD87+TJk1Ke+21V1pxxRUbHuvTp08aNmxYOuSQQ9q/lABA5wgL48aNyy0IEQpGjBiRBg0atORKBgB0zj4LPXv2TCeeeGL68MMPl0yJAIDO38Fx0003TS+//HL7lwYA6Bph4ZJLLknnnHNOuueee9Ibb7yRZs6c2egGAHTTPgtVI0eOzP/uv//+jYZQVodURr8GAKAbh4WJEye2f0kAgK4TFnbZZZf2LwkA0HXCQpg+fXr6/ve/n1544YV8f5NNNknHHnusGRwBoItZrA6OTz75ZFpnnXXSlVdemd599918Gz9+fF729NNPt38pAYC66VGJXolttNNOO6V11103XX/99Q1TPs+bNy8df/zxeUjlI488kpaGGHkRLRm7pgNSrx69l8p7AkBXMK/ycXoo3Z0vDNm3b9/2Pw0RLQu1QSG/UK9e6bzzzkvbbrvt4rwkANCVTkNEAnnttdeaLZ88eXJaaaWV2qNcAEBnDguHHXZYOu6449Jtt92WA0LcfvzjH+fTEEcccUT7lxIAqJvFOg3xH//xH3nypaOPPjr3VQi9e/dOJ598crrsssvau4wAQGfr4Fg1Z86c9NJLL+X/j5EQyy+/fFqadHAEgA7awbEqwsFmm232SV4CAOjgFissvP/++/l0wwMPPJCmTZuWFixY0OhxV6QEqn41ZVK9iwC0YOasBan/+mnJhYXoyPjwww+no446Kg0aNKjRxaQAgK5lscLCL3/5y/SLX/wiff7zn2//EgEAnX/oZP/+/dMqq6zS/qUBALpGWLj44ovTRRddlEdDAABdW6tPQ2y11VaN+ia8+OKLacCAAWnYsGF5joVaLiYFAN0wLBx44IFLtiQAQOcOC+PGjVuyJQEAOqRPNCnTU089lV544YX8/5tsskk+VQEAdC2LFRZiIqbDDz88PfTQQ2nllVfOy6ZPn5522223fEGp1VZbrb3LCQB0ptEQp512Wpo1a1Z6/vnn07vvvptvzz33XL5Ww+mnn97+pQQAOlfLwn333Zfuv//+tNFGGzUs23jjjdOECRPSiBEj2rN8AEBnbFmIa0E0HS4ZYlnT60QAAN0wLOy+++7pjDPOSFOmTGlY9o9//COdddZZafjw4e1ZPgCgM4aFa6+9NvdPiAmZ1llnnXxba6218rJrrrmm/UsJAHSuPgtDhgzJszRGv4U///nPeVn0X9hjjz3au3wAQGdqWXjwwQdzR8ZoQYipn/fcc888MiJu2223XZ5r4be//e2SKy0A0LHDwlVXXZXGjBmT+vbt2+yxfv36pRNPPDGNHz++PcsHAHSmsPDMM8+kvffee6GPx7DJmNURAOimYWHq1KktDpms6tWrV3rrrbfao1wAQGcMC2ussUaeqXFhnn322TRo0KD2KBcA0BnDwsiRI9OFF16YPvzww2aPffDBB/nKlPvuu297lg8AqLMelUql0pbTEFtvvXXq2bNnOvXUU9MGG2yQl8fwyZjqef78+XlI5YABA9LSEKMyomPlrumA1KvHwk+PAPXzqymT6l0EoAUzZy1I/dd/Oc2YMaPFgQuLPc9ChIDHHnssnXzyyWns2LGpmjNiGOVee+2VA8PSCgoAQAedlGno0KHp3nvvTe+991568cUXc2BYb731Uv/+/ZdMCQGAzjeDY4hwEBMxAQBd22JdGwIA6D6EBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCgSFgAAIqEBQCg44aFRx55JO23335p9dVXTz169Eh33XVXPYsDAHS0sPD++++nLbbYIk2YMKGexQAACnqlOtpnn33yDQDouOoaFtpq7ty5+VY1c+bMupYHALqDTtXB8dJLL039+vVruA0ZMqTeRQKALq9ThYWxY8emGTNmNNwmT55c7yIBQJfXqU5DLLvssvkGACw9naplAQDoZi0Ls2fPTi+++GLD/VdeeSVNmjQprbLKKmnNNdesZ9EAgI4QFp588sm02267Ndw/++yz87+jRo1KN910Ux1LBgB0iLCw6667pkqlUs8iAACLoM8CAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFAkLAAARcICAFDUK3VilUol/zsvfZzS//1foIOZOWtBvYsAtGDm7AWNjqVdNizMmjUr//u7dG+9iwIsRP/1610CYFHH0n79+hXX6VFpTaTooBYsWJCmTJmSVlpppdSjR49mj8+cOTMNGTIkTZ48OfXt27cuZexs1NniUW9tp84Wj3prO3XWsjj8R1BYffXV0zLLLNN1Wxbiww0ePHiR68XOYQdpG3W2eNRb26mzxaPe2k6dNbeoFoUqHRwBgCJhAQDovmFh2WWXTePGjcv/0jrqbPGot7ZTZ4tHvbWdOvvkOnUHRwBgyevSLQsAwCcnLAAARcICAFAkLAAARcIC0GrHHHNMni01br17904DBgxIe+65Z7rhhhvyjKpA1yQsAG2y9957pzfeeCO9+uqr6Ze//GXabbfd0hlnnJH23XffNG/evHoXD1gChAWgTWKs+sCBA9Maa6yRtt5663TBBReku+++OweHm266Ka8zfvz4tNlmm6UVVlghz8l/yimnpNmzZ+fH3n///Tzl7u23397ode+66668fvUCcUDHISwAn9juu++etthii/Szn/2s4botV199dXr++efTD37wg/Tggw+m8847Lz8WgeDwww9PN954Y6PXiPuHHnpovjAc0LGYlAloU5+F6dOn51aApiIAPPvss+lPf/pTs8eiFeGkk05Kb7/9dr7/xBNPpM997nP5KoCDBg1K06ZNyy0V999/f9pll12WymcBWk/LAtAu4ndH9VLxcdAfPnx4DgDRUnDUUUeld955J82ZMyc/vv3226dNNtkktzqEH/7wh2no0KFp5513rutnAFomLADt4oUXXkhrrbVW7vgYnR0333zzdMcdd6SnnnoqTZgwIa/z0UcfNax//PHHN/RxiFMQo0ePbggbQMciLACfWPRJ+N///d90yCGH5HAQwyi/+c1vph122CGtv/76acqUKc2ec+SRR6a///3vuW9DnLoYNWpUXcoOLFqvVqwD0GDu3LnpzTffTPPnz09Tp05N9913X7r00ktza8LRRx+dnnvuufTxxx+na665Ju23337p0UcfTdddd12z1+nfv386+OCD07nnnptGjBiRBg8eXJfPAyyalgWgTSIcRKfEYcOG5TkXJk6cmFsHYvhkz54986iIGDp5+eWXp0033TTdeuutOUy05LjjjsunJo499til/jmA1jMaAqibW265JZ111ln5NEWfPn3qXRxgIZyGAJa6GBURs0Bedtll6cQTTxQUoINzGgJY6q644oq04YYb5pkgx44dW+/iAIvgNAQAUKRlAQAoEhYAgCJhAQAoEhYAgCJhAQAoEhYAgCJhAQAoEhYAgCJhAQBIJf8H3zFUHXRFvNoAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pf.panelview(\n", + " data,\n", + " unit=\"user\",\n", + " time=\"day\",\n", + " treat=\"treat\",\n", + " collapse_to_cohort=True,\n", + " sort_by_timing=True,\n", + " ylab=\"Cohort\",\n", + " xlab=\"Day\",\n", + " title=\"Treatment Assignment Cohorts\",\n", + " figsize=(6, 5),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "a368b4ea", + "metadata": {}, + "source": [ + "We see that half of our users are treated after half the time. " + ] + }, + { + "cell_type": "markdown", + "id": "f9c7a8aa", + "metadata": {}, + "source": [ + "In the next step, we will look at three different ways to compute the average treatment effect of our feature on the outcome: \n", + "\n", + "- First, we will compute a standard \"difference in means\" estimator and conduct a two-sampled t-test for inference. We do this both by hand and by means of linear regression. While standard, we will show that this estimator is relatively inefficient \n", + " as it throws aways valuable information from the pre-treatment periods. \n", + "- In a second step, we improve on the difference in means estimator and include pre-treatment measures of the outcome variable to our regression model. In the tech blog world, this method \n", + " is often referred to as CUPED (which stands for \"Controlled-experiment Using Pre-Experiment Data\" as far as we know). We will demonstrate that CUPED uted leads to a significant reduction in the variance of our estimators. Pre-treatment averages are a \"good control\" variable because \n", + " a) under uncondional randomization, pre-treatment averages should be uncorrelated of the treatment assignment and b) pre-treatment averages should be highly predictive of the post-treatment outcome. The intuition here is simply that users will behave similarly \"after\" the experiment starts / when receiving the treatment than they did before. \n", + "- In a third step, we show that instead of using pre-treatment averages as a control variable, we can simply fit a panel model and control for user and time fixed effect. \n", + "\n", + "All three estimators identify the same average treatment effect, but the CUPED and panel estimator are much more efficient: they produce lower variances. \n", + "\n", + "But first, let's discuss why we are interested in reducing the variance of our estimators. " + ] + }, + { + "cell_type": "markdown", + "id": "d487603a", + "metadata": {}, + "source": [ + "## Variance of Estimators, Statistical Power and Sample Size Requirements\n", + "\n", + "In statistical experiments, we care about power to make sure that we detect a true effect. \n", + "\n", + "It depends on the **signal-to-noise ratio**:\n", + "\n", + "$$\n", + "\\text{Power} \\;\\sim\\; f\\!\\left(\\frac{|\\tau|}{\\text{SE}}\\right)\n", + "$$\n", + "\n", + "where\n", + "\n", + "- $\\tau$ = the true effect size \n", + "- $\\text{SE} = \\frac{\\sigma}{\\sqrt{n}}$ = the standard error of the estimate \n", + "- $\\sigma$ = standard deviation of the outcome \n", + "- $n$ = sample size per group (if balanced) \n", + "\n", + "So, anything that **increases the effect size**, **increases the sample size**, or **reduces outcome variance** will improve power. That's where our interest in variance reduction stems from. \n" + ] + }, + { + "cell_type": "markdown", + "id": "ef3befc5", + "metadata": {}, + "source": [ + "## Simple Difference in Means Estimator" + ] + }, + { + "cell_type": "markdown", + "id": "94686a5c", + "metadata": {}, + "source": [ + "The simplest way to analyse an AB test / estimate an average treatment effect is the difference in means estimator: \n", + "$$\n", + " \\tau = \\frac{1}{n_1} \\sum_{i=1}^{n_1} Y_{i,1} - \\frac{1}{n_0} \\sum_{i=1}^{n_0} Y_{i,0}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "b037d365", + "metadata": {}, + "source": [ + "We can compute it in a few lines of `pandas` by implementing three steps: \n", + "- First, we throw away all pre-experimental data. \n", + "- Then, we sum the post-treatment minutes watched into a single data point per user. \n", + "- Finally, we compute the difference of means of total minutes watched between the treated and control group.\n", + "Done! " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "5eeecd81", + "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", + "
usertotal_minutes_watched
treat
0NaNNaN
1-1387.6301510.761514
\n", + "
" + ], + "text/plain": [ + " user total_minutes_watched\n", + "treat \n", + "0 NaN NaN\n", + "1 -1387.630151 0.761514" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def _difference_in_means_pd(data):\n", + " # standard analyses: throw away pre-experimental data\n", + " data2 = data.copy()\n", + " data_post = data2[data2.day >= 15]\n", + " # collapse post-treatment minutes watched into a single data point per user\n", + " data_post_agg = (\n", + " data_post.groupby([\"user\", \"treat\"])\n", + " .agg({\"minutes_watched\": \"mean\"})\n", + " .reset_index()\n", + " .rename(columns={\"minutes_watched\": \"total_minutes_watched\"})\n", + " )\n", + " # compute difference of means estimator\n", + " return data_post_agg.groupby(\"treat\").mean().diff()\n", + "\n", + "\n", + "_difference_in_means_pd(data)" + ] + }, + { + "cell_type": "markdown", + "id": "c5d99260", + "metadata": {}, + "source": [ + "Because linear regression is just a fency way to compute and compare differences, we could also have estimated this via `pf.feols()`: " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "01495f78", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "###\n", + "\n", + "Estimation: OLS\n", + "Dep. var.: minutes_watched, Fixed effects: 0\n", + "Inference: iid\n", + "Observations: 1500000\n", + "\n", + "| Coefficient | Estimate | Std. Error | t value | Pr(>|t|) | 2.5% | 97.5% |\n", + "|:--------------|-----------:|-------------:|----------:|-----------:|-------:|--------:|\n", + "| Intercept | 0.357 | 0.015 | 24.163 | 0.000 | 0.328 | 0.386 |\n", + "| treat | 0.762 | 0.209 | 3.641 | 0.000 | 0.352 | 1.171 |\n", + "---\n", + "RMSE: 18.07 R2: 0.0 \n" + ] + } + ], + "source": [ + "def _difference_in_means(data):\n", + " data2 = data.copy()\n", + "\n", + " data_post = data2[data2.day >= 15]\n", + " fit = pf.feols(\"minutes_watched ~ treat\", data=data_post)\n", + "\n", + " return fit\n", + "\n", + "\n", + "fit_difference_in_means = _difference_in_means(data)\n", + "fit_difference_in_means.summary()" + ] + }, + { + "cell_type": "markdown", + "id": "c9f44820", + "metadata": {}, + "source": [ + "Note that no aggregation step is needed and that we get identical point estimates. As a bonus, we get standard errors and confidence intervals in the process. " + ] + }, + { + "cell_type": "markdown", + "id": "e972baa3", + "metadata": {}, + "source": [ + "## Cuped: Using Pre-Experimental Measures of the Outcomes Variable as Controls\n", + "\n", + "Sometimes, throwing away data can be a winning strategy (\"garbage in garbage out\"), but not in this example: we have high-quality pre-experimental measures of the behavior \n", + "of our users at hand, and we should use it in our favor! \n", + "\n", + "More specifically, instead of throwing away all of the pre-experimental measures, we could instead have used them as a control! \n", + "\n", + "This should help reduce the variance of our estimators because pre-experimental behavior is likely to be highly predicitive of what users do after the launch of the experiment. Consequently, including these baselines in our regression models should help us \"explain residual errors\" and thereby reduce the variance of our \n", + "estimators. \n", + "\n", + "If this works, it's great news, as it allows us to run the same experiment on a smaller number of users and still achieve the same level of power!" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "817bb864", + "metadata": {}, + "outputs": [], + "source": [ + "def _regression_cuped(data):\n", + " data = data.copy()\n", + " data[\"pre_experiment\"] = data[\"day\"] < 15\n", + "\n", + " # aggregate pre-post averages by user\n", + " agg = data.groupby([\"user\", \"pre_experiment\"], as_index=False).agg(\n", + " minutes_watched=(\"minutes_watched\", \"mean\")\n", + " )\n", + "\n", + " wide = (\n", + " agg.pivot(index=\"user\", columns=\"pre_experiment\", values=\"minutes_watched\")\n", + " .rename(columns={True: \"minutes_pre\", False: \"minutes_post\"})\n", + " .reset_index()\n", + " )\n", + "\n", + " wide = wide.merge(\n", + " data[[\"user\", \"ever_treated\"]].drop_duplicates(), on=\"user\", how=\"left\"\n", + " ).rename(columns={\"ever_treated\": \"treat\"})\n", + "\n", + " # center the pre metric\n", + " mu_pre = wide[\"minutes_pre\"].mean()\n", + " wide[\"minutes_pre_c\"] = wide[\"minutes_pre\"] - mu_pre\n", + "\n", + " fit_cuped = pf.feols(\n", + " \"minutes_post ~ treat + minutes_pre_c\", data=wide, vcov=\"hetero\"\n", + " )\n", + "\n", + " return fit_cuped, wide\n", + "\n", + "\n", + "fit_cuped, data_cuped = _regression_cuped(data)" + ] + }, + { + "cell_type": "markdown", + "id": "9ae49f77", + "metadata": {}, + "source": [ + "Per user, we now log the total minutes watched before and after the treatment and then fit a regression model of the following form: ^^" + ] + }, + { + "cell_type": "markdown", + "id": "6df64964", + "metadata": {}, + "source": [ + "$$\n", + " \\text{total minutes watched after treatment} = \\alpha + \\beta \\text{treat} + \\gamma (\\text{total minutes watched before treatment} - \\text{avg total minutes watched before treatment}) + \\epsilon\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "550da9fc", + "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", + "
userminutes_postminutes_pretreatminutes_pre_c
007.7478417.76717807.894779
1125.59592824.674535024.802136
22-32.224392-32.1303340-32.002733
33-12.269762-13.1424350-13.014835
44-1.448348-1.3920430-1.264442
\n", + "
" + ], + "text/plain": [ + " user minutes_post minutes_pre treat minutes_pre_c\n", + "0 0 7.747841 7.767178 0 7.894779\n", + "1 1 25.595928 24.674535 0 24.802136\n", + "2 2 -32.224392 -32.130334 0 -32.002733\n", + "3 3 -12.269762 -13.142435 0 -13.014835\n", + "4 4 -1.448348 -1.392043 0 -1.264442" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "data_cuped.head()" + ] + }, + { + "cell_type": "markdown", + "id": "56861ab2", + "metadata": {}, + "source": [ + "We can compare with the previous results: " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "67d8c3f6", + "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", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "\n", + "\n", + "
\n", + " minutes_watched\n", + " \n", + " minutes_post\n", + "
(1)(2)
coef
treat0.762***
(0.209)		0.192***
(0.029)
minutes_pre_c0.999***
(0.000)
Intercept0.357***
(0.015)		0.360***
(0.002)
stats
Observations1500000100000
S.E. typeiidhetero
R20.0000.999
Adj. R20.0000.999
Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001. Format of coefficient cell:\n", + "Coefficient \n", + " (Std. Error)
\n", + "\n", + "
\n", + " " + ], + "text/plain": [ + "GT(_tbl_data= level_0 level_1 0 1\n", + "0 coef treat 0.762***
(0.209) 0.192***
(0.029)\n", + "1 coef minutes_pre_c 0.999***
(0.000)\n", + "2 coef Intercept 0.357***
(0.015) 0.360***
(0.002)\n", + "3 stats Observations 1500000 100000\n", + "4 stats S.E. type iid hetero\n", + "5 stats R2 0.000 0.999\n", + "6 stats Adj. R2 0.000 0.999, _body=, _boxhead=Boxhead([ColInfo(var='level_0', type=, column_label='level_0', column_align='center', column_width=None), ColInfo(var='level_1', type=, column_label='level_1', column_align='center', column_width=None), ColInfo(var='0', type=, column_label='(1)', column_align='center', column_width=None), ColInfo(var='1', type=, column_label='(2)', column_align='center', column_width=None)]), _stub=, _spanners=Spanners([SpannerInfo(spanner_id='minutes_watched', spanner_level=1, spanner_label='minutes_watched', spanner_units=None, spanner_pattern=None, vars=['0'], built=None), SpannerInfo(spanner_id='minutes_post', spanner_level=1, spanner_label='minutes_post', spanner_units=None, spanner_pattern=None, vars=['1'], built=None)]), _heading=Heading(title=None, subtitle=None, preheader=None), _stubhead=None, _source_notes=['Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001. Format of coefficient cell:\\nCoefficient \\n (Std. Error)'], _footnotes=[], _styles=[], _locale=, _formats=[], _substitutions=[], _options=Options(table_id=OptionsInfo(scss=False, category='table', type='value', value=None), table_caption=OptionsInfo(scss=False, category='table', type='value', value=None), table_width=OptionsInfo(scss=True, category='table', type='px', value='auto'), table_layout=OptionsInfo(scss=True, category='table', type='value', value='fixed'), table_margin_left=OptionsInfo(scss=True, category='table', type='px', value='auto'), table_margin_right=OptionsInfo(scss=True, category='table', type='px', value='auto'), table_background_color=OptionsInfo(scss=True, category='table', type='value', value='#FFFFFF'), table_additional_css=OptionsInfo(scss=False, category='table', type='values', value=[]), table_font_names=OptionsInfo(scss=False, category='table', type='values', value=['-apple-system', 'BlinkMacSystemFont', 'Segoe UI', 'Roboto', 'Oxygen', 'Ubuntu', 'Cantarell', 'Helvetica Neue', 'Fira Sans', 'Droid Sans', 'Arial', 'sans-serif']), table_font_size=OptionsInfo(scss=True, category='table', type='px', value='16px'), table_font_weight=OptionsInfo(scss=True, category='table', type='value', value='normal'), table_font_style=OptionsInfo(scss=True, category='table', type='value', value='normal'), table_font_color=OptionsInfo(scss=True, category='table', type='value', value='#333333'), table_font_color_light=OptionsInfo(scss=True, category='table', type='value', value='#FFFFFF'), table_border_top_include=OptionsInfo(scss=False, category='table', type='boolean', value=True), table_border_top_style=OptionsInfo(scss=True, category='table', type='value', value='solid'), table_border_top_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_top_color=OptionsInfo(scss=True, category='table', type='value', value='#A8A8A8'), table_border_right_style=OptionsInfo(scss=True, category='table', type='value', value='none'), table_border_right_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_right_color=OptionsInfo(scss=True, category='table', type='value', value='#D3D3D3'), table_border_bottom_include=OptionsInfo(scss=False, category='table', type='boolean', value=True), table_border_bottom_style=OptionsInfo(scss=True, category='table', type='value', value='hidden'), table_border_bottom_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_bottom_color=OptionsInfo(scss=True, category='table', type='value', value='#A8A8A8'), table_border_left_style=OptionsInfo(scss=True, category='table', type='value', value='none'), table_border_left_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_left_color=OptionsInfo(scss=True, category='table', type='value', value='#D3D3D3'), heading_background_color=OptionsInfo(scss=True, category='heading', type='value', value=None), heading_align=OptionsInfo(scss=True, category='heading', type='value', value='center'), heading_title_font_size=OptionsInfo(scss=True, category='heading', type='px', value='125%'), heading_title_font_weight=OptionsInfo(scss=True, category='heading', type='value', value='initial'), heading_subtitle_font_size=OptionsInfo(scss=True, category='heading', type='px', value='85%'), heading_subtitle_font_weight=OptionsInfo(scss=True, category='heading', type='value', value='initial'), heading_padding=OptionsInfo(scss=True, category='heading', type='px', value='4px'), heading_padding_horizontal=OptionsInfo(scss=True, category='heading', type='px', value='5px'), heading_border_bottom_style=OptionsInfo(scss=True, category='heading', type='value', value='solid'), heading_border_bottom_width=OptionsInfo(scss=True, category='heading', type='px', value='2px'), heading_border_bottom_color=OptionsInfo(scss=True, category='heading', type='value', value='#D3D3D3'), heading_border_lr_style=OptionsInfo(scss=True, category='heading', type='value', value='none'), heading_border_lr_width=OptionsInfo(scss=True, category='heading', type='px', value='1px'), heading_border_lr_color=OptionsInfo(scss=True, category='heading', type='value', value='#D3D3D3'), column_labels_background_color=OptionsInfo(scss=True, category='column_labels', type='value', value=None), column_labels_font_size=OptionsInfo(scss=True, category='column_labels', type='px', value='100%'), column_labels_font_weight=OptionsInfo(scss=True, category='column_labels', type='value', value='normal'), column_labels_text_transform=OptionsInfo(scss=True, category='column_labels', type='value', value='inherit'), column_labels_padding=OptionsInfo(scss=True, category='column_labels', type='px', value='4px'), column_labels_padding_horizontal=OptionsInfo(scss=True, category='column_labels', type='px', value='5px'), column_labels_vlines_style=OptionsInfo(scss=True, category='table_body', type='value', value='none'), column_labels_vlines_width=OptionsInfo(scss=True, category='table_body', type='px', value='0px'), column_labels_vlines_color=OptionsInfo(scss=True, category='table_body', type='value', value='white'), column_labels_border_top_style=OptionsInfo(scss=True, category='column_labels', type='value', value='solid'), column_labels_border_top_width=OptionsInfo(scss=True, category='column_labels', type='px', value='2px'), column_labels_border_top_color=OptionsInfo(scss=True, category='column_labels', type='value', value='black'), column_labels_border_bottom_style=OptionsInfo(scss=True, category='column_labels', type='value', value='solid'), column_labels_border_bottom_width=OptionsInfo(scss=True, category='column_labels', type='px', value='0.5px'), column_labels_border_bottom_color=OptionsInfo(scss=True, category='column_labels', type='value', value='black'), column_labels_border_lr_style=OptionsInfo(scss=True, category='column_labels', type='value', value='none'), column_labels_border_lr_width=OptionsInfo(scss=True, category='column_labels', type='px', value='1px'), column_labels_border_lr_color=OptionsInfo(scss=True, category='column_labels', type='value', value='#D3D3D3'), column_labels_hidden=OptionsInfo(scss=False, category='column_labels', type='boolean', value=False), row_group_background_color=OptionsInfo(scss=True, category='row_group', type='value', value=None), row_group_font_size=OptionsInfo(scss=True, category='row_group', type='px', value='0px'), row_group_font_weight=OptionsInfo(scss=True, category='row_group', type='value', value='initial'), row_group_text_transform=OptionsInfo(scss=True, category='row_group', type='value', value='inherit'), row_group_padding=OptionsInfo(scss=True, category='row_group', type='px', value='0px'), row_group_padding_horizontal=OptionsInfo(scss=True, category='row_group', type='px', value='5px'), row_group_border_top_style=OptionsInfo(scss=True, category='row_group', type='value', value='solid'), row_group_border_top_width=OptionsInfo(scss=True, category='row_group', type='px', value='0.5px'), row_group_border_top_color=OptionsInfo(scss=True, category='row_group', type='value', value='black'), row_group_border_right_style=OptionsInfo(scss=True, category='row_group', type='value', value='none'), row_group_border_right_width=OptionsInfo(scss=True, category='row_group', type='px', value='1px'), row_group_border_right_color=OptionsInfo(scss=True, category='row_group', type='value', value='white'), row_group_border_bottom_style=OptionsInfo(scss=True, category='row_group', type='value', value='solid'), row_group_border_bottom_width=OptionsInfo(scss=True, category='row_group', type='px', value='0.5px'), row_group_border_bottom_color=OptionsInfo(scss=True, category='row_group', type='value', value='black'), row_group_border_left_style=OptionsInfo(scss=True, category='row_group', type='value', value='none'), row_group_border_left_width=OptionsInfo(scss=True, category='row_group', type='px', value='1px'), row_group_border_left_color=OptionsInfo(scss=True, category='row_group', type='value', value='white'), row_group_as_column=OptionsInfo(scss=False, category='row_group', type='boolean', value=False), table_body_hlines_style=OptionsInfo(scss=True, category='table_body', type='value', value='none'), table_body_hlines_width=OptionsInfo(scss=True, category='table_body', type='px', value='1px'), table_body_hlines_color=OptionsInfo(scss=True, category='table_body', type='value', value='#D3D3D3'), table_body_vlines_style=OptionsInfo(scss=True, category='table_body', type='value', value='none'), table_body_vlines_width=OptionsInfo(scss=True, category='table_body', type='px', value='0px'), table_body_vlines_color=OptionsInfo(scss=True, category='table_body', type='value', value='white'), table_body_border_top_style=OptionsInfo(scss=True, category='table_body', type='value', value='solid'), table_body_border_top_width=OptionsInfo(scss=True, category='table_body', type='px', value='0.5px'), table_body_border_top_color=OptionsInfo(scss=True, category='table_body', type='value', value='black'), table_body_border_bottom_style=OptionsInfo(scss=True, category='table_body', type='value', value='solid'), table_body_border_bottom_width=OptionsInfo(scss=True, category='table_body', type='px', value='2px'), table_body_border_bottom_color=OptionsInfo(scss=True, category='table_body', type='value', value='black'), data_row_padding=OptionsInfo(scss=True, category='data_row', type='px', value='4px'), data_row_padding_horizontal=OptionsInfo(scss=True, category='data_row', type='px', value='5px'), stub_background_color=OptionsInfo(scss=True, category='stub', type='value', value=None), stub_font_size=OptionsInfo(scss=True, category='stub', type='px', value='100%'), stub_font_weight=OptionsInfo(scss=True, category='stub', type='value', value='initial'), stub_text_transform=OptionsInfo(scss=True, category='stub', type='value', value='inherit'), stub_border_style=OptionsInfo(scss=True, category='stub', type='value', value='hidden'), stub_border_width=OptionsInfo(scss=True, category='stub', type='px', value='2px'), stub_border_color=OptionsInfo(scss=True, category='stub', type='value', value='#D3D3D3'), stub_row_group_background_color=OptionsInfo(scss=True, category='stub', type='value', value=None), stub_row_group_font_size=OptionsInfo(scss=True, category='stub', type='px', value='100%'), stub_row_group_font_weight=OptionsInfo(scss=True, category='stub', type='value', value='initial'), stub_row_group_text_transform=OptionsInfo(scss=True, category='stub', type='value', value='inherit'), stub_row_group_border_style=OptionsInfo(scss=True, category='stub', type='value', value='solid'), stub_row_group_border_width=OptionsInfo(scss=True, category='stub', type='px', value='2px'), stub_row_group_border_color=OptionsInfo(scss=True, category='stub', type='value', value='#D3D3D3'), source_notes_padding=OptionsInfo(scss=True, category='source_notes', type='px', value='4px'), source_notes_padding_horizontal=OptionsInfo(scss=True, category='source_notes', type='px', value='5px'), source_notes_background_color=OptionsInfo(scss=True, category='source_notes', type='value', value=None), source_notes_font_size=OptionsInfo(scss=True, category='source_notes', type='px', value='90%'), source_notes_border_bottom_style=OptionsInfo(scss=True, category='source_notes', type='value', value='none'), source_notes_border_bottom_width=OptionsInfo(scss=True, category='source_notes', type='px', value='2px'), source_notes_border_bottom_color=OptionsInfo(scss=True, category='source_notes', type='value', value='#D3D3D3'), source_notes_border_lr_style=OptionsInfo(scss=True, category='source_notes', type='value', value='none'), source_notes_border_lr_width=OptionsInfo(scss=True, category='source_notes', type='px', value='2px'), source_notes_border_lr_color=OptionsInfo(scss=True, category='source_notes', type='value', value='#D3D3D3'), source_notes_multiline=OptionsInfo(scss=False, category='source_notes', type='boolean', value=True), source_notes_sep=OptionsInfo(scss=False, category='source_notes', type='value', value=' '), row_striping_background_color=OptionsInfo(scss=True, category='row', type='value', value='rgba(128,128,128,0.05)'), row_striping_include_stub=OptionsInfo(scss=False, category='row', type='boolean', value=False), row_striping_include_table_body=OptionsInfo(scss=False, category='row', type='boolean', value=False), container_width=OptionsInfo(scss=False, category='container', type='px', value='auto'), container_height=OptionsInfo(scss=False, category='container', type='px', value='auto'), container_padding_x=OptionsInfo(scss=False, category='container', type='px', value='0px'), container_padding_y=OptionsInfo(scss=False, category='container', type='px', value='10px'), container_overflow_x=OptionsInfo(scss=False, category='container', type='overflow', value='auto'), container_overflow_y=OptionsInfo(scss=False, category='container', type='overflow', value='auto'), quarto_disable_processing=OptionsInfo(scss=False, category='quarto', type='logical', value=False), quarto_use_bootstrap=OptionsInfo(scss=False, category='quarto', type='logical', value=False)), _has_built=False)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pf.etable([fit_difference_in_means, fit_cuped])" + ] + }, + { + "cell_type": "markdown", + "id": "34dc9d20", + "metadata": {}, + "source": [ + "Point estimates for the difference-in-means estimator is 0.762, while it is 0.192 for the Cuped estimator. But most importantly, the standard errors are much smaller for the CUPED estimator (0.03 vs 0.209). Because we have fitted a regression model with covariates, we have used heteroskedasticity-robust standard errors, which should be more conservative than the iid errors used in the difference-in-means regression." + ] + }, + { + "cell_type": "markdown", + "id": "2899a0fe", + "metadata": {}, + "source": [ + "## Panel Estimator\n", + "\n", + "Instead of collapsing all pre-and post information into a single average (and thereby losing information), we can as well be lazy and just use a panel estimator.\n", + "\n", + "Via `pyfixest`, that's one line of code: " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "0e6998a8", + "metadata": {}, + "outputs": [], + "source": [ + "fit_panel = pf.feols(\n", + " \"minutes_watched ~ treat | user + day\",\n", + " data=data,\n", + " vcov=\"hetero\",\n", + " demeaner_backend=\"rust\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "52670ecf", + "metadata": {}, + "source": [ + "We can compare all results: " + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "8cf20732", + "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", + " \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", + " minutes_watched\n", + " \n", + " minutes_post\n", + " \n", + " minutes_watched\n", + "
(1)(2)(3)
coef
treat0.762***
(0.209)		0.192***
(0.029)		0.191***
(0.012)
minutes_pre_c0.999***
(0.000)
Intercept0.357***
(0.015)		0.360***
(0.002)
fe
day--x
user--x
stats
Observations15000001000003000000
S.E. typeiidheterohetero
R20.0000.9990.999
R2 Within--0.000
Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001. Format of coefficient cell:\n", + "Coefficient \n", + " (Std. Error)
\n", + "\n", + "
\n", + " " + ], + "text/plain": [ + "GT(_tbl_data= level_0 level_1 0 1 \\\n", + "0 coef treat 0.762***
(0.209) 0.192***
(0.029) \n", + "1 coef minutes_pre_c 0.999***
(0.000) \n", + "2 coef Intercept 0.357***
(0.015) 0.360***
(0.002) \n", + "3 fe day - - \n", + "4 fe user - - \n", + "5 stats Observations 1500000 100000 \n", + "6 stats S.E. type iid hetero \n", + "7 stats R2 0.000 0.999 \n", + "8 stats R2 Within - - \n", + "\n", + " 2 \n", + "0 0.191***
(0.012) \n", + "1 \n", + "2 \n", + "3 x \n", + "4 x \n", + "5 3000000 \n", + "6 hetero \n", + "7 0.999 \n", + "8 0.000 , _body=, _boxhead=Boxhead([ColInfo(var='level_0', type=, column_label='level_0', column_align='center', column_width=None), ColInfo(var='level_1', type=, column_label='level_1', column_align='center', column_width=None), ColInfo(var='0', type=, column_label='(1)', column_align='center', column_width=None), ColInfo(var='1', type=, column_label='(2)', column_align='center', column_width=None), ColInfo(var='2', type=, column_label='(3)', column_align='center', column_width=None)]), _stub=, _spanners=Spanners([SpannerInfo(spanner_id='minutes_watched', spanner_level=1, spanner_label='minutes_watched', spanner_units=None, spanner_pattern=None, vars=['0', '2'], built=None), SpannerInfo(spanner_id='minutes_post', spanner_level=1, spanner_label='minutes_post', spanner_units=None, spanner_pattern=None, vars=['1'], built=None)]), _heading=Heading(title=None, subtitle=None, preheader=None), _stubhead=None, _source_notes=['Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001. Format of coefficient cell:\\nCoefficient \\n (Std. Error)'], _footnotes=[], _styles=[], _locale=, _formats=[], _substitutions=[], _options=Options(table_id=OptionsInfo(scss=False, category='table', type='value', value=None), table_caption=OptionsInfo(scss=False, category='table', type='value', value=None), table_width=OptionsInfo(scss=True, category='table', type='px', value='auto'), table_layout=OptionsInfo(scss=True, category='table', type='value', value='fixed'), table_margin_left=OptionsInfo(scss=True, category='table', type='px', value='auto'), table_margin_right=OptionsInfo(scss=True, category='table', type='px', value='auto'), table_background_color=OptionsInfo(scss=True, category='table', type='value', value='#FFFFFF'), table_additional_css=OptionsInfo(scss=False, category='table', type='values', value=[]), table_font_names=OptionsInfo(scss=False, category='table', type='values', value=['-apple-system', 'BlinkMacSystemFont', 'Segoe UI', 'Roboto', 'Oxygen', 'Ubuntu', 'Cantarell', 'Helvetica Neue', 'Fira Sans', 'Droid Sans', 'Arial', 'sans-serif']), table_font_size=OptionsInfo(scss=True, category='table', type='px', value='16px'), table_font_weight=OptionsInfo(scss=True, category='table', type='value', value='normal'), table_font_style=OptionsInfo(scss=True, category='table', type='value', value='normal'), table_font_color=OptionsInfo(scss=True, category='table', type='value', value='#333333'), table_font_color_light=OptionsInfo(scss=True, category='table', type='value', value='#FFFFFF'), table_border_top_include=OptionsInfo(scss=False, category='table', type='boolean', value=True), table_border_top_style=OptionsInfo(scss=True, category='table', type='value', value='solid'), table_border_top_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_top_color=OptionsInfo(scss=True, category='table', type='value', value='#A8A8A8'), table_border_right_style=OptionsInfo(scss=True, category='table', type='value', value='none'), table_border_right_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_right_color=OptionsInfo(scss=True, category='table', type='value', value='#D3D3D3'), table_border_bottom_include=OptionsInfo(scss=False, category='table', type='boolean', value=True), table_border_bottom_style=OptionsInfo(scss=True, category='table', type='value', value='hidden'), table_border_bottom_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_bottom_color=OptionsInfo(scss=True, category='table', type='value', value='#A8A8A8'), table_border_left_style=OptionsInfo(scss=True, category='table', type='value', value='none'), table_border_left_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_left_color=OptionsInfo(scss=True, category='table', type='value', value='#D3D3D3'), heading_background_color=OptionsInfo(scss=True, category='heading', type='value', value=None), heading_align=OptionsInfo(scss=True, category='heading', type='value', value='center'), heading_title_font_size=OptionsInfo(scss=True, category='heading', type='px', value='125%'), heading_title_font_weight=OptionsInfo(scss=True, category='heading', type='value', value='initial'), heading_subtitle_font_size=OptionsInfo(scss=True, category='heading', type='px', value='85%'), heading_subtitle_font_weight=OptionsInfo(scss=True, category='heading', type='value', value='initial'), heading_padding=OptionsInfo(scss=True, category='heading', type='px', value='4px'), heading_padding_horizontal=OptionsInfo(scss=True, category='heading', type='px', value='5px'), heading_border_bottom_style=OptionsInfo(scss=True, category='heading', type='value', value='solid'), heading_border_bottom_width=OptionsInfo(scss=True, category='heading', type='px', value='2px'), heading_border_bottom_color=OptionsInfo(scss=True, category='heading', type='value', value='#D3D3D3'), heading_border_lr_style=OptionsInfo(scss=True, category='heading', type='value', value='none'), heading_border_lr_width=OptionsInfo(scss=True, category='heading', type='px', value='1px'), heading_border_lr_color=OptionsInfo(scss=True, category='heading', type='value', value='#D3D3D3'), column_labels_background_color=OptionsInfo(scss=True, category='column_labels', type='value', value=None), column_labels_font_size=OptionsInfo(scss=True, category='column_labels', type='px', value='100%'), column_labels_font_weight=OptionsInfo(scss=True, category='column_labels', type='value', value='normal'), column_labels_text_transform=OptionsInfo(scss=True, category='column_labels', type='value', value='inherit'), column_labels_padding=OptionsInfo(scss=True, category='column_labels', type='px', value='4px'), column_labels_padding_horizontal=OptionsInfo(scss=True, category='column_labels', type='px', value='5px'), column_labels_vlines_style=OptionsInfo(scss=True, category='table_body', type='value', value='none'), column_labels_vlines_width=OptionsInfo(scss=True, category='table_body', type='px', value='0px'), column_labels_vlines_color=OptionsInfo(scss=True, category='table_body', type='value', value='white'), column_labels_border_top_style=OptionsInfo(scss=True, category='column_labels', type='value', value='solid'), column_labels_border_top_width=OptionsInfo(scss=True, category='column_labels', type='px', value='2px'), column_labels_border_top_color=OptionsInfo(scss=True, category='column_labels', type='value', value='black'), column_labels_border_bottom_style=OptionsInfo(scss=True, category='column_labels', type='value', value='solid'), column_labels_border_bottom_width=OptionsInfo(scss=True, category='column_labels', type='px', value='0.5px'), column_labels_border_bottom_color=OptionsInfo(scss=True, category='column_labels', type='value', value='black'), column_labels_border_lr_style=OptionsInfo(scss=True, category='column_labels', type='value', value='none'), column_labels_border_lr_width=OptionsInfo(scss=True, category='column_labels', type='px', value='1px'), column_labels_border_lr_color=OptionsInfo(scss=True, category='column_labels', type='value', value='#D3D3D3'), column_labels_hidden=OptionsInfo(scss=False, category='column_labels', type='boolean', value=False), row_group_background_color=OptionsInfo(scss=True, category='row_group', type='value', value=None), row_group_font_size=OptionsInfo(scss=True, category='row_group', type='px', value='0px'), row_group_font_weight=OptionsInfo(scss=True, category='row_group', type='value', value='initial'), row_group_text_transform=OptionsInfo(scss=True, category='row_group', type='value', value='inherit'), row_group_padding=OptionsInfo(scss=True, category='row_group', type='px', value='0px'), row_group_padding_horizontal=OptionsInfo(scss=True, category='row_group', type='px', value='5px'), row_group_border_top_style=OptionsInfo(scss=True, category='row_group', type='value', value='solid'), row_group_border_top_width=OptionsInfo(scss=True, category='row_group', type='px', value='0.5px'), row_group_border_top_color=OptionsInfo(scss=True, category='row_group', type='value', value='black'), row_group_border_right_style=OptionsInfo(scss=True, category='row_group', type='value', value='none'), row_group_border_right_width=OptionsInfo(scss=True, category='row_group', type='px', value='1px'), row_group_border_right_color=OptionsInfo(scss=True, category='row_group', type='value', value='white'), row_group_border_bottom_style=OptionsInfo(scss=True, category='row_group', type='value', value='solid'), row_group_border_bottom_width=OptionsInfo(scss=True, category='row_group', type='px', value='0.5px'), row_group_border_bottom_color=OptionsInfo(scss=True, category='row_group', type='value', value='black'), row_group_border_left_style=OptionsInfo(scss=True, category='row_group', type='value', value='none'), row_group_border_left_width=OptionsInfo(scss=True, category='row_group', type='px', value='1px'), row_group_border_left_color=OptionsInfo(scss=True, category='row_group', type='value', value='white'), row_group_as_column=OptionsInfo(scss=False, category='row_group', type='boolean', value=False), table_body_hlines_style=OptionsInfo(scss=True, category='table_body', type='value', value='none'), table_body_hlines_width=OptionsInfo(scss=True, category='table_body', type='px', value='1px'), table_body_hlines_color=OptionsInfo(scss=True, category='table_body', type='value', value='#D3D3D3'), table_body_vlines_style=OptionsInfo(scss=True, category='table_body', type='value', value='none'), table_body_vlines_width=OptionsInfo(scss=True, category='table_body', type='px', value='0px'), table_body_vlines_color=OptionsInfo(scss=True, category='table_body', type='value', value='white'), table_body_border_top_style=OptionsInfo(scss=True, category='table_body', type='value', value='solid'), table_body_border_top_width=OptionsInfo(scss=True, category='table_body', type='px', value='0.5px'), table_body_border_top_color=OptionsInfo(scss=True, category='table_body', type='value', value='black'), table_body_border_bottom_style=OptionsInfo(scss=True, category='table_body', type='value', value='solid'), table_body_border_bottom_width=OptionsInfo(scss=True, category='table_body', type='px', value='2px'), table_body_border_bottom_color=OptionsInfo(scss=True, category='table_body', type='value', value='black'), data_row_padding=OptionsInfo(scss=True, category='data_row', type='px', value='4px'), data_row_padding_horizontal=OptionsInfo(scss=True, category='data_row', type='px', value='5px'), stub_background_color=OptionsInfo(scss=True, category='stub', type='value', value=None), stub_font_size=OptionsInfo(scss=True, category='stub', type='px', value='100%'), stub_font_weight=OptionsInfo(scss=True, category='stub', type='value', value='initial'), stub_text_transform=OptionsInfo(scss=True, category='stub', type='value', value='inherit'), stub_border_style=OptionsInfo(scss=True, category='stub', type='value', value='hidden'), stub_border_width=OptionsInfo(scss=True, category='stub', type='px', value='2px'), stub_border_color=OptionsInfo(scss=True, category='stub', type='value', value='#D3D3D3'), stub_row_group_background_color=OptionsInfo(scss=True, category='stub', type='value', value=None), stub_row_group_font_size=OptionsInfo(scss=True, category='stub', type='px', value='100%'), stub_row_group_font_weight=OptionsInfo(scss=True, category='stub', type='value', value='initial'), stub_row_group_text_transform=OptionsInfo(scss=True, category='stub', type='value', value='inherit'), stub_row_group_border_style=OptionsInfo(scss=True, category='stub', type='value', value='solid'), stub_row_group_border_width=OptionsInfo(scss=True, category='stub', type='px', value='2px'), stub_row_group_border_color=OptionsInfo(scss=True, category='stub', type='value', value='#D3D3D3'), source_notes_padding=OptionsInfo(scss=True, category='source_notes', type='px', value='4px'), source_notes_padding_horizontal=OptionsInfo(scss=True, category='source_notes', type='px', value='5px'), source_notes_background_color=OptionsInfo(scss=True, category='source_notes', type='value', value=None), source_notes_font_size=OptionsInfo(scss=True, category='source_notes', type='px', value='90%'), source_notes_border_bottom_style=OptionsInfo(scss=True, category='source_notes', type='value', value='none'), source_notes_border_bottom_width=OptionsInfo(scss=True, category='source_notes', type='px', value='2px'), source_notes_border_bottom_color=OptionsInfo(scss=True, category='source_notes', type='value', value='#D3D3D3'), source_notes_border_lr_style=OptionsInfo(scss=True, category='source_notes', type='value', value='none'), source_notes_border_lr_width=OptionsInfo(scss=True, category='source_notes', type='px', value='2px'), source_notes_border_lr_color=OptionsInfo(scss=True, category='source_notes', type='value', value='#D3D3D3'), source_notes_multiline=OptionsInfo(scss=False, category='source_notes', type='boolean', value=True), source_notes_sep=OptionsInfo(scss=False, category='source_notes', type='value', value=' '), row_striping_background_color=OptionsInfo(scss=True, category='row', type='value', value='rgba(128,128,128,0.05)'), row_striping_include_stub=OptionsInfo(scss=False, category='row', type='boolean', value=False), row_striping_include_table_body=OptionsInfo(scss=False, category='row', type='boolean', value=False), container_width=OptionsInfo(scss=False, category='container', type='px', value='auto'), container_height=OptionsInfo(scss=False, category='container', type='px', value='auto'), container_padding_x=OptionsInfo(scss=False, category='container', type='px', value='0px'), container_padding_y=OptionsInfo(scss=False, category='container', type='px', value='10px'), container_overflow_x=OptionsInfo(scss=False, category='container', type='overflow', value='auto'), container_overflow_y=OptionsInfo(scss=False, category='container', type='overflow', value='auto'), quarto_disable_processing=OptionsInfo(scss=False, category='quarto', type='logical', value=False), quarto_use_bootstrap=OptionsInfo(scss=False, category='quarto', type='logical', value=False)), _has_built=False)" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pf.etable([fit_difference_in_means, fit_cuped, fit_panel])" + ] + }, + { + "cell_type": "markdown", + "id": "eb2bab0a", + "metadata": {}, + "source": [ + "The panel estimator almost exactly matches CUPED, and we do even better in terms of variance. \n", + "\n", + "However, because we are working with panel data, the error terms are likely correlated over time within each user. If we ignore this dependence, our standard errors may be underestimated, which in turn can lead to over-rejecting the null hypothesis of no effect. A common way to address this issue is to compute cluster-robust standard errors at the user level, which account for arbitrary heteroskedasticity and autocorrelation within users." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "1342b84e", + "metadata": {}, + "outputs": [], + "source": [ + "fit_panel_crv = pf.feols(\n", + " \"minutes_watched ~ treat | user + day\",\n", + " data=data,\n", + " vcov={\"CRV1\": \"user\"},\n", + " demeaner_backend=\"rust\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9ba7a9c2", + "metadata": {}, + "source": [ + "Comparing all results, the panel estimator still does quite well in terms of the size of the estimated standard errors. " + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "d99943b1", + "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", + " \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", + " minutes_watched\n", + " \n", + " minutes_post\n", + " \n", + " minutes_watched\n", + "
(1)(2)(3)(4)
coef
treat0.7615***
(0.2092)		0.1918***
(0.0288)		0.1913***
(0.0118)		0.1913***
(0.0287)
minutes_pre_c0.9992***
(0.0001)
Intercept0.3574***
(0.0148)		0.3602***
(0.0021)
fe
day--xx
user--xx
stats
Observations150000010000030000003000000
S.E. typeiidheteroheteroby: user
R20.00000.99870.99850.9985
R2 Within--0.00010.0001
Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001. Format of coefficient cell:\n", + "Coefficient \n", + " (Std. Error)
\n", + "\n", + "
\n", + " " + ], + "text/plain": [ + "GT(_tbl_data= level_0 level_1 0 \\\n", + "0 coef treat 0.7615***
(0.2092) \n", + "1 coef minutes_pre_c \n", + "2 coef Intercept 0.3574***
(0.0148) \n", + "3 fe day - \n", + "4 fe user - \n", + "5 stats Observations 1500000 \n", + "6 stats S.E. type iid \n", + "7 stats R2 0.0000 \n", + "8 stats R2 Within - \n", + "\n", + " 1 2 3 \n", + "0 0.1918***
(0.0288) 0.1913***
(0.0118) 0.1913***
(0.0287) \n", + "1 0.9992***
(0.0001) \n", + "2 0.3602***
(0.0021) \n", + "3 - x x \n", + "4 - x x \n", + "5 100000 3000000 3000000 \n", + "6 hetero hetero by: user \n", + "7 0.9987 0.9985 0.9985 \n", + "8 - 0.0001 0.0001 , _body=, _boxhead=Boxhead([ColInfo(var='level_0', type=, column_label='level_0', column_align='center', column_width=None), ColInfo(var='level_1', type=, column_label='level_1', column_align='center', column_width=None), ColInfo(var='0', type=, column_label='(1)', column_align='center', column_width=None), ColInfo(var='1', type=, column_label='(2)', column_align='center', column_width=None), ColInfo(var='2', type=, column_label='(3)', column_align='center', column_width=None), ColInfo(var='3', type=, column_label='(4)', column_align='center', column_width=None)]), _stub=, _spanners=Spanners([SpannerInfo(spanner_id='minutes_watched', spanner_level=1, spanner_label='minutes_watched', spanner_units=None, spanner_pattern=None, vars=['0', '2', '3'], built=None), SpannerInfo(spanner_id='minutes_post', spanner_level=1, spanner_label='minutes_post', spanner_units=None, spanner_pattern=None, vars=['1'], built=None)]), _heading=Heading(title=None, subtitle=None, preheader=None), _stubhead=None, _source_notes=['Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001. Format of coefficient cell:\\nCoefficient \\n (Std. Error)'], _footnotes=[], _styles=[], _locale=, _formats=[], _substitutions=[], _options=Options(table_id=OptionsInfo(scss=False, category='table', type='value', value=None), table_caption=OptionsInfo(scss=False, category='table', type='value', value=None), table_width=OptionsInfo(scss=True, category='table', type='px', value='auto'), table_layout=OptionsInfo(scss=True, category='table', type='value', value='fixed'), table_margin_left=OptionsInfo(scss=True, category='table', type='px', value='auto'), table_margin_right=OptionsInfo(scss=True, category='table', type='px', value='auto'), table_background_color=OptionsInfo(scss=True, category='table', type='value', value='#FFFFFF'), table_additional_css=OptionsInfo(scss=False, category='table', type='values', value=[]), table_font_names=OptionsInfo(scss=False, category='table', type='values', value=['-apple-system', 'BlinkMacSystemFont', 'Segoe UI', 'Roboto', 'Oxygen', 'Ubuntu', 'Cantarell', 'Helvetica Neue', 'Fira Sans', 'Droid Sans', 'Arial', 'sans-serif']), table_font_size=OptionsInfo(scss=True, category='table', type='px', value='16px'), table_font_weight=OptionsInfo(scss=True, category='table', type='value', value='normal'), table_font_style=OptionsInfo(scss=True, category='table', type='value', value='normal'), table_font_color=OptionsInfo(scss=True, category='table', type='value', value='#333333'), table_font_color_light=OptionsInfo(scss=True, category='table', type='value', value='#FFFFFF'), table_border_top_include=OptionsInfo(scss=False, category='table', type='boolean', value=True), table_border_top_style=OptionsInfo(scss=True, category='table', type='value', value='solid'), table_border_top_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_top_color=OptionsInfo(scss=True, category='table', type='value', value='#A8A8A8'), table_border_right_style=OptionsInfo(scss=True, category='table', type='value', value='none'), table_border_right_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_right_color=OptionsInfo(scss=True, category='table', type='value', value='#D3D3D3'), table_border_bottom_include=OptionsInfo(scss=False, category='table', type='boolean', value=True), table_border_bottom_style=OptionsInfo(scss=True, category='table', type='value', value='hidden'), table_border_bottom_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_bottom_color=OptionsInfo(scss=True, category='table', type='value', value='#A8A8A8'), table_border_left_style=OptionsInfo(scss=True, category='table', type='value', value='none'), table_border_left_width=OptionsInfo(scss=True, category='table', type='px', value='2px'), table_border_left_color=OptionsInfo(scss=True, category='table', type='value', value='#D3D3D3'), heading_background_color=OptionsInfo(scss=True, category='heading', type='value', value=None), heading_align=OptionsInfo(scss=True, category='heading', type='value', value='center'), heading_title_font_size=OptionsInfo(scss=True, category='heading', type='px', value='125%'), heading_title_font_weight=OptionsInfo(scss=True, category='heading', type='value', value='initial'), heading_subtitle_font_size=OptionsInfo(scss=True, category='heading', type='px', value='85%'), heading_subtitle_font_weight=OptionsInfo(scss=True, category='heading', type='value', value='initial'), heading_padding=OptionsInfo(scss=True, category='heading', type='px', value='4px'), heading_padding_horizontal=OptionsInfo(scss=True, category='heading', type='px', value='5px'), heading_border_bottom_style=OptionsInfo(scss=True, category='heading', type='value', value='solid'), heading_border_bottom_width=OptionsInfo(scss=True, category='heading', type='px', value='2px'), heading_border_bottom_color=OptionsInfo(scss=True, category='heading', type='value', value='#D3D3D3'), heading_border_lr_style=OptionsInfo(scss=True, category='heading', type='value', value='none'), heading_border_lr_width=OptionsInfo(scss=True, category='heading', type='px', value='1px'), heading_border_lr_color=OptionsInfo(scss=True, category='heading', type='value', value='#D3D3D3'), column_labels_background_color=OptionsInfo(scss=True, category='column_labels', type='value', value=None), column_labels_font_size=OptionsInfo(scss=True, category='column_labels', type='px', value='100%'), column_labels_font_weight=OptionsInfo(scss=True, category='column_labels', type='value', value='normal'), column_labels_text_transform=OptionsInfo(scss=True, category='column_labels', type='value', value='inherit'), column_labels_padding=OptionsInfo(scss=True, category='column_labels', type='px', value='4px'), column_labels_padding_horizontal=OptionsInfo(scss=True, category='column_labels', type='px', value='5px'), column_labels_vlines_style=OptionsInfo(scss=True, category='table_body', type='value', value='none'), column_labels_vlines_width=OptionsInfo(scss=True, category='table_body', type='px', value='0px'), column_labels_vlines_color=OptionsInfo(scss=True, category='table_body', type='value', value='white'), column_labels_border_top_style=OptionsInfo(scss=True, category='column_labels', type='value', value='solid'), column_labels_border_top_width=OptionsInfo(scss=True, category='column_labels', type='px', value='2px'), column_labels_border_top_color=OptionsInfo(scss=True, category='column_labels', type='value', value='black'), column_labels_border_bottom_style=OptionsInfo(scss=True, category='column_labels', type='value', value='solid'), column_labels_border_bottom_width=OptionsInfo(scss=True, category='column_labels', type='px', value='0.5px'), column_labels_border_bottom_color=OptionsInfo(scss=True, category='column_labels', type='value', value='black'), column_labels_border_lr_style=OptionsInfo(scss=True, category='column_labels', type='value', value='none'), column_labels_border_lr_width=OptionsInfo(scss=True, category='column_labels', type='px', value='1px'), column_labels_border_lr_color=OptionsInfo(scss=True, category='column_labels', type='value', value='#D3D3D3'), column_labels_hidden=OptionsInfo(scss=False, category='column_labels', type='boolean', value=False), row_group_background_color=OptionsInfo(scss=True, category='row_group', type='value', value=None), row_group_font_size=OptionsInfo(scss=True, category='row_group', type='px', value='0px'), row_group_font_weight=OptionsInfo(scss=True, category='row_group', type='value', value='initial'), row_group_text_transform=OptionsInfo(scss=True, category='row_group', type='value', value='inherit'), row_group_padding=OptionsInfo(scss=True, category='row_group', type='px', value='0px'), row_group_padding_horizontal=OptionsInfo(scss=True, category='row_group', type='px', value='5px'), row_group_border_top_style=OptionsInfo(scss=True, category='row_group', type='value', value='solid'), row_group_border_top_width=OptionsInfo(scss=True, category='row_group', type='px', value='0.5px'), row_group_border_top_color=OptionsInfo(scss=True, category='row_group', type='value', value='black'), row_group_border_right_style=OptionsInfo(scss=True, category='row_group', type='value', value='none'), row_group_border_right_width=OptionsInfo(scss=True, category='row_group', type='px', value='1px'), row_group_border_right_color=OptionsInfo(scss=True, category='row_group', type='value', value='white'), row_group_border_bottom_style=OptionsInfo(scss=True, category='row_group', type='value', value='solid'), row_group_border_bottom_width=OptionsInfo(scss=True, category='row_group', type='px', value='0.5px'), row_group_border_bottom_color=OptionsInfo(scss=True, category='row_group', type='value', value='black'), row_group_border_left_style=OptionsInfo(scss=True, category='row_group', type='value', value='none'), row_group_border_left_width=OptionsInfo(scss=True, category='row_group', type='px', value='1px'), row_group_border_left_color=OptionsInfo(scss=True, category='row_group', type='value', value='white'), row_group_as_column=OptionsInfo(scss=False, category='row_group', type='boolean', value=False), table_body_hlines_style=OptionsInfo(scss=True, category='table_body', type='value', value='none'), table_body_hlines_width=OptionsInfo(scss=True, category='table_body', type='px', value='1px'), table_body_hlines_color=OptionsInfo(scss=True, category='table_body', type='value', value='#D3D3D3'), table_body_vlines_style=OptionsInfo(scss=True, category='table_body', type='value', value='none'), table_body_vlines_width=OptionsInfo(scss=True, category='table_body', type='px', value='0px'), table_body_vlines_color=OptionsInfo(scss=True, category='table_body', type='value', value='white'), table_body_border_top_style=OptionsInfo(scss=True, category='table_body', type='value', value='solid'), table_body_border_top_width=OptionsInfo(scss=True, category='table_body', type='px', value='0.5px'), table_body_border_top_color=OptionsInfo(scss=True, category='table_body', type='value', value='black'), table_body_border_bottom_style=OptionsInfo(scss=True, category='table_body', type='value', value='solid'), table_body_border_bottom_width=OptionsInfo(scss=True, category='table_body', type='px', value='2px'), table_body_border_bottom_color=OptionsInfo(scss=True, category='table_body', type='value', value='black'), data_row_padding=OptionsInfo(scss=True, category='data_row', type='px', value='4px'), data_row_padding_horizontal=OptionsInfo(scss=True, category='data_row', type='px', value='5px'), stub_background_color=OptionsInfo(scss=True, category='stub', type='value', value=None), stub_font_size=OptionsInfo(scss=True, category='stub', type='px', value='100%'), stub_font_weight=OptionsInfo(scss=True, category='stub', type='value', value='initial'), stub_text_transform=OptionsInfo(scss=True, category='stub', type='value', value='inherit'), stub_border_style=OptionsInfo(scss=True, category='stub', type='value', value='hidden'), stub_border_width=OptionsInfo(scss=True, category='stub', type='px', value='2px'), stub_border_color=OptionsInfo(scss=True, category='stub', type='value', value='#D3D3D3'), stub_row_group_background_color=OptionsInfo(scss=True, category='stub', type='value', value=None), stub_row_group_font_size=OptionsInfo(scss=True, category='stub', type='px', value='100%'), stub_row_group_font_weight=OptionsInfo(scss=True, category='stub', type='value', value='initial'), stub_row_group_text_transform=OptionsInfo(scss=True, category='stub', type='value', value='inherit'), stub_row_group_border_style=OptionsInfo(scss=True, category='stub', type='value', value='solid'), stub_row_group_border_width=OptionsInfo(scss=True, category='stub', type='px', value='2px'), stub_row_group_border_color=OptionsInfo(scss=True, category='stub', type='value', value='#D3D3D3'), source_notes_padding=OptionsInfo(scss=True, category='source_notes', type='px', value='4px'), source_notes_padding_horizontal=OptionsInfo(scss=True, category='source_notes', type='px', value='5px'), source_notes_background_color=OptionsInfo(scss=True, category='source_notes', type='value', value=None), source_notes_font_size=OptionsInfo(scss=True, category='source_notes', type='px', value='90%'), source_notes_border_bottom_style=OptionsInfo(scss=True, category='source_notes', type='value', value='none'), source_notes_border_bottom_width=OptionsInfo(scss=True, category='source_notes', type='px', value='2px'), source_notes_border_bottom_color=OptionsInfo(scss=True, category='source_notes', type='value', value='#D3D3D3'), source_notes_border_lr_style=OptionsInfo(scss=True, category='source_notes', type='value', value='none'), source_notes_border_lr_width=OptionsInfo(scss=True, category='source_notes', type='px', value='2px'), source_notes_border_lr_color=OptionsInfo(scss=True, category='source_notes', type='value', value='#D3D3D3'), source_notes_multiline=OptionsInfo(scss=False, category='source_notes', type='boolean', value=True), source_notes_sep=OptionsInfo(scss=False, category='source_notes', type='value', value=' '), row_striping_background_color=OptionsInfo(scss=True, category='row', type='value', value='rgba(128,128,128,0.05)'), row_striping_include_stub=OptionsInfo(scss=False, category='row', type='boolean', value=False), row_striping_include_table_body=OptionsInfo(scss=False, category='row', type='boolean', value=False), container_width=OptionsInfo(scss=False, category='container', type='px', value='auto'), container_height=OptionsInfo(scss=False, category='container', type='px', value='auto'), container_padding_x=OptionsInfo(scss=False, category='container', type='px', value='0px'), container_padding_y=OptionsInfo(scss=False, category='container', type='px', value='10px'), container_overflow_x=OptionsInfo(scss=False, category='container', type='overflow', value='auto'), container_overflow_y=OptionsInfo(scss=False, category='container', type='overflow', value='auto'), quarto_disable_processing=OptionsInfo(scss=False, category='quarto', type='logical', value=False), quarto_use_bootstrap=OptionsInfo(scss=False, category='quarto', type='logical', value=False)), _has_built=False)" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pf.etable(\n", + " [fit_difference_in_means, fit_cuped, fit_panel, fit_panel_crv],\n", + " digits=4,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "d78bd300", + "metadata": {}, + "source": [ + "The panel errors produces point estimates and standard errors for the treatment variable that are very close to CUPED. " + ] + }, + { + "cell_type": "markdown", + "id": "9f4e3bc4", + "metadata": {}, + "source": [ + "The sceptical reader might now object that this is all far from overwhelming evidence. After all , we fit each model only once! So it might well be \n", + "that we were just lucky and our results are driven by sampling error. In the next section, we will argue that this is not the case by means of a monte carlo simulation. " + ] + }, + { + "cell_type": "markdown", + "id": "bc9cf093", + "metadata": {}, + "source": [ + "### Digression: Heterogeneous Effects over Time via Panel Estimators" + ] + }, + { + "cell_type": "markdown", + "id": "bcf15cb8", + "metadata": {}, + "source": [ + "One example of the panel estimator is that it allows us to track treatment effects over time. This allows us to track novelty effects where in the beginning, \n", + "users interact a lot with a new feature, i.e. they spend a lot of time playing a gamee, but over time lose interest. " + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "97d0a67c", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAJOCAYAAACqS2TfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbQBJREFUeJzt3QuczOX+wPHvYrHuSlgs61IuudUqYV3qLKETkvincqlTUkKi6GJ3K5Uu2KJsFLo4kZyubEqupZSii1tJYa1Lp1jW5rLm//o+nZlmdmeZ2Z3fXD/v12uY3+/3zDzPzM488/v+nluUzWazCQAAAAAA8LkSvn9KAAAAAACgCLoBAAAAALAIQTcAAAAAABYh6AYAAAAAwCIE3QAAAAAAWISgGwAAAAAAixB0AwAAAABgEYJuAAAAAAAsQtANAAAAAIBFCLoBAIClVq5cKVFRUeZ/f5o7d67J95dffvFrvgAAOCPoBgAgX5BW2O3zzz8PdBHls88+k5SUFDl06JDHj3nvvfekc+fOUr16dSlXrpw0aNBA+vfvLxkZGY40e/fuNc+7ceNGCXZaTue/i76mZs2ayYMPPijZ2dk+yWP+/Pkybdo0nzwXACCylQp0AQAACDYPP/yw1K9fv8D+Ro0aSTAE3ampqTJkyBCpUqXKWdM//fTTMm7cOBN0T5gwwQSoP/30k3z88cfyxhtvSPfu3R1Btz5vfHy8tG7dWkLBCy+8IBUqVJCjR4/KsmXLZNKkSfLJJ5/Ip59+aoLx4gbd33//vYwePdpn5QUARCaCbgAA8unRo4e0adNGQt2pU6fkkUceka5du5qgNL8DBw5IKOvXr59Uq1bN3L/99tvl2muvlcWLF5seCe3atQt08QAAMOheDgCAF06ePCnnnHOODB06tMAx7dpctmxZGTt2rGPf8ePHJTk52bSSlylTRuLi4uTee+81+51py+yIESPk7bfflubNm5u0F154oUsXcO1Wra3WSlvi7d2rCxuz/Ntvv5kydejQwe1x7W6udKz1JZdcYu7r67I/r3a3V9r6rS3r+XXp0sXcnO3Zs0f69Okj5cuXN89/9913F3it+n5ER0fLwYMHCzznbbfdZlrw//zzT/HWFVdcYf7fuXPnGdM9//zz5r3V97hWrVpy5513unTX19f0wQcfyK+//up4L/Q9AACgKGjpBgAgn8OHD5uA1ZkGXueee64JFq+55hrTopqeni6lS5d2pNGAWQPM//u//zPbp0+fll69esnatWtNMNm0aVP57rvvZOrUqbJ9+3aT3pmm0+e94447pGLFivLss8+a1ttdu3aZvPv27Wse9+9//9s8h72V97zzznP7OjTojYmJMWO677rrLnOxwB0tl3apnzhxoilnx44dzf727dt79b7l5ubKP/7xD1PekSNHmoD21VdfNV2+nd10000mvwULFpgLDXYnTpyQRYsWmdesFy+8tWPHDvO/vleF0QsX2o0+KSlJhg8fLtu2bTPd1L/88kvTLV3/vg888ID5DOgFBH2flXZjBwCgSGwAAMCYM2eOTX8a3d3KlCnjSPfhhx+afe+9957L43v27Glr0KCBY/vVV1+1lShRwrZmzRqXdDNnzjSP//TTTx37dLt06dK2n376ybFv06ZNZv9zzz3n2PfUU0+ZfTt37vToNU2cONGkL1++vK1Hjx62SZMm2TZs2FAg3ZdffmnS6XuQX7169WyDBw8usL9z587mZjdt2jTzHAsXLnTsy8nJsTVq1MjsX7FihWN/u3btbG3btnV5vsWLFxdI505ycrJJt23bNtvBgwfNe5Genm7+RjVq1DB5Ov897e/VgQMHzHvcrVs3W15enuP5pk+fbtK9/PLLjn1XXXWVed0AABQX3csBAMhnxowZ8tFHH7ncli5d6tKNWVuZtaXW7o8//jDpBgwY4Nj35ptvmlbkJk2amJZz+83eDXrFihUu+Wrra8OGDR3bLVu2lEqVKsnPP/9c5Neirbo6KdhFF10kH374oWnFTUhIkIsvvli2bNkivrRkyRKJjY01Y63tdOI2bT3Pb9CgQfLFF184WqfV66+/brrf66RvnmjcuLFp5deu9sOGDTNd+LVbuObpjk4ep63pOjlaiRJ/nwLdeuut5n3WxwIA4Gt0LwcAIJ9LL730jBOplSpVynSB1mBWu5Pr2GDtFq7jvZ2D7h9//NEEtoV1/84/kVndunULpKlataoJ6Ivj+uuvNzcd362Bro7V1rJfffXVZobuonTldkfHQGvgm3/mcA2O89P3SYNfDbS1W7t2537//ffNGHBPZx5/6623TLCsXcLr1KnjcsGisPK5K48OEdBl1OzHAQDwJYJuAACKQMdt65hubQHXicMWLlxoWrRbtWrlSKNjulu0aCFTpkxx+xzaquusZMmSbtP91fu8+DRA1ZnM9aaB6rx580wQfraW5cKC4Ly8vELLfDZ6MeGf//ynI+jWsdx6AePGG2/0+Dk6derkGNcOAECwIugGAKAINODTrtTaxTwxMdFMFqZdt51py+umTZvM5GLFXTfazlfPoy35GnRnZWWd9Xk1QHae3dtOW4a1hdiuXr16puVcLxI4P59OVuaOdjHv3bu3mcRMg2/tAq+ziltFy2cvj3O5tcu5zniu3ft9/T4DAMCYbgAAikDHBOvYZZ0ZXGfo1jWxnbuWq/79+0tmZqbMmjXL7UzfOTk5XuerS3Epd0FwfseOHZN169a5PWYfo27van2m59WLB7r2tQandtoVfPfu3S7pevbsKXv37jWt1s5lePHFFwtdD11bqidPniyrVq3yqpW7KDSo1q7kOiu8c++Bl156yXRvv+qqqxz79P3QfQAAFBct3QAAuAlIt27dWmC/LqHl3EKqQfZzzz1n1p3WbuQ6aVr+pbG02/ntt99uJk3T9bK1S7Y+t+7Xic3ONHbcHZ0ETWmrunZx127iOjbbHjQ704BXy3zZZZdJ9+7dTXd2Dap1qbI1a9aYbvHaumwPrHV97JkzZ5rlyvT52rZtayYp+9e//mUCaX0OvZCgk5+99tprBcZQ64Rk06dPNy3YGzZsMD0B9IJEYRObadn1NehjtJu6jju3ko6tnzBhgplcTl+LLuemrd66breuU+4c9Ov7rL0YxowZY47pkmH6PgMA4LViz38OAEAELBnmbjmt06dP2+Li4syxRx991O1znjhxwjZ58mTbhRdeaJa0qlq1qi0hIcGWmppqO3z4sCOdPsedd97p0XJdjzzyiK127dpmObIzLR928uRJ26xZs2x9+vQxz6P5lytXznbRRReZpceOHz/ukv6dd96xNWvWzFaqVKkCr/eZZ54xeepzdOjQwfbVV18VWDJM/frrr7ZevXqZfKpVq2YbNWqULSMjo9ClwNavX2+O6TJenrIvGabLhZ1J/iXDnJcIa9KkiS06OtosMTZ8+HDbH3/84ZLm6NGjtoEDB9qqVKlinoPlwwAARRWl/3gfqgMAABSfjnlv3bq1vPLKK6ZnAAAA4YYx3QAAIGB0vLt23e7bt2+giwIAgCUY0w0AAPxOJ6DbvHmzmWRtxIgRbsekAwAQDuheDgAA/C4+Pl72798vV155pZlsTSdvAwAgHBF0AwAAAABgEcZ0AwAAAABgEYJuAAAAAAAsEnETqZ0+fVr27t1rxo5FRUUFujgAAAAAgBCkI7WPHDkitWrVkhIlCm/PjrigWwPuuLi4QBcDAAAAABAGdu/eLXXq1Cn0eMQF3fbZUfWNqVSpUqCLAwAAAAAIQdnZ2aZB92wrcERc0G3vUq4BN0E3AAAAAKA4zjZsmYnUAAAAAACwCEE3AAAAAAAWIegGAAAAAMAiBN0AAAAAAFiEoBsAAAAAAIsQdAMAAAAAYBGCbgAAAAAALELQDQAAAACARQi6AQAAAACwCEE3AAAAAAAWIegGAAAAAMAiBN0AAAAAAFiEoBsAAAAAAIsQdAMAAAAAYBGCbgAAAAAALELQDQAAAACARQi6AQAAAACwCEE3AAAAAAAWIegGAAAAAMAipax6YgAAACCcZWUdkayso14/Lja2gsTGVrSkTACCD0E3AAAAUATp6RskNXWV149LTu4sKSldLCkTgOBD0A0AAAAUwbBhCdKrV2OXfbm5JyUxcY65v3btUImJiXbb0g0gchB0AwAAAEWgXcTzdxPPyTnhuN+6dU0pX750AEoGIJgwkRoAAAAAABYh6AYAAAAAwCIE3QAAAAAAWISgGwAAAAAAixB0AwAAAABgEYJuAAAAAAAsQtANAAAAAIBFCLoBAAAAALAIQTcAAAAAABYh6AYAAAAAwCIE3QAAAAAAWISgGwAAAAAAixB0AwAAAABgEYJuAAAAAAAsQtANAAAAAIBFCLoBAAAAALAIQTcAAAAAABYh6AYAAAAAwCIE3QAAAAAAWKSUVU8MAAAA+EtW1hHJyjrq9eNiYytIbGxFS8oEAIqgGwAAACEvPX2DpKau8vpxycmdJSWliyVlAgBF0A0AAICQN2xYgvTq1dhlX27uSUlMnGPur107VGJiot22dAOAlQi6AQAAEPK0i3j+buI5OScc91u3rinly5cOQMkARDomUgMAAAAAwCIE3QAAAAAAWISgGwAAAAAAixB0AwAAAABgEYJuAAAAAAAsQtANAAAAAIBFCLoBAAAAALAI63QDAADAZ7KyjkhW1lGvHxcbW6HAOtsAEA4IugEAAOAz6ekbJDV1ldePS07uLCkpXSwpEwAEEkE3AAAAfGbYsATp1auxy77c3JOSmDjH3F+7dqjExES7bekGgHBE0A0AAACf0S7i+buJ5+SccNxv3bqmlC9fOgAlA4AInUhtxowZEh8fL2XLlpW2bdvK+vXrz5j+0KFDcuedd0psbKyUKVNGLrjgAlmyZInfygsAAAAAQEi0dC9YsEDGjBkjM2fONAH3tGnT5Morr5Rt27ZJ9erVC6Q/ceKEdO3a1RxbtGiR1K5dW3799VepUqVKQMoPAAAAAEDQBt1TpkyRW2+9VYYOHWq2Nfj+4IMP5OWXX5bx48cXSK/7f//9d/nss88kOvqvsUDaSg4AAAAAQDAKWPdybbXesGGDJCUl/V2YEiXM9rp169w+5t1335V27dqZ7uU1atSQ5s2by2OPPSZ5eXl+LDkAAAAAAEHe0v3bb7+ZYFmDZ2e6vXXrVreP+fnnn+WTTz6RG264wYzj/umnn+SOO+6QkydPSnJystvHHD9+3NzssrOzffxKAAAAAAAI0onUvHH69GkznvvFF1+UhIQEGTBggDzwwAOmW3phHn/8calcubLjFhcX59cyAwAAIHLk5Z123F+9+leXbQCRKWBBd7Vq1aRkyZKyf/9+l/26XbNmTbeP0RnLdbZyfZxd06ZNZd++faa7ujsTJkyQw4cPO267d+/28SsBAAAARBYv3iLNmj3v2O7Zc77Ex6eZ/QAiV8CC7tKlS5vW6uXLl7u0ZOu2jtt2p0OHDqZLuaaz2759uwnG9fnc0WXFKlWq5HIDAAAAfEkD6379Fkpm5hGX/ZmZ2WY/gTcQuQLavVyXC5s1a5bMmzdPtmzZIsOHD5ecnBzHbOaDBg0yLdV2elxnLx81apQJtnWmc51ITSdWAwAAAAJBu5CPGpUhNlvBY/Z9o0dn0NUciFABXTJMx2QfPHhQJk6caLqIt27dWjIyMhyTq+3atcvMaG6n47E//PBDufvuu6Vly5ZmnW4NwO+7774AvgoAAABEsjVrdsmePYVP1quB9+7d2SZdly4sdwtEmoAG3WrEiBHm5s7KlSsL7NOu559//rkfSgYAAACcXVbWEZ+mAxBeQmr2cgAAACDYxMZW9Gk6AOGFoBsAAAAoho4d60qdOpUkKsr9cd0fF1fJpAMQeQi6AQAAgGIoWbKEpKV1N/fzB9727WnTupt0ACIP33wAAACgmPr2bSqLFvWXWrVcu5BrC7ju1+MAIlPAJ1IDAAAAwoEG1klJ9aVy5clme8mSgdKtW0NauIEIRw0AAAAA+IhzgN2pUz0CbgAE3QAAAAAAWIWgGwAAAAAAixB0AwAAAABgESZSAwAACENZWUckK+uo14+Lja0gsbGuM3ADAIqOoBsAACAMpadvkNTUVV4/Ljm5s6SkdLGkTAAQiQi6AQAAwtCwYQnSq1djl325uSclMXGOub927VCJiYl229INAPAdgm4AAIAwpF3E83cTz8k54bjfunVNKV++dABKBgCRhYnUAAAAAACwCEE3AAAAAAAWIegGAAAAAMAiBN0AAAAAAFiEoBsAAAAAAIsQdAMAACAs5eWddtxfvfpXl20A8BeCbgAAAISdxYu3SLNmzzu2e/acL/HxaWY/APgTQTcAAADCigbW/fotlMzMIy77MzOzzX4CbwD+RNANAACAsKFdyEeNyhCbreAx+77RozPoag7Abwi6AQAAEDbWrNkle/ZkF3pcA+/du7NNOgDwB4JuAAAAhI2srCM+TQcAxUXQDQAAgLARG1vRp+kAoLgIugEAABA2OnasK3XqVJKoKPfHdX9cXCWTDgD8gaAbAAAAYaNkyRKSltbd3M8feNu3p03rbtIBgD9Q2wAAACCs9O3bVBYt6i+1arl2IdcWcN2vxwHAXwi6AQAAYCnn5blWr/7VL8t1aWC9efMdju0lSwbKzp2jCLgB+B1BNwAAACyzePEWadbsecd2z57zJT4+zey3mnMX8k6d6tGlHEBAUPMAAADAEhpY9+u3UDIzXZfnyszMNvv9EXgDQKARdAMAAMDntAv5qFEZYrMVPGbfN3p0hl+6mgNAIBF0AwAAwOfWrNkle/ZkF3pcA+/du7NNOgAIZwTdAAAA8LmsrCM+TQcAoYqgGwAAAD4XG1vRp+kAIFQRdAMAAMDnOnasa9bFjopyf1z3x8VVMukAIJwRdAMAAMDndHmutLTu5n7+wNu+PW1ad5bxAhD2qOUAAABgib59m8qiRf2lVi3XLuTaAq779TgAhLtSgS4AAAAAwpcG1klJ9aVy5clme8mSgdKtW0NauAFEDGo7AAAAWMo5wO7UqR4BN4CIQo0HAAAAAIBFCLoBAAAAALAIQTcAAAAAAIGcSC07O9vjJ6xUqVJxygMAAAAAQGQF3VWqVJGo/Ass5mOz2UyavLw8X5UNAAAACFpZWUckK+uoy77c3JOO+xs37pOYmOgCj4uNrSCxsa7LqAGI8KB7xYoV1pcEAAAACCHp6RskNXVVoccTE+e43Z+c3FlSUrpYWDIAwcSjoLtz587WlwQAAAAIIcOGJUivXo29fpy2dAOIHB4F3e4cO3ZMdu3aJSdOnHDZ37JlS1+UCwAAAAhq2kWcbuIAfB50Hzx4UIYOHSpLly51e5wx3QAAAAAAFHHJsNGjR8uhQ4fkiy++kJiYGMnIyJB58+bJ+eefL++++663TwcAAAAAQNjyuqX7k08+kXfeeUfatGkjJUqUkHr16knXrl3NUmGPP/64XHXVVdaUFAAAAMWSl3facX/16l+lW7eGUrKk120wAAAveF3L5uTkSPXq1c39qlWrmu7mqkWLFvL11197+3QAAADwg8WLt0izZs87tnv2nC/x8WlmPwAgiFq6GzduLNu2bZP4+Hhp1aqVpKenm/szZ86U2NhYa0oJAAAQRms5e8KXazlrYN2v30Kx2Vz3Z2Zmm/2LFvWXvn2b+iQvAEAxg+5Ro0ZJVlaWuZ+cnCzdu3eX119/XUqXLi1z58719ukAAAAiei3nwvhqLWftUj5qVEaBgFvpvqgonbMnQ3r3bkxXcwAIhqD7xhtvdNxPSEiQX3/9VbZu3Sp169aVatWq+bp8AAAAYbeWc27uSUlMnGPur107VGJioi1by3nNml2yZ092occ18N69O9uk69Il3id5AgCKEXSvXbtWEhMTHdvlypWTiy++2NunAQAAiNi1nHNyTjjut25dU8qXL21p93ZfpgMAeMfrPkRXXHGF1K9fX+6//37ZvHmztw8HAACAH3k6LtxX48cBAMUMuvfu3Sv33HOPrFq1Spo3by6tW7eWp556Svbs2ePtUwEAAMBiHTvWlTp1Kpmx2+7o/ri4SiYdACAIgm4dtz1ixAj59NNPZceOHXLdddfJvHnzzAzm2goOAACA4KGTo6WldTf38wfe9u1p07oziRoAWKRYtat2Mx8/frw88cQTZp1ubf0GAABAcNHlwHRZsFq1XLuQaws4y4UBQJAG3drSfccdd5i1uQcOHGi6mn/wwQe+LR0AAAB8QgPrzZvvcGwvWTJQdu4cRcANAME2e7m2bC9YsMCM7e7ataukpaVJ7969zSzmAAAACF7OXcg7dapHl3IACMage82aNTJu3Djp378/63IDAAAAAHAGXl3ePHnypDRu3Fh69OhBwA0AAAAAgC+D7ujoaHnrrbe8eQgAAAAAABHL64E8ffr0kbffftua0gAAAAAAEMljus8//3x5+OGHzezlCQkJUr58eZfjI0eO9GX5AAAAAACInKD7pZdekipVqsiGDRvMzVlUVBRBNwAAAAAARQ26d+7c6e1DAAAAAACISF4H3XYnTpwwAXjDhg2lVKkiPw0AAABQbFlZRyQr66jLvtzck477Gzfuk5iY6AKPi42tILGxFf1SRgCRyeto+dixY3LXXXfJvHnzzPb27dulQYMGZl/t2rVl/PjxVpQTAAAAKFR6+gZJTV1V6PHExDlu9ycnd5aUlC4WlgxApPM66J4wYYJs2rRJVq5cKd27d3fsT0pKkpSUFIJuAAAA+N2wYQnSq1djrx+nLd0AEFRBty4XtmDBArnsssvMxGl2F154oezYscPX5QMAAADOSruI000cQFis033w4EGpXr16gf05OTkuQTgAAAAAAJHO66C7TZs28sEHHzi27YH27NmzpV27dr4tHQAAAAAAkdS9/LHHHpMePXrI5s2b5dSpU5KWlmbuf/bZZ7JqVeGTVwAAAAAAEGm8bulOTEyUjRs3moC7RYsWsmzZMtPdfN26dZKQkGBNKQEAAAAACEFFWmBb1+aeNWuW70sDAAAAAEAkt3SXLFlSDhw4UGD/f//7X3MMAAAAAAAUsaXbZrO53X/8+HEpXbq0t08HAACAMJKVdUSyso667MvNPem4v3HjPomJiXa7XjZLfgGI6KD72WefdcxWrjOVV6hQwXEsLy9PVq9eLU2aNLGmlAAAAAgJ6ekbJDW18Ml1ExPnuN2fnNxZUlK6WFgyAAjyoHvq1KmOlu6ZM2e6dCXXFu74+HizHwAAAJFr2LAE6dWrsdeP05ZuAIjooHvnzp3m/8svv1wWL14sVatWtbJcAAAACEHaRZxu4gBQjDHdK1as8PYhAAAAAABEJK+Dbh2/PXfuXFm+fLmZxfz06dMuxz/55BNflg8AAMCyCb48wQRfAAC/Bt2jRo0yQfdVV10lzZs3NxOrAQAAhPoEX4Vhgi8AgF+D7jfeeEMWLlwoPXv2LFbGAAAAgZ7gS5eyss+mvXbt0EKXsgIAwG9Bt85U3qhRoyJnCAAAECwTfOXknHDcb926ppQvXzoAJQMAhLMS3j7gnnvukbS0NLN0GAAAAAAA8GFL99q1a80M5kuXLpULL7xQoqNdu2HpcmIAAAAAAKAIQXeVKlXkmmuusaY0AAAAESAv7+/VX1av/lW6dWsoJUt63QERABCOQfecOX9NNgIAAADvLV68RUaOXOrY7tlzvtSpU0nS0rpL375NA1o2AIDvcUkVAADAjwF3v34LJTPziMv+zMxss1+PAwAitKX7oosu8mhN7q+//rq4ZQIAAAjLLuWjRmWIu7lodZ+eZo0enSG9ezemqzkARGLQ3adPH2tLAgAAEMbWrNkle/ZkF3pcA+/du7NNui5d4v1aNgBAEATdycnJFhYDAAAgvGVlHfFpOgBAaKDvEgAAgB/Exlb0aToAQGgg6AYAAPCDjh3rmlnKC5siR/fHxVUy6QAAEbxkGAAAALynk6PpsmA6S7kG2M4TqtkD8WnTuvtsEjXtpp6VddRlX27uScf9jRv3SUxMdIHHxcZWoLUdAMKtpXvGjBkSHx8vZcuWlbZt28r69es9etwbb7xhZlRnkjcAABAKdB3uRYv6S61arkGttoDrfl+u052evkESEl50uSUmznEc1/v5j+tNHwcACJKW7j///NMEysWxYMECGTNmjMycOdME3NOmTZMrr7xStm3bJtWrVy/0cb/88ouMHTtWOnbsWKz8AQAA/EkD66Sk+lK58mSzvWTJQOnWraHPlwkbNixBevVq7PXjtKUbABDAoPv06dMyadIkEyTv379ftm/fLg0aNJCHHnrItFbfcsstXj3flClT5NZbb5WhQ4eabX3eDz74QF5++WUZP36828fk5eXJDTfcIKmpqbJmzRo5dOiQty8DAAAgYJwD7E6d6lmyLrd2EaebOAAEntc1/KOPPipz586VJ598UkqXLu3Y37x5c5k9e7ZXz3XixAnZsGGDJCUl/V2gEiXM9rp16wp93MMPP2xawb0N8AEAAAAACOqg+5VXXpEXX3zRtDSXLFnSsb9Vq1aydetWr57rt99+M63WNWrUcNmv2/v27XP7mLVr18pLL70ks2bN8iiP48ePS3Z2tssNAAAAAICgDLozMzOlUaNGbrudnzz594yYVjhy5IjcdNNNJuCuVq2aR495/PHHpXLlyo5bXFycpWUEAAAAAKDIY7qbNWtmxlHXq1fPZf+iRYvkoosu8uq5NHDW1nIdG+5Mt2vWrFkg/Y4dO8wEaldffbVLsK9KlSplJl9r2LChy2MmTJhgJmqz05ZuAm8AAAAAQFAG3RMnTpTBgwebFm8NeBcvXmyCXe12/v7773v1XDomPCEhQZYvX+5Y9kufU7dHjBhRIH2TJk3ku+++c9n34IMPmhbwtLQ0t8F0mTJlzA0AAAAAgKAPunv37i3vvfeemcysfPnyJgi/+OKLzb6uXbt6XQBthdYgvk2bNnLppZeaJcNycnIcs5kPGjRIateubbqJ6/JkOmGbsypVqpj/8+8HAAAAACAk1+nWtbE/+ugjnxRgwIABcvDgQRO86+RprVu3loyMDMfkart27TIzmgMAAAAAEGqibDabzZsHfPnll6YLeNu2bV32f/HFF2Z8trZYBzMd060Tqh0+fFgqVaoU6OIAAIAAysk5IRUqPG7uHz06QcqXLx3W+QIA/B9bet2EfOedd8ru3bsL7Ncx3noMAAAAAAAUMejevHmzGcOdn85crscAAAAAAEARg26dCTz/El8qKyvLLNsFAAAAAACKGHR369bNrH2t/dbtDh06JPfff3+RZi8HAAAAACBced00/fTTT0unTp2kXr16pku52rhxo5lt/NVXX7WijAAAAAAAREbQrWtmf/vtt/L666/Lpk2bJCYmxqypff3110t0dLQ1pQQAAAAAIAQVaRB2+fLl5bbbbvN9aQAAAAAAiPSg+8cff5QVK1bIgQMHzJrdziZOnOirsgEAAAAAEFlB96xZs2T48OFSrVo1qVmzpkRFRTmO6X2CbgAAECry8v5uPFi9+lfp1q2hlCzp9TyzAAD4Luh+9NFHZdKkSXLfffd5+1AAAICgsXjxFhk5cqlju2fP+VKnTiVJS+suffs2DWjZAADhw+tLuX/88Ydcd9111pQGAADATwF3v34LJTPziMv+zMxss1+PAwAQkKBbA+5ly5b5JHMAAIBAdCkfNSpDbLaCx+z7Ro/OcOl6DgCA37qXN2rUSB566CH5/PPPpUWLFgWWCRs5cmSRCwMAAGC1NWt2yZ492YUe18B79+5sk65Ll3i/lg0AEH68DrpffPFFqVChgqxatcrcnOlEagTdAAAgmGVlHfFpOgAAfBp079y509uHAAAAuASzWVlHvX5cbGwFiY2tWOz8PX0OX+QFAECR1ulWJ06cMAF4w4YNpVSpIj8NAACIMOnpGyQ11bW3nCeSkztLSkqXYuffsWNdM0u5Tprmbly3roaqxzUdAADF5XW0fOzYMbnrrrtk3rx5Znv79u3SoEEDs6927doyfvz4YhcKAACEr2HDEqRXr8Yu+3JzT0pi4hxzf+3aoRIT4zpnjL2l2xd0HW5dFkxnKdcA2znw1m01bVp31usGAAQm6J4wYYJs2rRJVq5cKd27d3fsT0pKkpSUFIJuAABw1m7b+btu5+SccNxv3bqmlC9f2tIy6Drcixb1N+t0Oy8bpi3cGnCzTjcAIGBB99tvvy0LFiyQyy67zEycZnfhhRfKjh07fFYwAAAAK2lgnZRUXypXnmy2lywZKN26NaSFGwDgU17/qhw8eFCqV69eYH9OTo5LEA4AABDsnAPsTp3qEXADAHzO61+WNm3ayAcffODYtgfas2fPlnbt2vm2dAAAAAAARFL38scee0x69OghmzdvllOnTklaWpq5/9lnnxVYtxsAAAAAgEjmdUt3YmKibNy40QTcLVq0kGXLlpnu5uvWrZOEhARrSgkAAAAAQAgq0gLbujb3rFmzfF8aAAAAAAAiLejOzs72+AkrVapUnPIAAAAAABBZQXeVKlU8npk8Ly+vuGUCAAAAACBygu4VK1Y47v/yyy8yfvx4GTJkiGO2ch3PPW/ePHn88cetKykAAEAIyso6IllZR1325eaedNzfuHGfxMREF3hcbGwFiY2t6JcyAgACHHR37tzZcf/hhx+WKVOmyPXXX+/Y16tXLzOp2osvviiDBw+2pqQAAAAhKD19g6SmFr7CS2LiHLf7k5M7S0pKFwtLBgAIyonUtFV75syZbtfv/te//uWrcgEAAISFYcMSpFevxl4/Tlu6AQARGHTHxcWZmcuffPJJl/2zZ882xwAAAPA37SJON3EAiFxeB91Tp06Va6+9VpYuXSpt27Y1+9avXy8//vijvPXWW1aUEQAAAACAkFTC2wf07NlTtm/fLldffbX8/vvv5qb3dZ8eAwAAAAAARWzpVtqN/LHHHivKQwEAAAAAiBgeBd3ffvutNG/eXEqUKGHun0nLli19VTYAAAAAAMI/6G7durXs27dPqlevbu5HRUWJzWYrkE735+XlWVFOAAAAAADCM+jeuXOnnHfeeY77AAAAAADAR0H3NddcI8uXL5eqVavKvHnzZOzYsVKuXDlPHgoAAAAAQMTyaPbyLVu2SE5OjrmfmpoqR48etbpcAAAAAABEzpjuoUOHSmJiohnL/fTTT0uFChXcpp04caKvywgAAAAAQPgG3XPnzpXk5GR5//33zWRpS5culVKlCj5UjxF0AwAAAADgRdDduHFjeeONN8x9XTZMx3frTOYAAAAAAKCYQbez06dPe/sQAAAAAAAiktdBt/rxxx9lxYoVcuDAgQJBON3LAQAAAAAoYtA9a9YsGT58uFSrVk1q1qxpxnHbMaYbAAAURV7e3xfxV6/+Vbp1ayglS3q0yAoAAOEVdD/66KMyadIkue+++6wpEQAAiCiLF2+RkSOXOrZ79pwvdepUkrS07tK3b9OAlg0AgOLy+hLyH3/8Idddd12xMwYAANCAu1+/hZKZecRlf2ZmttmvxwEAiKigWwPuZcuWWVMaAAAQUV3KR43KEJut4DH7vtGjM1y6ngMAEPbdyxs1aiQPPfSQfP7559KiRQuJjo52OT5y5Ehflg8AAISpNWt2yZ492YUe18B79+5sk65Ll3i/lg0AgIAF3S+++KJUqFBBVq1aZW7OdCI1gm4AAOCJrKwjPk0HAEBYBN07d+60piQAACCixMZW9Gk6AADCZp1uO9v/Blw5LxsGAADgiY4d65pZynXSNHfjuvX0Qo9rOl/QFvOsrKMu+3JzTzrub9y4T2JiXIfNqdjYCgT+AAD/Bt2vvPKKPPXUU/Ljjz+a7QsuuEDGjRsnN910U9FLAgAAIoquw63Lguks5RpgOwfe9uv506Z199l63enpGyQ11XVonLPExDlu9ycnd5aUlC4+KQMAIPJ4HXRPmTLFTKQ2YsQI6dChg9m3du1auf322+W3336Tu+++24pyAgCAMKTrcC9a1N+s0+28bJi2cGvA7ct1uocNS5BevRp7/Tht6QYAoKiibPY+4h6qX7++pKamyqBBg1z2z5s3T1JSUoJ+zHd2drZUrlxZDh8+LJUqVQp0cQAAgPl9/lMqV55s7i9ZMlC6dWvosxZuAAACGVt6/WuWlZUl7du3L7Bf9+kxAAAAbzkH2J061SPgBgCEjRJFWad74cKFBfYvWLBAzj//fF+VCwAAAACAyBvTrV3LBwwYIKtXr3aM6f70009l+fLlboNxAAAAAAAildct3ddee6188cUXUq1aNXn77bfNTe+vX79errnmGmtKCQAAAABApCwZlpCQIK+99prvSwMAAAAAQCS2dO/du1fGjh1rZmjLT2dr03W69+/f7+vyAQAAAAAQ/kG3rs+tAbe7qdB1mvQjR46YNAAAAAAAwMugOyMjo8Da3M702Pvvv+/p0wEAAAAAEPY8Drp37twpdevWLfR4nTp15JdffvFVuQAAAAAAiJygOyYm5oxBtR7TNAAAAAAAwMugu23btvLqq68WevyVV16RSy+91NOnAwAAAAAg7Hm8ZJjOXN61a1czaZrOVF6jRg2zX2csf/LJJ2Xu3LmybNkyK8sKAAAAAEB4Bt2XX365zJgxQ0aNGiVTp041s5hHRUWZ5cKio6PlueeekyuuuMLa0gIAAAAAEI5Btxo2bJj885//lIULF8pPP/0kNptNLrjgAunXr5+ZSA0AAAAAABQx6Fa1a9eWu+++29uHAQAAAAAQcTyeSA0AAAAAAHiHoBsAAAAAAIsQdAMAAAAAYBGCbgAAAAAAgmUiNWfff/+9rFq1SvLy8qRDhw6SkJDgu5IBAAAAABCpLd26Zvc//vEPE3SvWLHCrNE9adIk35YOAAAAAIBIaOnevXu3xMXFObanT58uP/zwg1SrVs1sr1u3Tnr16iUPPPCANSUFAAAAACBcW7qTkpIkLS1NbDab2T733HMlIyNDjh8/LkeOHJGPP/5YzjvvPCvLCgAAAABAeAbdX375pWzbtk3atm0rGzdulBdffFGmTp0qMTExUqVKFVmwYIHMmzfP2tICAAAAABCO3csrVaokzz//vHz22WcyZMgQM4Z7zZo1ZhI1vWngDQAAAAAAijGRWvv27eWrr76SqlWrykUXXSSrV68m4AYAAAAAoDgt3adOnTJdyrds2SKtWrWS+++/XwYMGCC33367zJ0710ysVqNGDU+fDgAARKisrCOSlXXUZV9u7knH/Y0b90lMTHSBx8XGVpDY2Ip+KSMAAH4Pum+55RYzrltnKJ8zZ458++238uyzz8onn3wiL730krRr107GjRsnw4cP91nhAABA+ElP3yCpqasKPZ6YOMft/uTkzpKS0sXCkgEA4HtRNvt05GehXch1WbCmTZvKsWPHpEWLFrJjxw7H8QMHDsjo0aNl/vz5Esyys7OlcuXKcvjwYTNOHQAABL6l2xO0dAMAQjG29LilW7uOL1u2TBo2bGhat3XJMGfVq1cP+oAbAAAEngbOBM8AgEjhcdCtY7ZvuOEGGTNmjMTGxsrChQutLRkAAAAAAJESdHft2lX2798vv/32m5x33nnWlgoAAAAAgEgKulVUVJQJuLXP+r59+8y+mjVrmn7sAAAAAACgGOt0z549W5o1aybnnHOO+d/5vs5gDgAAAAAAitDS/dRTT0lKSoqMHDlSrrzySsea3NrlXCdYGzVqlPzxxx8yduxYT58SAAAAAICw5vGSYfXq1TOBd//+/d0eX7BggVmne9euXRLMWDIMAAAAAOCv2NLj7uW6DreuzV0YPaaTrAEAAAAAAC+D7ksuuUSeeOIJOXXqVIFjeXl5MnnyZJMGAAAAAAAUYZ1uHcuts5V36tTJZUz36tWrpXTp0mZsNwAAAAAA8HJMtzpy5Ii89tpr8vnnn7ssGdauXTsZOHBgSIyRZkw3AAAAAMBfsaVXQXc4IOgGAAAAAATdRGoAAAAAAMA7Pgu6k5KSpEGDBr56OgAAAAAAImcitbO55pprWDIMAAAAAAArgu4777zTV08FAAAAAEBYCIox3TNmzJD4+HgpW7astG3bVtavX19o2lmzZknHjh2latWq5qbd2s+UHgAAAACAoA+6Dxw44LK9ceNGGTx4sHTo0EH69esnK1euLFIBFixYIGPGjJHk5GT5+uuvpVWrVmY98Pz52Wk+119/vaxYsULWrVsncXFx0q1bN8nMzCxS/gAAAAAAWMXjJcNKliwpWVlZUr16dfnss8+kS5cu0r59e7n00ktNAK5B8PLly6VTp05eFUBbti+55BKZPn262T59+rQJpO+66y4ZP378WR+fl5dnWrz18YMGDTprepYMAwAAAAAUl6expcdjup1j85SUFLnpppvkpZdecuwbPXq0pKammsDbUydOnJANGzbIhAkTHPtKlChhuoxrK7Ynjh07JidPnpRzzjnH43wBAAAAAAjaMd3ff/+93HrrrS77dPvbb7/16nl0tnNtqa5Ro4bLft3et2+fR89x3333Sa1atUyg7s7x48fNFQjnGwAAAAAAQRd0HzlyxAStOuFZmTJlXI7pPm119qcnnnhC3njjDfnPf/5j8nfn8ccfN03+9pt2XQcAAAAAIOiC7gsuuMCMn/7ll1/kq6++cjn2ww8/mBZnb1SrVs2MFd+/f7/Lft2uWbPmGR/79NNPm6B72bJl0rJly0LTadd17WNvv+3evdurMgIAAAAAUFQej+nWidKcxcbGumzv3LlTbrvtNq8yL126tCQkJJhx4H369HFMpKbbI0aMKPRxTz75pEyaNEk+/PBDadOmzRnz0Bb5/K3yAAAAAAAEVdDduXPnMx4fNWpUkQqgy4Xp0mMaPOtM6NOmTZOcnBwZOnSoOa4zkteuXdt0E1eTJ0+WiRMnyvz5883a3vax3xUqVDA3AAAAAABCLuh2tmvXLrN8mM403qBBAzn33HOLXIABAwbIwYMHTSCtAXTr1q0lIyPDMbma5qX52L3wwgtm1nNdG9yZrvOts6oDAAAAABBy63Sr559/3rQ079mzx2V/u3btJC0tzXQVD3as0w0AAAAA8Fds6fFEajpxmY6jHjdunKSnp0vjxo1Ny/IHH3xgWrs7depUYHI1AAAAAAAimcct3fXr1zct3T169DDb27dvl/bt25su4aVKlTJjurds2WJmEw9mtHQDAPC3rKwjkpV11OvHxcZWkNjYipaUCQCAUOBpbOnxmO4DBw5I06ZNHdvnn3++eXIdj60zmd98882SmJhY/JIDAAC/SU/fIKmpq7x+XHJyZ0lJ6WJJmQAACCelvFmj+6OPPpJbb73VsYSYLvllX0+7bNmyEhUVZV1JAQCAzw0bliC9ejV22Zebe1ISE+eY+2vXDpWYmGi3Ld0AAMCHQfeECRPkxhtvlI8//tgE2IsXL5aRI0c6Au2VK1dK8+bNPX06AAAQBLSLeP5u4jk5Jxz3W7euKeXLlw5AyQAACA8eT6TWv39/eeedd8z4bV1He8qUKY61s5Uu4fXee+9ZVU4AAAAAAMJ7nW6dRM0+kVp+xVmrGwAAAACAiG7pBgAAAAAA3iHoBgAAAADAIgTdAAAAAABYhKAbAAAAAACLEHQDAAAAABDI2cv79u3r8RPq+t0AAAAAAMDDlu7KlSs7bpUqVZLly5fLV1995Ti+YcMGs0+PAwAAAAAAL1q658yZ47h/3333Sf/+/WXmzJlSsmRJsy8vL0/uuOMOE5ADAAAAAIAijul++eWXZezYsY6AW+n9MWPGmGMAAAAAAMCLlm5np06dkq1bt0rjxo1d9uu+06dPe/t0gF9kZR2RrKyjXj8uNraCxMZWtKRMAAAAAMKf10H30KFD5ZZbbpEdO3bIpZdeavZ98cUX8sQTT5hjQDBKT98gqamrvH5ccnJnSUnpYkmZAAAAAIQ/r4Pup59+WmrWrCnPPPOMZGVlmX2xsbEybtw4ueeee6woI1Bsw4YlSK9err0zcnNPSmLiX/MVrF07VGJiot22dIcaWvUBAACAEA66S5QoIffee6+5ZWdnm31MoIZgp8Fk/oAyJ+eE437r1jWlfPnSEg5o1QcAAABCOOh2RrANBJ9IatUHAAAAwi7o3r9/v5m9XNflPnDggNhsNpfjunwYgMCJpFZ9AAAAIOyC7iFDhsiuXbvkoYceMmO5o6KirCkZAAAIiLy8v1cjWb36V+nWraGULOn1KqMAAKAoQffatWtlzZo10rp1a2tKBIQJJjQDEIoWL94iI0cudWz37Dlf6tSpJGlp3aVv36YBLRsAABERdMfFxRXoUg6gICY0AxCKAXe/fgsl/898Zma22b9oUX8CbwAArA66p02bJuPHj5f09HSJj4/39uFAxGBCMwCh1qV81KiMAgG30n06mmz06Azp3bsxXc0BALAy6B4wYIAcO3ZMGjZsKOXKlZPoaNeg4ffff/f2KYGwxIRmAELJmjW7ZM+ev5YCdUcD7927s026Ll246A4AgKUt3QAAIPzmofBlOgAAUMSge/Dgwd4+BAAABDlPJ3BkokcAACwOup39+eefcuLE391lVaVKlYrzlAAAIAA6dqxrZinXSdPcjevWMd16XNMBAADPeT0TSk5OjowYMUKqV68u5cuXl6pVq7rcAABA6NHJ0XRZMHuA7cy+PW1adyZRAwDAS17/ct57773yySefyAsvvCBlypSR2bNnS2pqqtSqVUteeeUVb58OAAAECV0OTJcFq1XLtQu5tnCzXBgAAH7qXv7ee++Z4LpLly4ydOhQ6dixozRq1Ejq1asnr7/+utxwww1FLAoAAAg0DayTkupL5cqTzfaSJQOlW7eGtHADAFBEXv+C6pJgDRo0cIzfti8RlpiYKKtXry5qOQAAQJBwDrA7dapHwA0AQDF4/SuqAffOnTvN/SZNmsjChQsdLeBVqlQpTlkAAAAAAIjs7uXapXzTpk3SuXNnGT9+vFx99dUyffp0OXnypEyZMsWaUgLAGdYMzso66vXjYmMrsPQRAAAAgi/ovvvuux33k5KSZOvWrbJhwwYzrrtly5a+Lh+AEBGo4Dc9fYOkpq7y+nHJyZ0lJaVLkfMFAAAALF+nW+kEanoDENkCFfwOG5YgvXo1dtmXm3tSEhPnmPtr1w6VmJhot8F+cdDCDgAAAL8E3QAQyOBXA9j8QWxOzgnH/data0r58qXF12hhBwAAgCcIugH4RKCC30i7yAAAAIDQQtANAEUQaRcZAAAAUDQsvAkAAAAAQDC1dO/YsUPmzJlj/k9LS5Pq1avL0qVLpW7dunLhhRf6vpTwOSaBAgAAAIAgDLpXrVolPXr0kA4dOsjq1atl0qRJJujWtbtfeuklWbRokTUlRVhMAkWwDwAAACCSeB10jx8/Xh599FEZM2aMVKz4dxB0xRVXyPTp031dPoTZJFDM+AwAAAAgkngddH/33Xcyf/78Avu1tfu3337zVbkQppNAMeMzAAAAgEjiddBdpUoVycrKkvr167vs/+abb6R27dq+LBvCEDM+AwAAhCaGCQJ+Crr/7//+T+677z558803JSoqSk6fPi2ffvqpjB07VgYNGlTEYgAAAAAIZgwTBPwUdD/22GNy5513SlxcnOTl5UmzZs3M/wMHDpQHH3ywiMUAAJwNLQwAgEBimCDgp6C7dOnSMmvWLJk4caIZ33306FG56KKL5Pzzzy9iEQAAnqCFAQAQSAwTBPwUdD/88MOmK7m2dOvNLjc3V5566ikTjAMAfI8WBgAAgAgIulNTU+X222+XcuXKuew/duyYOUbQ7R26iwLwFC0M8NfvkF7Msdu4cV+hF3P4HQIAwIKg22azmQnU8tu0aZOcc8453j5dxKO7KAC4x0XJ4PgdsvekyI/fIQAAfBx0V61a1QTbervgggtcAm+dSE3HdmsLOLxDd1EAcI+LkoH7HfIEv0MAAPg46J42bZpp5b755ptNN/LKlSu7TK4WHx8v7dq18/Tp8D90FwUQ7ALV4sxFycD9DgEAgAAE3YMHDzb/169fX9q3by/R0QVPdAAA4SdQLc5clAQARCKGV4Ufr8d0d+7c2XH/zz//lBMn/j4BUpUqVfJNyQAAQYEWZwBAJApU8MvwqvDjddCts5Tfe++9snDhQvnvf/9b4LiO7wYAhA9anAEguNAS6h+BCn652B1+vA66x40bJytWrJAXXnhBbrrpJpkxY4ZkZmZKenq6PPHEE9aUEgAAP+FkFkCwoyXUPwIV/HKxO/x4HXS/99578sorr0iXLl1k6NCh0rFjR2nUqJHUq1dPXn/9dbnhhhusKSkAAH7AySyAYEdLqH8Q/CJgQffvv/8uDRo0cIzf1m2VmJgow4cP91nBAAAIBE5mAQQ7gkEgzINuDbh37twpdevWlSZNmpix3ZdeeqlpAa9SpYo1pQQAwE84mQUAAL5UwtsHaJfyTZs2mfvjx483Y7rLli0rd999txnvDQAAAAAAitjSrcG1XVJSkmzdulU2bNhgxnW3bNnS26cDAAAAACBseR1056cTqOkNAAAAAAD4IOj+8ssvzbJhBw4ckNOnT7scmzJlSlGeEgCAiMZSZQAAhCevg+7HHntMHnzwQWncuLHUqFFDoqKiHMec7wPBLi/v7wtGq1f/Kt26NZSSJb2e5gAAfIKlygAACE9eB91paWny8ssvy5AhQ6wpEeAHixdvkZEjlzq2e/acL3XqVJK0tO7St29TCTdcYACCH0uVAYB79ARCxAXdJUqUkA4dOlhTGsBPAXe/fgvFZnPdn5mZbfYvWtQ/rALvSLvAAIQqlioDAPfoCYSInL1clwmbNm2aNSUCLG7xHTUqo0DArXSfjpAYPTpDevduHBYtwZF2gQEAAIQfegIh4oLusWPHylVXXSUNGzaUZs2aSXS06wd88eLFviwf4FNr1uySPXuyCz2uwenu3dkmXZcu8RLKIu0CAwAACE/0BELEBd0jR440M5dffvnlcu655zJ5GkJuTJAv0wWzSLrAAAAAgKJj3HyQBd3z5s2Tt956y7R2A6HG00ohHCqPSLrAAAAAgKJj3HyQBd3nnHOO6VoOhKKOHeuaScR0TLO7btfacUOPa7pQF0kXGAAAAFB0jJsPsqA7JSVFkpOTZc6cOVKuXDlrSgVYRMcu66zdOomYBtjOgbd9pMS0ad0tG+Psz6W7IukCAwAAAIqOcfPW8vps/9lnn5WlS5dKjRo1pEWLFnLxxRe73IBgp7N166zdtWq5ViwagFo5m7fOJN6s2fMuS3fFx6eZ/VZeYFD5p17wxwUGAAAAAEVo6e7Tp481JQH8SAPrpKT6UrnyZLO9ZMlAS1udA7V0l/0Cg67TnZl5xOUCgwbcLBcGAAAABFnQrV3LgXDgHGB36lTP0i7lgVy6y98XGAAAAAD8jbNuIIiW7gr1CwwAAAAAitDSrTOWb9++XapVqyZVq1Y949rcv//+uydPCUQMlu4CAAAAIpdHQffUqVOlYsWKjvtnCroBuGLpLgAAACByeRR0Dx482HF/yJAhVpYHCDss3QUAAAC47+mZlXVUvKXrg4dSg5XXE6mVLFlSsrKypHr16i77//vf/5p9eXl5viwfEPICvTZ4JPLneugAAAAomvT0DZKausrrxyUnd5aUlC4StkG3zV1TnYgcP35cSpdmwXTAHZbu8h9dnk3fZ+f10PV91gsfVr/PBPsAAACeGzYsQXr1auyyLzf3pCQmzjH3164dKjEx0W5bukOJx0H3s88+a/7X8dyzZ8+WChX+fqHaur169Wpp0qSJNaUEwgBLd1kvUOuhBzrYBwAACEWxsRULdBPPyTnhuN+6dU0pXz70G3Y9Drp1AjV7S/fMmTNNN3M7beGOj483+wEULtKW7vJny28g10MPZLAPAACAMAm6d+7caf6//PLLZfHixWbpMAAIlpZfb9ZD79IlPiyCfQAAAAQ/r88AV6xY4RJwa9fyjRs3yh9//OHrsgEIUfaWX+fx684tv3o8XNZD9ybYB4rbW8R5GwAAhGnQPXr0aHnppZccAXenTp3k4osvlri4OFm5cqUVZYSfcHIHf7T8Km359fXnK1DroQcq2I9UkVRP6cWpZs2ed+ktEh+fZslFKwAAEERB95tvvimtWrUy99977z355ZdfZOvWrXL33XfLAw88YEUZ4Qec3CHUW37t66Hbl2HLT/fHxfl+PfRABfuRKJD1lL+D/UD0FgEAAEESdOt63DVr1jT3lyxZItddd51ccMEFcvPNN8t3331nRRlhMU7uEA4tv/b10FX+wNvK9dADFexHmkDWU/4O9gPVWwQAAFjD67PPGjVqyObNm03X8oyMDOnatavZf+zYMZcZzREaguHkLpK6i0aCQLb82tdDr1XL9bk1KLZqBvFABfuRJJD1VCCCfeYJAAAgvHh9Fjh06FDp37+/NG/e3KzZnZSUZPZ/8cUXrNMdggJ9cke39vAT6JZfDaw3b77Dsa3roe/cOcrSJbsCEexH0kWrQNVTgQr2mScAQCiJhN8hwO9Bd0pKisyePVtuu+02+fTTT6VMmTJmv7Zyjx8/vtgFgn8rr0Ce3NGtPTwFQ8tvINZDD0SwHyljnCNtZnrmCQAQKmg8ATxTpDPRfv36mYnT6tSp49g3ePBg6d27d1GeDgGsvAJ1chcM3dohYdnyG0j+DPYjaYxzpM1MH+jeIgDgCRpPAM95fEbYs2dPOXz4sGP7iSeekEOHDrlMsNasWTMvskYwVF6BOrkLdLd2SFi2/EaKSBvjHGkz0wdDbxEAOBMaTwDvePyL/eGHH8rx48cd24899pj8/vvvju1Tp07Jtm3bvMwega68AnVyx5jFyBCIbt6RINLGOEfizPSR2lsEQGig8QTwjsdnKLZ8Z1n5txG6lVcgTu4YswgUXaSNcY7UmenpLQIgWNF4AniHZqcgEejKy98nd4xZBIou0sY4R+rM9IreIgCCEY0ngHc8/vXW5cH0ln+fL8yYMUPi4+OlbNmy0rZtW1m/fv0Z07/55ptmeTJN36JFC1myZImEumCovPx5chfoFiQglEXaGOdInJkeAIIZjSeAhd3LhwwZIn379jW3P//8U26//XbH9s033yxFsWDBAhkzZowkJyfL119/La1atZIrr7xSDhw44Db9Z599Jtdff73ccsst8s0330ifPn3M7fvvv5dQFomVV6BbkIBQFYljnAONFmcA+BuNJ4B3PP4m6JJg1atXl8qVK5vbjTfeKLVq1XJs67FBgwZ5mb3IlClT5NZbb5WhQ4ea2c9nzpwp5cqVk5dfftlt+rS0NOnevbuMGzdOmjZtKo888ohcfPHFMn36dAllkVp50YIEFE0kjnEGAAQPGk8Az5XyNOGcOXPE106cOCEbNmyQCRMmOPaVKFFCkpKSZN26dW4fo/u1ZdyZtoy//fbbXuWdk5MjJUuWLLBf92m3ded0hdGyxsTEFCntsWPHCkxGd+WVdeW113rJuHGfyN69Rx37a9euKJMnX26O589Du/jrRQq73NxcOX268JmDy5cvX2janJwT+ldxvJby5Us7jmnPhry8PI+e92xptbz2oQk6I/6ff+Y68k1IqPa/bfdpdZb8wuj7q++z/bN18uTJM6b926n/va/u0+vnwf5ZOdvzOqfVdJremfN7/NdrKV1oWmdlypSRUqVKOR7nvJJAfqVLl5bo6GiXtPn/tvbX6pxW/2b6tyuMptP03qYVOX3G99g5rX4e9XNZGH0P9L1Q+v3R75E7f71efX9LnTWtt9/7M6XN/z6XKJHn0zoi//deT2iSkupL5cqP6l5ZvLi//OMf8SboteflyzrCuZ4aO/YTycryrJ4qbh3h/L3P/x6XKxdtWR3hnPav11TwO+TrOsKZ/bP+l7wz1lPFrSM8SWtVHXG2770VdYS3aX1VR1h9HmGX/3vvTdri1hGFpbXyPOJM33srzyM8TWtlHWH/3mva7GydT0OfO0qWLdvm+D2wso64+upGkpR0h1SuPNn83i9e3K/A75BVdcTfvwd/X+j1Rx3h7pzKX3VEYb9DVtURf7/W0n6tI3LcvMfBWkec6e/nwhZAmZmZ+te1ffbZZy77x40bZ7v00kvdPiY6Oto2f/58l30zZsywVa9e3W36P//803b48GHHbffu3SbPwm49e/Z0eXy5cuUKTdu5c2eXtNWqVSs0bZs2bVzS1qtXr9C0jRs3s4mkmNuSJdttTZteWGhafR5nmk9habV8zrT8haXV1+1M35czvW/O+vXrd8a0R48edaQdPHjwGdMeOHDAkfaOO+44Y9qdO3c60o4dO/aMab///nvb0aPH//c+F/4+6G39+vWO533yySfPmHbFihWOtNOnTz9j2kWL/uNIO2fOnDOmXbhwoSOt3j9TWn0uu/fff/+MabWMdlr2M6XV126n78mZ0iYnJ5t0f73HZ/676d/KTv+GZ0qrnwE7/WycKa1IK/P31TLoZ+5MafUz6yzY64hmzZo50v31Hp/n9zoiJqa8Sz3Vo0f41RF2+nkORB2h39+/66neYVlHKH2vA1FH6Gfr7+9R+NYRSrcLS8t5RGjXEXZ33jnDJnK3o17+66bbTS2vI/6upwYHpI4QaW/K4L86IsomEm8Taf6//6P8UkccPpxrE2nskq8/6giRyo7zKUUd4b6O0DjzTDxu6Q5Vjz/+uKSmpkooiTLfob/HDjpvA0Aw1lMhPsIHAELW4sVbZMaMgyJSKd8R3e4vIgsDVLJwpF3mdZhVZad9hyUz87Dlf+ORI5eKyPUu+YpkiMgWS/OGb0Rp5C0Bos3y2qS/aNEiMxma8/jxQ4cOyTvvvFPgMXXr1jXdy0ePHu3Yp5OwaffyTZs2FUiv3QScu8VkZ2dLXFyc7N27VypVqhSU3cKOHTsp1aunmftHj06QqKhTfukWpl05atR42tzfv3+sVK9e1S/dwg4fPuaSr3O3diu7heXmnpIKFR433ZD37x/jkq+V3cvtr/XQoQekcuVyfute7u49trrrqOZbocKk/73Hrn9bq7uX16gxxXQv1++QdkH2V/dy5/e5YsWylncd/es9/uvCorv32Kquo3+91mfNfX2PS5Y87Zeuo/nf4/POq+KXrqOHDuUUWk9Z2XX0+PHT/6un8mT//rsLrafoXv4Xupf/he7lkdO9PCqqhMTHp8mePdlu0+hL1+E/P/00QsqUKW1JHXHihO1/9dRp2b9/dKH1lBXdy/+ql0vI0aMPmXytrCP+/e+NcuON70r+r5J9ThPncey+rCOWLv1F+vVbWGi+OuSrd+8LLOte/td7XNr81ut77K/u5TXy/eYGax2hsaXOc3b48GG3saVdQFu69cuUkJAgy5cvdwTd+ofR7REjRrh9TLt27cxx56D7o48+Mvvd0S+L67i4v+gfzvmPVxhP0hQlrfMHvaATXqR15Truw9u0WsGWdvtanCuHs/Emrf5typePcsn3TCeV7v6WhX22/h5TfDalzphvUZ9XfzDsP1pOex2v1X6CXHjaQkpbqpTLYz1L6/q3dfda9UfA08+wN2n/GmtV2qP3WCsyT59XK8jC0+rrLeVh2oKKntb1fY6JKW1RHZHf2b8/vqkjnEX7rY5w/d67vsfOS1daVUdour/+dGd/n4tfR7g6ftz+e1DS43qqaHXE2VlVR3jzvfddHVH0tCoY0npTR/jvPKJw/qsjfJPW2zrCirSe1BErV/5SaMCtNI7as+eIrFu3V7p0ibekjvj7wkAJj+sp39QRf/8enD1t8b73p0/b5L77VhYIfJXu05+i0aMzpHfvxmY8u6++93l5p2XUqIwz5jt+/EoZMKBVgUlMfVNHFHyP/VNHRJ/xNzeY6ogzBf/OAj7FrLZaz5o1S+bNmydbtmyR4cOHmys+Opu50hnRnSdaGzVqlGRkZMgzzzwjW7dulZSUFPnqq68KDdIBAACAcJSVdcSn6eDemjW7znpxY/fubJMuHPKF7wV8TPeAAQPk4MGDMnHiRNm3b5+0bt3aBNU1atQwx3ft2uVoxlft27eX+fPny4MPPij333+/nH/++aZrefPmzQP4KgAAAAD/io2t6NN0CK6LG1xUCR8BD7qVtlIX1lK9cuXKAvuuu+46cwMAAAAiVceOdc262JmZ2W67IGv3Yz2u6RB6Fze4qBI+At69HAAAAID3dBxvWlp3l4m17Ozb06Z1LzDeF0W7uJH/PbbT/XFxvr+4Eah84Xt8AwEAAIAQpTNm68zZtWq5tnZqsOY8ozZC7+IGF1XCB38hAAAAIIRpYL158x2O7SVLBsrOnaPCOuDWmb3tVq/+1WU7nC5ucFElPBB0AwAAACHOubWzU6d6Yd36uXjxFmnW7HnHds+e88165bo/HC9uROJFlXATvt9GAAAAAGFFA+t+/RZKZqbrjN06mZzutzrwDtTFjUi6qBKO+GsBAAAACHrahXzUqAy3M7Xb940enWF5V3PAWwTdAAAAAILemjW7ZM+e7EKPa+C9e3e2SQcEk6BYpxsAgEiXlXVEsrKOuuzLzT3puL9x4z6JiYku8LjY2Aqs0QogYupJX6YD/IWgGwCAIJCevkFSU1cVejwxcY7b/cnJnSUlpYuFJQOA4ODpBUYuRCLYEHQDABAEhg1LkF69Gnv9OG3pBoBI0LFjXbNUlk6a5m5ct65drcc1HRBMCLoBAAgC2jJD6wwAFE5n7E5L625mKdcA2znw1m01bVp3ZvZG0OETCQAAACAk6NrUixb1l1q1XC9Sagu37mftat9xngV+9epfmRW+GAi6AQAAAIQMDaw3b77Dsb1kyUDZuXMUAbcP6XrnzZo979ju2XO+xMenWb4Oergi6AYAAAAQUpy7kHfqVI8u5T6kgbV24c/MdJ0FXsfS634Cb+8xphsRIZKW4omk1woAAADf0S7ko0ZluJ2oTvfp2PnRozOkd+/GXOjwAkE3IkIkLcUTSa8VAAAAvrNmzS7Zsye70OMaeO/enW3SdekS79eyhTKCbkSESFqKJ5JeKwAAAHzbY9KX6fAXgm5EhEhaiieSXisAAAB8x9NzSM41vUPQDQCAE+ZFAABEqo4d65rl13TSNHfjunVMtx7XdPAcQTcAAE6YFwEAEKl0crS0tO5mlnINsJ0Db91W06Z1ZxI1LxF0AwDghHkRAACRTNc7X7Sov4wcudRl2TBt4daAm/XQvUfQDQCAE+ZFAABEOg2sk5LqS+XKk832kiUDpVu3hn5p4c7LO+24v3r1r37L10oE3QCAoMTYagAAAsc50O3UqZ5fAt/Fi7eYFna7nj3nmxZ27fIeyi3sBN2ARQgYgOJhbDUAAJFj8eItZix5/gncdFI33a9d3kM18CboBixCwAAUD2OrAQCIDHl5p2XUqAy3M6brPp3EbfToDOndu3FIdjUn6AYsQsAAFA9jqwEAiAxr1uySPXuyCz2ugffu3dkmXZcu8RJqCLoBixAwAAAAAJ4Ny/RlumBD0A0gpAVq7Dxj9gEAAHwj1sNzo1A9hyLoBhDSAjV2njH7AAAAvtGxY10zS7lOmuZuXLeO6dbjmi4UEXQDCGmBGjvPmH0AAADfKFmyhFkWTGcp1wDbOfDWbTVtWveQnERNEXQDCGmBGjvPmH0AAADf6du3qVkWTNfpzsz8e+y2tnBrwB2qy4Upgu4IxXhUAJ6ivgAAAP7Qt29TSUqqL5UrTzbbS5YMlG7dGoZsC7cdQXeEYjwqfI3ALHxRXwAAAH8p6RRgd+pUL+QDbkXQHaGBCuNR4WsEZuGL+gIAAKDoCLojNFBhPCp8jcAsfC/SUV8AAAAUHUF3gBGoIFwQmFmP3gQAgEBiKBlQNATdAUagAsBTXKQDAAQSF3+BoiHohl9xhRQoOi7SAQACiYu/QNEQdMOvuEIKAAAQmrj4CxQNQTf8iiukAAAAACIJQTf8iiukAAAAACJJ6K80DgAAAABAkCLoBgAAAADAIgTdAAAAAABYhKAbAAAAAACLEHQDAAAAAGARZi8HAAAAQkhW1hHJyjrqsi8396Tj/saN+yQmJtrtEqysIgP4H0E3AAAAEELS0zdIauqqQo8nJs5xuz85ubOkpHSxsGQA3CHoBgAAAELIsGEJ0qtXY68fpy3d8Bw9CuArBN0AAABACNGAjqAufHsUEOyHH4JuAAAAAAiSHgUMHwg/BN0AAAAAECQ9Chg+EH4IugEAAAAgSDB8IPywTjcAAAAAABYh6AYAAAAAwCIE3QAAAAAAWIQx3QAAAACCFktoIdQRdAMAAAAIWiyhhVBH0A0AAAAgaLGEFkIdQTcAAACAoMUSWgh1TKQGAAAAAIBFCLoBAAAAALAIQTcAAAAAABYh6AYAAAAAwCIE3QAAAAAAWISgGwAAAAAAixB0AwAAAABgEYJuAAAAAAAsQtANAAAAAIBFSln1xAAAAACA4JeVdUSyso667MvNPem4v3HjPomJiS7wuNjYChIbW9EvZQxlBN0AAAAAEMHS0zdIauqqQo8nJs5xuz85ubOkpHSxsGThgaAbAAAAACLYsGEJ0qtXY68fpy3dODuCbgAAAACIYNpFnG7i1mEiNQAAAAAALELQDQAAAACARQi6AQAAAACwCEE3AAAAAAAWIegGAAAAAMAiBN0AAAAAAFiEoBsAAAAAAIuwTjcAAAAAwO+yso5IVtZRl325uScd9zdu3CcxMdEFHhcbWyGk1hUn6AYAAAAA+F16+gZJTV1V6PHExDlu9ycnd5aUlC4SKgi6AQAAAAB+N2xYgvTq1djrx2lLdygh6AYAAAAA+F1sbMWQ6iZeVEykBgAAAACARQi6AQAAAACwCEE3AAAAAAAWIegGAAAAAMAiBN0AAAAAAFiEoBsAAAAAAIsQdAMAAAAAYBGCbgAAAAAALELQDQAAAACARQi6AQAAAACwCEE3AAAAAAAWIegGAAAAAMAiBN0AAAAAAFiEoBsAAAAAAIsQdAMAAAAAYBGCbgAAAAAALELQDQAAAACARQi6AQAAAACwCEE3AAAAAAAWKSURxmazmf+zs7MDXRQAAAAAQIiyx5T2GLMwERd0HzlyxPwfFxcX6KIAAAAAAMIgxqxcuXKhx6NsZwvLw8zp06dl7969UrFiRYmKipJgvmqiFwZ2794tlSpVIt8wyjeSXiv58pkiX/IN9jzJl3zDJU/yJd9wyTOQ+XpLQ2kNuGvVqiUlShQ+cjviWrr1zahTp46ECv2QBeKDRr7hmSf5hne+kfRayTe8842k10q+4Z1vJL1W8g3vfCPptXrrTC3cdkykBgAAAACARQi6AQAAAACwCEF3kCpTpowkJyeb/8k3vPKNpNdKvuGbJ/mSb7jkSb7kGy55ki/5hkuegczXKhE3kRoAAAAAAP5CSzcAAAAAABYh6AYAAAAAwCIE3QAAAAAAWISgO0jNmDFD4uPjpWzZstK2bVtZv369pfmtXr1arr76arOwe1RUlLz99ttitccff1wuueQSqVixolSvXl369Okj27ZtszzfF154QVq2bOlY969du3aydOlS8bcnnnjCvNejR4+2NJ+UlBSTj/OtSZMm4g+ZmZly4403yrnnnisxMTHSokUL+eqrryzNU783+V+v3u68807L8szLy5OHHnpI6tevb15nw4YN5ZFHHhF/TJlx5MgR8xmqV6+eybt9+/by5Zdf+rV+0Nc5ceJEiY2NNWVISkqSH3/80fJ8Fy9eLN26dTOfLz2+cePGYud5tnxPnjwp9913n/ksly9f3qQZNGiQ7N2719J87d9l/e5qvlWrVjXv8xdffGFpnkOGDCnwXerevXux8vQkX7Vlyxbp1auXWf9UX7P+XuzatcvSfN3VHXp76qmnLM336NGjMmLECKlTp475DjVr1kxmzpxZrDw9yXf//v3mb6zHy5UrZ/62xf3uevLb/ueff5o6Wb+7FSpUkGuvvdaUxep8X3zxRenSpYv57df349ChQ8XK05N8f//9d7nrrrukcePG5m9bt25dGTlypBw+fNjSfNWwYcPM75Hme95550nv3r1l69atlubp/LvQo0cPn5xTepKv/l3zf29vv/12y/NV69atkyuuuMLUU/rZ6tSpk+Tm5lqW7y+//FJoXfXmm29a+nr37dsnN910k9SsWdO83osvvljeeustS/PcsWOHXHPNNeYzrO9v//79i11fBAJBdxBasGCBjBkzxszY9/XXX0urVq3kyiuvlAMHDliWZ05OjslHg31/WbVqlfnR/fzzz+Wjjz4yJ7N6Aq1lsZKe1GjAu2HDBhMAakWpP0Q//PCD+IsGRenp6Sb494cLL7xQsrKyHLe1a9danucff/whHTp0kOjoaHNRY/PmzfLMM8+YYMHq99b5tepnS1133XWW5Tl58mRzMWf69OkmUNDtJ598Up577jmx2r/+9S/zGl999VX57rvvzHdIgzG94OGv+kFf67PPPmuCBA0C9YdY6yw9sbYyXz2emJho3m9fOlO+x44dM/WyXmTR/zXw1xMEDQ6tzFddcMEF5jOmf2f9DusFJv17Hzx40LI8lQZizt+pf//730XOz9N89SRL/7Z6kWHlypXy7bffmvdcL0Rbma/z69Tbyy+/bE5kNSi0Ml/9zc/IyJDXXnvN1CF6IU2D8HfffdeyfDUo0hPcn3/+Wd555x355ptvzMU7rT+K8zvsyW/73XffLe+9954JEDS9XrTq27dvkfP0NF/9/urn+f777y9WXt7kq69Nb08//bR8//33MnfuXPO3vuWWWyzNVyUkJMicOXPMZ+rDDz80f3NNoxeKrcrTbtq0aea74wue5nvrrbe6fH/1t8nqfDXg1s+U7tcGMj0H0e9uiRIlLMs3Li6uQF2VmppqLmDphQ4rX69eZNbfPK2b9LdIv7caBGv9YUWeOTk5Zls/S5988ol8+umncuLECXMx8fTp0xJSdPZyBJdLL73Udueddzq28/LybLVq1bI9/vjjfslfPxb/+c9/bP524MABk/eqVav8nnfVqlVts2fP9kteR44csZ1//vm2jz76yNa5c2fbqFGjLM0vOTnZ1qpVK5u/3XfffbbExERboOn727BhQ9vp06cty+Oqq66y3XzzzS77+vbta7vhhhtsVjp27JitZMmStvfff99l/8UXX2x74IEH/FI/6Ptas2ZN21NPPeXYd+jQIVuZMmVs//73vy3L19nOnTvN8W+++cZn+XmSr9369etNul9//dWv+R4+fNik+/jjjy3Lc/DgwbbevXv75Pm9yXfAgAG2G2+80e/55qev/YorrrA83wsvvND28MMPW/o9zp/vtm3bzL7vv//e5XzjvPPOs82aNcuy33atH6Kjo21vvvmmI82WLVtMmnXr1lmWr7MVK1aYY3/88YfP8vMkX7uFCxfaSpcubTt58qRf8920aZNJ89NPP1map9bFtWvXtmVlZVlyTukuX3+cT7nLt23btrYHH3zQ7/nm17p16wLnIVbkW758edsrr7ziku6cc87xWZ1xIF+eH374oa1EiRLm985O65CoqChzHh1KaOkOMnr1Rltg9UqznV4t0229mhbO7F2tzjnnHL/lqVd733jjDXMlTbuZ+4Ne0bvqqqtc/sZW0+6C2n2wQYMGcsMNNxS7i6Yn9CpomzZtTAuzdhm66KKLZNasWeLv75O2HN18880+u+LujnbpXr58uWzfvt1sb9q0ybREFueKsydOnTplPsP5W/+0K6E/ejOonTt3mu5mzp9n7Q6sw2LCvc5yrrv081WlShW/fra1q6y+19qSaSVtadbvsHaPHT58uPz3v/+1ND9tvfjggw9My772mNC89fPkj2FPzrT7opajuC2SntYhWmdqDxWNj1esWGHqE23hscrx48fN/871h55v6Jq4vqw/8v+26zmOtmY51xnao0G7XfuyzgjEOYWn+Woa7SZbqlQpv+Wr5zna6q3DoLSl1Ko8tTfBwIEDTe8K7YJshcJe6+uvvy7VqlWT5s2by4QJE0xZrMxXe6Bq7y6to/Q7XKNGDencubPPf3/P9rfV75QOr/J1XeUuX32d2iNXh01oXa3n0NqrTbv3W5Hn8ePHze+r81rdWmdpXeWv8xyfCXTUD1eZmZnmCs9nn33msn/cuHGmBTxcW7r16rq2Fnbo0MEv+X377bfmap22ElauXNn2wQcf+CVfbflr3ry5LTc3129XZpcsWWKuqusV7oyMDFu7du1sdevWtWVnZ1uar7Z06m3ChAm2r7/+2paenm4rW7asbe7cuTZ/WbBggfkb6/fK6s+vtuzrlddSpUqZ/x977DGbP+jfUz9H+hpPnTple/XVV81V4QsuuMAv9cOnn35q9u3du9cl3XXXXWfr37+/ZfkGS0u3fpe1RXLgwIF+yfe9994zdZd+xrQHlLayW5mn1lnvvPOOqTP1WNOmTW2XXHKJ+axZla+9daxcuXK2KVOmmL+r9vTS17xy5UrL8s1v8uTJpheUvb62Mt8///zTNmjQIHNM6xBtBZ03b56l+Z44ccL8Fuh39ffff7cdP37c9sQTT5h03bp1s+y3/fXXXzevLz/9XN17772W5euPlm5PzmUOHjxo3vf777/fL/nOmDHD1Bn6ehs3buyzVu7C8rzttttst9xyi2XnlIXlq+cYeo6jddVrr71mWtqvueYaS/PVnhn6+rSl9+WXXzbnOqNHjzaf7+3bt1uWb37Dhw83dbMvFZavfme0frDXVZUqVTKt0VbleeDAAZOHnivn5OTYjh49ahsxYoTJXz9roYSgO8hEatB9++232+rVq2fbvXu3X/LTk4sff/zR9tVXX9nGjx9vq1atmu2HH36wNM9du3bZqlevboJfO38E3flphakVmNXd6bX7oAaEzu666y7bZZddZvMX/WH45z//aXk+GpjUqVPH/K8/+Nr1Sn+E/XGBQU+gOnXqZL63eoFBT1y1W3uTJk0syY+g2+YStFx99dW2iy66yKXrm5X56gmH1l16sqddCePj42379++3NE9nO3bs8GmXdnf52n8Hr7/+epd0+l7/3//9n2X55qcBip7c+Zq7fHV4hl4oe/fdd81vxHPPPWerUKGCT7tPustXfwN1+JG9/rjyyittPXr0sHXv3t2y33Z/BN1nO6ewKug+W75aT+i5nL6/Wn/4I1/tiqsBoHbX1e+QXiT0xYUkd3nqBbpGjRqZYXRWnVN6er64fPlyn3ald5ev/fdPGxectWjRwpxbWpVv/mFm2nj09NNP+yS/s+WrdaJ+hvU3YOPGjbaUlBSTv577WJXnhx9+aGvQoIG58Kr1lA490s+xpg8lBN1BRoNB/UDlr6D0CnivXr3CMujW8esasPz888+2QPnHP/5h+RUzfU/tJzb2m27bKxFfthydTZs2bXz2g1AYvZLvfLVbPf/886Z1zh9++eUX0+L79ttvW56Xfn6nT5/usu+RRx4xJ+3+osGYPfDVYLdnz55+qR/sQVj+gFcvBIwcOdKyfAMddOsJc58+fWwtW7a0/fbbb37LNz89wfVVrwpP89SLlDNnzvRJnu7y1d9BbUHR75AzDcjat29vWb7OVq9ebY7rSaWv5c9XT5r1ImX+uRm0/tQg2Kp88wdm2qKk9IT6jjvusOy33R4M5Q949TdDezZYla/VQffZ8tXeZXohWs83fNl7wptzKP1uaQ+S+fPnW5KnNiLYz2mcz3P0t1gbGfz5WvU3UfPW1m+r8tVtzUN7mDnT32Bf9H7y5PXqRX6tP+zfX18oLF+9gJF/Hgiln+lhw4ZZkmf+XiL272yNGjVsTz75pC2UMKY7yJQuXdrMNqnjQ+10zIRu+2vMsb/oOYDO8Pif//zHzEio44wCRd9j+/g2q/zjH/8wMz3quBv7Tcc86xhrvV+yZEnxB12aRmcG1uWdrKQzl+df9kHHKOrsuP6gY9d0nJWOn7eajhvLP1Op/j39ObOmzhiuf1OdNV5nqdUZ+f1Bv7c6bs+5zsrOzjbj3MKtzrLT8ag6W6vOlfDxxx+bZY/Cue5ytmfPHjOm28r6Q38HdQmZQNYfL730kvkttnq8vP3zpLdA1iE6N4Aux6OfaV3Vozj1x9l+2/V91VUtnOsM/VvrXCPFqTMCdU7hSb5aJ+r4fP1s69j94s7C72m+7h6jt6LWGWfLc/z48WalAefzHDV16lTzm+zP12rPuzh11dny1RUkdL4cX9dV3rxerat09Qz9/hbX2fK1j5H3ZV1l8+K16nh9nTtF0+l4el+sGuJXgY76UdAbb7xhxsJq19TNmzebFtgqVarY9u3bZ1me2hVIW4n0ph8L+zg6X87G624MinZJ0TF6OobPftOr/lbSFl7tZqWtY9odRrf1yuyyZcts/uaP7uX33HOPeY/19WpXqKSkJNNS5curou7oWFNtrZo0aZLpDqtdCvUKu461spqOC9JWEx1n7Q86w7OOH9OWKn2fFy9ebN5jX3WVPBO9ir906VJzdVg/w9pVVGdT9WXXxbPVDzoOVOso+9hfnfG5fv36xW7NOVu+//3vf822zsmgx7Xu1G2tR6zKV99X7XWkV+S1FdS57tJWJKvy1VYb7cKo3cq1F4d2Cx46dKj5rcjf6uCrPPXY2LFjTZ76udbuhNqlT1df0DHIVr1Wpd8hbb158cUXTf2h3a211WzNmjWW5mvvAqx11QsvvFCsvLzJV38LdAZzbYHV7/KcOXPMHBjaO8jKfHW+D81Te6xoryDt2qkrL1j9267dQrWO/uSTT8xnWVuA8w9HsiJf3dbXrzMt6/uhPRp0W+sSq/LVz5PWydrlWFsKndMUp4fb2fLVv6n2gtH3V//e+vuv3ct16FNRh6QU5bzNF70nz5avvq86+7++Vq2r9LdIuyNrjysr81VTp041Q/Z0Nn6tq3Qmc/3uFqdbu6fvs+an5696DuALZ8tXf/+0d1XHjh1tX3zxhXmN2q1dy1DUuZGGe/Badby8/g5pftqrQD/DY8aMsYUagu4gpScY+oOk4560q9fnn39uaX72rlb5bxpMWMVdfnrTkw0r6ThIPbHQ91aXRtFuMYEIuP0VdOvSO7Gxseb1amCo274a43Q2OumTThyngYGOMdYTaH/Q8T/6WdIlcfxBuw3q31G/s/pjqz/2utRPcYMwTyeL0/z076tLd2kXLe0q6s/6QZcNe+ihh0x3L/1b63fKF+/92fLVusLdcV0mz6p87V3Z3d30cVblqxcwdFIgHZ6hf2v9TmvwX9yJ1M6Up5706LwIWk9qAKz15q233uqTC8Ce/Oa89NJL5gRPv1N6MckXQ0U8yVcnZIqJifHp9+hs+epJ5pAhQ8zfV1+vDk155plnir3U4dnyTUtLMxeQ9O+r9ZcGDMWttzz5bdfPs3Zh14nq9AKHfraLe7HMk3y1bvD1ecfZ8i3sb6A3rU+sylfnRtDx+TqXjP599e+sXZ63bt1qWZ5WBd1ny1fnzNEAW4Mx/Q3SekPnQiruXBuevl6d6FHfX/0s68Wj4l4c9DRfvRAbFxdnGhp8wZN8dX4AvTCnnyt9vTrEKv8SYr7O87777jPnF/o51ou+vqgbAyFK/wl0azsAAAAAAOGIMd0AAAAAAFiEoBsAAAAAAIsQdAMAAAAAYBGCbgAAAAAALELQDQAAAACARQi6AQAAAACwCEE3AAAAAAAWIegGAAAAAMAiBN0AAESYffv2SdeuXaV8+fJSpUqVQvdFRUXJ22+/7dFzpqSkSOvWrS0tNwAAoYigGwCAIKLB71133SUNGjSQMmXKSFxcnFx99dWyfPlyn+UxdepUycrKko0bN8r27dsL3afbPXr08Og5x44d69Myqrlz5zouAAAAEKpKBboAAADgL7/88ot06NDBBJpPPfWUtGjRQk6ePCkffvih3HnnnbJ161af5LNjxw5JSEiQ888//4z7atas6fFzVqhQwdwAAIArWroBAAgSd9xxh+nSvX79ern22mvlggsukAsvvFDGjBkjn3/+uUmza9cu6d27twlwK1WqJP3795f9+/e7PM8777wjF198sZQtW9a0mKempsqpU6fMsfj4eHnrrbfklVdeMXkNGTLE7T533cv37Nkj119/vZxzzjmmG3qbNm3kiy++KLR7+ezZs6Vp06amHE2aNJHnn3/e5QKDPv/ixYvl8ssvl3LlykmrVq1k3bp15vjKlStl6NChcvjwYZNOb5oHAAChhpZuAACCwO+//y4ZGRkyadIkE9Dmp63fp0+fdgTcq1atMoG0toAPGDDABKlqzZo1MmjQIHn22WelY8eOpgX7tttuM8eSk5Plyy+/NMc1YE9LS5OYmBg5ceJEgX35HT16VDp37iy1a9eWd99917SCf/3116ZM7rz++usyceJEmT59ulx00UXyzTffyK233mpe2+DBgx3pHnjgAXn66adNC7ve16D+p59+kvbt28u0adPMc2zbts2kpSUdABCKCLoBAAgCGmjabDbTIlwYHTP93Xffyc6dO81Yb6Wt09oarsH0JZdcYlq1x48f7whstaX7kUcekXvvvdcE3eedd54ZK66BtXP3cXf7nM2fP18OHjxo8tGWbtWoUaNCy6p5PfPMM9K3b1+zXb9+fdm8ebOkp6e7BN06Fvyqq64y97Xs+lr0vdD3oXLlyqaF25tu7gAABBuCbgAAgoAG3GezZcsWE2zbA27VrFkz0wquxzTo3rRpk3z66aemxdwuLy9P/vzzTzl27Jjpxl0UOsGatljbA+4zycnJMS3st9xyi2ndttOWeQ2knbVs2dJxPzY21vx/4MCBM158AAAglBB0AwAQBLR7tbbqFneyNO0Gri3G9hZmZzq2uqjcdTk/UxnUrFmzpG3bti7HSpYs6bIdHR3tuK+vXxXWZR0AgFBE0A0AQBDQFuQrr7xSZsyYISNHjiwwrvvQoUNmUrLdu3ebm721W7ts6zFt8VY6gZqOgT5T1++i0BZpnRhNx56frbW7Ro0aUqtWLfn555/lhhtuKHKepUuXNq30AACEMmYvBwAgSGjArUHmpZdeamYT//HHH023cZ0UrV27dpKUlGSWEdNAVicx01nOdQI0neBMZxJXOvGYjvPW1u4ffvjBPP6NN96QBx98sFhl0wnOdGx1nz59TPd1Dai1jPbZxvPT/B9//HFTdl33W8eiz5kzR6ZMmeJxnjqruraa61j23377zXSPBwAg1BB0AwAQJHTSMw2mdQmte+65R5o3by5du3Y1QecLL7xgul/rcmBVq1aVTp06mSBcH7NgwQLHc2hr+fvvvy/Lli0zY7wvu+wymTp1qtSrV69YZdNWZ33O6tWrS8+ePU3w/8QTTxToLm73r3/9y7SMa6CtafXCwNy5c82Eap7SGcxvv/12Mzu7TgD35JNPFus1AAAQCFE2T2ZuAQAAAAAAXqOlGwAAAAAAixB0AwAAAABgEYJuAAAAAAAsQtANAAAAAIBFCLoBAAAAALAIQTcAAAAAABYh6AYAAAAAwCIE3QAAAAAAWISgGwAAAAAAixB0AwAAAABgEYJuAAAAAAAsQtANAAAAAIBY4/8BYDZu56MyDVsAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fit_panel_time = pf.feols(\n", + " \"minutes_watched ~ i(day, ever_treated, ref = 14) | user + day\",\n", + " vcov={\"CRV1\": \"user\"},\n", + " data=data,\n", + " demeaner_backend=\"rust\",\n", + ")\n", + "fit_panel_time.iplot(\n", + " coord_flip=False,\n", + " plot_backend=\"matplotlib\",\n", + " cat_template=\"{value}\",\n", + " title=\"Event Study Plot\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "bf038777", + "metadata": {}, + "source": [ + "Initially, we observe a large treatment effect that reverts to a null effect on day 23: a clear novelty effect pattern. Clearly, we should not accept this \n", + "experiment as it only adds technical complexity, but has no impact on user behavior beyond the initial effect. When ignoring time heterogeneity (as with the difference in means, cuped, and ATE panel estimators above), we completely miss that the effect fades out. With a dynamic panel estimator, we get it almost for free. \n", + "\n", + "For more examples, see the paper by [Lal et al](https://arxiv.org/abs/2410.09952), which also shows how to estimate very large panel models in [SQL/duckdb](https://github.com/py-econometrics/duckreg) (yep, that's not a joke!). " + ] + }, + { + "cell_type": "markdown", + "id": "41afbbe7", + "metadata": {}, + "source": [ + "## Monte Carlo Simulation" + ] + }, + { + "cell_type": "markdown", + "id": "2a4d8a02", + "metadata": {}, + "source": [ + "To rule out that the examples above were purely do to look, we simply repeat them a couple of times using a monte carlo simulation. " + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "edd54613", + "metadata": {}, + "outputs": [], + "source": [ + "def _variance_monte_carlo(data):\n", + " fit_dim = _difference_in_means(data)\n", + " fit_cuped, _ = _regression_cuped(data)\n", + " fit_panel = pf.feols(\n", + " \"minutes_watched ~ treat | user + day\",\n", + " data=data,\n", + " vcov={\"CRV1\": \"user\"},\n", + " demeaner_backend=\"rust\",\n", + " )\n", + "\n", + " dim_se = fit_dim.tidy().xs(\"treat\")[\"Std. Error\"]\n", + " cuped_se = fit_cuped.tidy().xs(\"treat\")[\"Std. Error\"]\n", + " panel_se = fit_panel.tidy().xs(\"treat\")[\"Std. Error\"]\n", + "\n", + " return dim_se, cuped_se, panel_se\n", + "\n", + "\n", + "def _run_simulation(N, sigma_unit, n_sim=100):\n", + " dim_se_arr = np.zeros(n_sim)\n", + " cuped_se_arr = np.zeros(n_sim)\n", + " panel_se_arr = np.zeros(n_sim)\n", + "\n", + " for i in tqdm(range(n_sim)):\n", + " data = get_sharkfin(\n", + " num_units=N,\n", + " num_periods=30,\n", + " num_treated=int(N / 2),\n", + " treatment_start=15,\n", + " seed=i,\n", + " sigma_unit=sigma_unit,\n", + " )\n", + " data.rename(\n", + " columns={\"Y\": \"minutes_watched\", \"unit\": \"user\", \"year\": \"day\"},\n", + " inplace=True,\n", + " )\n", + "\n", + " dim_se_arr[i], cuped_se_arr[i], panel_se_arr[i] = _variance_monte_carlo(data)\n", + "\n", + " return pd.Series(\n", + " {\n", + " \"Difference-in-Means\": np.mean(dim_se_arr),\n", + " \"Cuped\": np.mean(cuped_se_arr),\n", + " \"Panel\": np.mean(panel_se_arr),\n", + " }\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b0fd4eb", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 20/20 [01:30<00:00, 4.53s/it]\n" + ] + }, + { + "data": { + "text/plain": [ + "Difference-in-Means 0.029457\n", + "Cuped 0.004184\n", + "Panel 0.004185\n", + "dtype: float64" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# note that in a proper scientific simulation, we might want to set the\n", + "# number of iterations to be much higher than 10\n", + "se_sim_18 = _run_simulation(N=100_000, sigma_unit=18, n_sim=20)\n", + "display(se_sim_18)" + ] + }, + { + "cell_type": "markdown", + "id": "c24d0c3e", + "metadata": {}, + "source": [ + "So we manage to reduce our standard errors by around 6x! This is pretty fantastic news, as it implies that we can run our experiment on a much smaller sample, but still achieve the same level of power! This is what we want to verify in the last step. " + ] + }, + { + "cell_type": "markdown", + "id": "671ca0d9", + "metadata": {}, + "source": [ + "## Power Simulation \n", + "\n", + "In all simulations, we conduct a two-sided t-test with size $\\alpha = 0.01$. " + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "e04bee49", + "metadata": {}, + "outputs": [], + "source": [ + "def _power(nobs, n_sim):\n", + " res_df = pd.DataFrame()\n", + "\n", + " for N in nobs:\n", + " dim_reject_null_arr = np.zeros(n_sim)\n", + " cuped_reject_null_arr = np.zeros(n_sim)\n", + " panel_reject_null_arr = np.zeros(n_sim)\n", + "\n", + " for i in tqdm(range(n_sim)):\n", + " data = get_sharkfin(\n", + " num_units=N,\n", + " num_periods=30,\n", + " num_treated=int(N / 2),\n", + " treatment_start=15,\n", + " seed=i,\n", + " sigma_unit=18,\n", + " )\n", + " data.rename(\n", + " columns={\"Y\": \"minutes_watched\", \"unit\": \"user\", \"year\": \"day\"},\n", + " inplace=True,\n", + " )\n", + "\n", + " fit_dim = _difference_in_means(data)\n", + " fit_cuped, _ = _regression_cuped(data)\n", + " fit_panel = pf.feols(\n", + " \"minutes_watched ~ treat | user + day\", data=data, vcov={\"CRV1\": \"user\"}\n", + " )\n", + "\n", + " dim_reject_null_arr[i] = fit_dim.pvalue().xs(\"treat\") < 0.01\n", + " cuped_reject_null_arr[i] = fit_cuped.pvalue().xs(\"treat\") < 0.01\n", + " panel_reject_null_arr[i] = fit_panel.pvalue().xs(\"treat\") < 0.01\n", + "\n", + " dim_reject_null_mean = np.mean(dim_reject_null_arr)\n", + " cuped_reject_null_mean = np.mean(cuped_reject_null_arr)\n", + " panel_reject_null_mean = np.mean(panel_reject_null_arr)\n", + "\n", + " res = pd.DataFrame(\n", + " {\n", + " \"N\": [N],\n", + " \"Difference-in-Means\": [dim_reject_null_mean],\n", + " \"Cuped\": [cuped_reject_null_mean],\n", + " \"Panel\": [panel_reject_null_mean],\n", + " }\n", + " )\n", + "\n", + " res_df = pd.concat([res_df, res], axis=0, ignore_index=True)\n", + "\n", + " return res_df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e9265641", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 10/10 [00:30<00:00, 3.01s/it]\n", + "100%|██████████| 10/10 [00:07<00:00, 1.41it/s]\n", + "100%|██████████| 10/10 [00:06<00:00, 1.46it/s]\n", + "100%|██████████| 10/10 [00:10<00:00, 1.02s/it]\n", + "100%|██████████| 10/10 [00:12<00:00, 1.26s/it]\n", + "100%|██████████| 10/10 [00:53<00:00, 5.31s/it]\n" + ] + } + ], + "source": [ + "# TODO: set n_sim higher\n", + "power_df = _power(nobs=[100, 500, 1000, 5_000, 10_000, 100_000], n_sim=10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6211f87d", + "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", + "
NDifference-in-MeansCupedPanel
01000.50.10.1
15000.40.90.9
210000.71.01.0
350000.51.01.0
4100000.61.01.0
51000000.91.01.0
\n", + "
" + ], + "text/plain": [ + " N Difference-in-Means Cuped Panel\n", + "0 100 0.5 0.1 0.1\n", + "1 500 0.4 0.9 0.9\n", + "2 1000 0.7 1.0 1.0\n", + "3 5000 0.5 1.0 1.0\n", + "4 10000 0.6 1.0 1.0\n", + "5 100000 0.9 1.0 1.0" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "display(power_df)" + ] + }, + { + "cell_type": "markdown", + "id": "85cecc4d", + "metadata": {}, + "source": [ + "CUPED and the panel estimator achieve 90% power with 500 observations, while the difference in means-estimator requires a sample that is orders of magnitude larger. \n", + "To conclude, we should mention that these effects stem from a stylized exmaple data set - in practice, the gains from CUPED / panel methods for variance reduction are much lower - \n", + "companies with access to high quality panel data report reductions of 10% to 50% (see e.g. this paper by [Microsoft Bing](https://exp-platform.com/Documents/2013-02-CUPED-ImprovingSensitivityOfControlledExperiments.pdf))." + ] + }, + { + "cell_type": "markdown", + "id": "9f2394bc", + "metadata": {}, + "source": [ + "## When does CUPED not work? \n", + "\n", + "CUPED only works when we explain a lot of \"residual\" variation, which in our example data set is controlled by the `sigma_unit` parameter. Let's see what happens if we drastically reduce it - we now change the parameter from `18` to `4`." + ] + }, + { + "cell_type": "markdown", + "id": "8a776941", + "metadata": {}, + "source": [ + "We repeat the standard error monte carlo simulation, but set `sigma_unit = 4`." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "9e308933", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 10/10 [00:48<00:00, 4.89s/it]\n" + ] + } + ], + "source": [ + "se_sim_1 = _run_simulation(N=100_000, sigma_unit=4, n_sim=10)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "16388f84", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Difference-in-Means 0.006852\n", + "Cuped 0.004174\n", + "Panel 0.004186\n", + "dtype: float64" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "display(se_sim_1)" + ] + }, + { + "cell_type": "markdown", + "id": "00fdf6ac", + "metadata": {}, + "source": [ + "All of a sudden, the difference-in-means estimator has lower variance!" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "f9b04373", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Difference-in-Means 0.029457\n", + "Cuped 0.004184\n", + "Panel 0.004185\n", + "dtype: float64" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "display(se_sim_18)" + ] + }, + { + "cell_type": "markdown", + "id": "d78b86e9", + "metadata": {}, + "source": [ + "The reason we do not achieve variance reduction this time is that the addition of pre-treatment covariates / unit fixed effects simply explains less unobserved variance than before. " + ] + }, + { + "cell_type": "markdown", + "id": "e6f98862", + "metadata": {}, + "source": [ + "## Other useful links\n", + "\n", + "- [Matteo Courthoud's great blog post on CUPED](https://matteocourthoud.github.io/post/cuped/), goes through some theory, compares CUPED with other estimators, runs simulations studies, and summarizes the literature. \n", + "- [Lal et al on panel experiments](https://arxiv.org/abs/2410.09952) explains why using panel estimators might be a good idea when analysing AB tests, and shows how it can be done efficiently in SQL. " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "docs", + "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.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}