{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"(Disclaimer: I work for Pro Football Focus, one of the data provider for this competition and hence not eligible for the prize)\n\n## What contribute to successful field goal?\n\n\nMultiple factors can affect field goal percentage in NFL and some of them might not be what you think. The following use [PyMC](https://github.com/pymc-devs/pymc) to construct Bayesian model to estimate field goal success percentage across different factors, with [Jax](https://github.com/google/jax) to speed up inference time and make Bayesian model more accessible.\n\nWith the previous notebook https://www.kaggle.com/s903124/bayesian-field-goal-model-with-pymc , it has established that field goal data with field goal distance and angle perform the best, and the below expand beyond the base model.\n\nThe variable field goal success is estimated by Binomail distribution\n\n$$ y_i \\sim Binomial(n, p_i)$$\n\nwhere i is success rate of individual field goal. In base model field goal success depend on distance and angle only, and therefore\n\n$$ p_i \\sim InverseLogit( \\theta \\cdot angle_i + d \\cdot distance_i + intercept)$$\n\nwhere \n\n$$ \\theta \\sim Normal(0, 10)$$\n$$ d \\sim Normal(0, 10)$$\n$$ intecept \\sim Normal(0, 1)$$\n","metadata":{}},{"cell_type":"code","source":"!pip install numpyro\n!git clone https://github.com/pymc-devs/pymc/ && cd pymc && pip install .","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-01-05T12:32:26.786251Z","iopub.execute_input":"2022-01-05T12:32:26.786882Z","iopub.status.idle":"2022-01-05T12:33:42.49392Z","shell.execute_reply.started":"2022-01-05T12:32:26.786767Z","shell.execute_reply":"2022-01-05T12:33:42.492853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport pymc as pm\nimport pymc.sampling_jax\nimport arviz as az\nfrom aesara import tensor as aet\n\npd.options.display.max_columns = 999\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\nfrom patsy import dmatrix\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-01-05T12:39:10.483526Z","iopub.execute_input":"2022-01-05T12:39:10.484354Z","iopub.status.idle":"2022-01-05T12:39:17.719871Z","shell.execute_reply.started":"2022-01-05T12:39:10.4843Z","shell.execute_reply":"2022-01-05T12:39:17.719275Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"play_data = pd.read_csv('/kaggle/input/nfl-big-data-bowl-2022/plays.csv')\nfield_goal_data = play_data[play_data.specialTeamsPlayType == 'Field Goal'][['gameId','playId','absoluteYardlineNumber','specialTeamsResult','playDescription']]\npff_data = pd.read_csv('/kaggle/input/nfl-big-data-bowl-2022/PFFScoutingData.csv')\nstadium_data = pd.read_csv('../input/weather-data/stadium_coordinates.csv')\nweather_data = pd.read_csv('../input/weather-data/games_weather.csv')\ngame_data = pd.read_csv('../input/weather-data/games.csv')\n\nfield_goal_data = pd.merge(field_goal_data,game_data,left_on='gameId',right_on='game_id')\nfield_goal_data = pd.merge(field_goal_data,stadium_data)","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-01-05T12:39:17.721312Z","iopub.execute_input":"2022-01-05T12:39:17.721637Z","iopub.status.idle":"2022-01-05T12:39:18.072842Z","shell.execute_reply.started":"2022-01-05T12:39:17.72161Z","shell.execute_reply":"2022-01-05T12:39:18.071996Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Load tracking data\n\ntracking_data = []\n\nfor year in range(2018,2021):\n    data = pd.read_csv('/kaggle/input/nfl-big-data-bowl-2022/tracking'+str(year) + '.csv')\n    data = data[data.event == 'field_goal_attempt']\n    tracking_data.append(data)\ntracking_data = pd.concat(tracking_data)\ndel data\n\ntracking_data = pd.merge(tracking_data,field_goal_data)\ntracking_data.loc[tracking_data.playDirection == 'left','x'] = 120-tracking_data['x']\ntracking_data.loc[tracking_data.playDirection == 'left','y'] = 53.33-tracking_data['y']","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-01-05T12:39:18.074337Z","iopub.execute_input":"2022-01-05T12:39:18.074621Z","iopub.status.idle":"2022-01-05T12:41:04.170959Z","shell.execute_reply.started":"2022-01-05T12:39:18.074586Z","shell.execute_reply":"2022-01-05T12:41:04.169212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Merge with weather data\n\nweather_data['TimeMeasure'] = pd.to_datetime(weather_data['TimeMeasure'])\ntracking_data['time'] = pd.to_datetime(tracking_data['time'])\ntracking_list = []\nfor game_id in tracking_data['gameId'].unique():\n    tracking_list.append(pd.merge_asof(tracking_data[tracking_data.gameId == game_id],weather_data[weather_data.game_id == game_id],left_on='time',right_on='TimeMeasure',direction='nearest'))\ntracking_data = pd.concat(tracking_list)\nfield_goal_ball_df = tracking_data[tracking_data.team == 'football']","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-05T12:41:04.17512Z","iopub.execute_input":"2022-01-05T12:41:04.175509Z","iopub.status.idle":"2022-01-05T12:41:10.312785Z","shell.execute_reply.started":"2022-01-05T12:41:04.175456Z","shell.execute_reply":"2022-01-05T12:41:10.311793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Calculate distance and angle\n\nba = np.array(np.array([120,23.583])-field_goal_ball_df[['x','y']])\nbc = np.array(np.array([120,29.75])-field_goal_ball_df[['x','y']])\n\nfield_goal_ball_df['angle'] =  np.degrees(np.arccos(np.array([np.dot(a,b) for a,b in zip(ba,bc)])/(np.linalg.norm(ba,axis=1)  * np.linalg.norm(bc,axis=1) )))\nfield_goal_ball_df['fg_dist'] = ((field_goal_ball_df['x'] - 120)**2 + (field_goal_ball_df['y'] - 26.33)**2)**0.5\nfield_goal_ball_df = field_goal_ball_df[field_goal_ball_df.fg_dist <= 70]\nfield_goal_ball_df['fg_make'] = np.array(field_goal_ball_df[\"specialTeamsResult\"] == 'Kick Attempt Good').astype(int)\n\nfield_goal_ball_df['kickerName'] = field_goal_ball_df['playDescription'].str.split('.',expand=True)[0].str.rsplit(' ',n=1,expand=True)[1].astype(str) + '.' + field_goal_ball_df['playDescription'].str.split('.',expand=True)[1].str.split(' ',expand=True)[0]","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-01-05T12:41:10.313971Z","iopub.execute_input":"2022-01-05T12:41:10.314196Z","iopub.status.idle":"2022-01-05T12:41:10.371763Z","shell.execute_reply.started":"2022-01-05T12:41:10.314172Z","shell.execute_reply":"2022-01-05T12:41:10.370917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Merge with stadium data\n\ndirection_df = pd.DataFrame({'CompassDirection':['N','NE','E','SE','S','SW','W','NW','N'],\"adjusted_StadiumAzimuthAngle\":[0.0,45,90,135,180,225,270,315,360]})\ndirection_df = direction_df.sort_values(by='adjusted_StadiumAzimuthAngle')\n\nfield_goal_ball_df['adjusted_StadiumAzimuthAngle'] = field_goal_ball_df['StadiumAzimuthAngle']\nfield_goal_ball_df.loc[field_goal_ball_df.playDirection == 'right','adjusted_StadiumAzimuthAngle'] += 180\nfield_goal_ball_df.loc[field_goal_ball_df.adjusted_StadiumAzimuthAngle > 360,'adjusted_StadiumAzimuthAngle'] -= 360\n\nfield_goal_ball_df = field_goal_ball_df.sort_values(by='adjusted_StadiumAzimuthAngle')\nfield_goal_ball_df = pd.merge_asof(field_goal_ball_df,direction_df,direction='nearest')\n\nfield_goal_ball_df['StadiumDirection'] = field_goal_ball_df['StadiumName'] + '_' + field_goal_ball_df['CompassDirection']\nfield_goal_ball_df = field_goal_ball_df.sort_values(by='StadiumDirection')","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-01-05T12:41:10.373169Z","iopub.execute_input":"2022-01-05T12:41:10.373561Z","iopub.status.idle":"2022-01-05T12:41:10.414125Z","shell.execute_reply.started":"2022-01-05T12:41:10.373519Z","shell.execute_reply":"2022-01-05T12:41:10.413465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Factorize random effect\n\nplayer_idxs, players = pd.factorize(field_goal_ball_df['kickerName'])\nplayer_idxs = player_idxs.astype('int32')\n\nstadium_idxs, stadiums = pd.factorize(field_goal_ball_df['StadiumName'])\nstadium_idxs = stadium_idxs.astype('int32')\n\nstadium_direction_idxs, stadiums_direction = pd.factorize(field_goal_ball_df['StadiumDirection'])\nstadium_direction_idxs = stadium_direction_idxs.astype('int32')","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-01-05T12:41:10.415975Z","iopub.execute_input":"2022-01-05T12:41:10.417008Z","iopub.status.idle":"2022-01-05T12:41:10.426268Z","shell.execute_reply.started":"2022-01-05T12:41:10.416965Z","shell.execute_reply":"2022-01-05T12:41:10.425338Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = np.array(field_goal_ball_df['fg_make'])\nn = np.ones_like(field_goal_ball_df['fg_make'])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-05T12:41:10.428111Z","iopub.execute_input":"2022-01-05T12:41:10.428508Z","iopub.status.idle":"2022-01-05T12:41:10.438466Z","shell.execute_reply.started":"2022-01-05T12:41:10.428468Z","shell.execute_reply":"2022-01-05T12:41:10.437647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"First we use [weather data provided by Thomas Bliss](https://www.kaggle.com/tombliss/weather-data) to see if whether would affect field goal percentage","metadata":{}},{"cell_type":"code","source":"field_goal_ball_df['Temperature'] = field_goal_ball_df['Temperature'].fillna(60)\nfield_goal_ball_df['WindSpeed'] = field_goal_ball_df['WindSpeed'].fillna(7)\nfield_goal_ball_df['Precipitation'] = field_goal_ball_df['Precipitation'].fillna(0)\nfield_goal_ball_df['Pressure'] = field_goal_ball_df['Pressure'].fillna(30)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-05T12:41:10.439557Z","iopub.execute_input":"2022-01-05T12:41:10.439771Z","iopub.status.idle":"2022-01-05T12:41:10.453513Z","shell.execute_reply.started":"2022-01-05T12:41:10.439748Z","shell.execute_reply":"2022-01-05T12:41:10.452682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fg_dist = np.array(field_goal_ball_df[\"fg_dist\"])\nangle = np.array(field_goal_ball_df[\"angle\"])\npressure = np.array(field_goal_ball_df['Pressure'])\ntemperature = np.array(field_goal_ball_df['Temperature'])\nprecipitation = np.array(field_goal_ball_df['Precipitation'])\nwindspeed = np.array(field_goal_ball_df['WindSpeed'])","metadata":{"execution":{"iopub.status.busy":"2022-01-05T12:41:10.455686Z","iopub.execute_input":"2022-01-05T12:41:10.455933Z","iopub.status.idle":"2022-01-05T12:41:10.468953Z","shell.execute_reply.started":"2022-01-05T12:41:10.45588Z","shell.execute_reply":"2022-01-05T12:41:10.467981Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nwith pm.Model() as model_weather:\n\n    intercept = pm.Normal(\"intercept\", mu=0, sd=1)\n    d = pm.Normal(\"d\", mu=0, sd=10)\n    θ = pm.Normal(\"θ\", mu=0, sd=10)\n    P = pm.Normal(\"Pressure\", mu=0, sd=100)\n    T = pm.Normal(\"Temperature\", mu=0, sd=100)\n    ppt = pm.Normal(\"Precipitation\", mu=0, sd=100)\n    ws = pm.Normal(\"Windspeed\", mu=0, sd=100)\n    \n    z = intercept + pm.math.dot(fg_dist, d) + pm.math.dot(angle, θ) + pm.math.dot(pressure, P) + pm.math.dot(temperature, T) + pm.math.dot(precipitation, ppt) + pm.math.dot(windspeed, ws) \n\n    p = pm.Deterministic(\"p\", pm.math.invlogit(z))\n\n    y_obs = pm.Binomial(\"y_obs\", n=n, p=p, observed=y)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-12-30T13:34:31.136079Z","iopub.execute_input":"2021-12-30T13:34:31.137525Z","iopub.status.idle":"2021-12-30T13:34:31.236203Z","shell.execute_reply.started":"2021-12-30T13:34:31.13747Z","shell.execute_reply":"2021-12-30T13:34:31.235338Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with model_weather:\n    logit_weather_trace_jax = pm.sampling_jax.sample_numpyro_nuts(\n        2000, tune=2000, target_accept=.9,chains=2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-12-30T13:34:33.798762Z","iopub.execute_input":"2021-12-30T13:34:33.799086Z","iopub.status.idle":"2021-12-30T13:51:29.981598Z","shell.execute_reply.started":"2021-12-30T13:34:33.799053Z","shell.execute_reply":"2021-12-30T13:51:29.979395Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"az.summary(logit_weather_trace_jax)","metadata":{"execution":{"iopub.status.busy":"2021-12-30T13:51:29.983682Z","iopub.status.idle":"2021-12-30T13:51:29.984628Z","shell.execute_reply.started":"2021-12-30T13:51:29.984309Z","shell.execute_reply":"2021-12-30T13:51:29.984341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\naz.style.use(\"arviz-darkgrid\")\naz.plot_posterior(logit_weather_trace_jax, var_names=('d','θ','Pressure','Temperature','Precipitation','Windspeed'))","metadata":{"execution":{"iopub.status.busy":"2021-12-30T13:10:26.337389Z","iopub.execute_input":"2021-12-30T13:10:26.337631Z","iopub.status.idle":"2021-12-30T13:10:30.01397Z","shell.execute_reply.started":"2021-12-30T13:10:26.337603Z","shell.execute_reply":"2021-12-30T13:10:30.013269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As shown above, there are some small effect on how weather affect field goal success (e.g. field goal percentage decrease with high wind speed, high precipitation condition) but the effect is not too strong that zeros are inside the 94% credible interval for all weather effect. Thus for simplicity, the model below would not account for weather effect.\n\nNext, we would look at random effect by players and stadium, since some player may be better than others and it is more easier to covert field goal in some stadium. For stadium and player random effect:\n\n$$ p_i \\sim InverseLogit( \\theta \\cdot angle_i + d \\cdot distance_i + player_{j|i|} + stadium_{j|i|}  + intercept)$$\n\nAnd for player and stadium level:\n\n$$ player_j \\sim Normal(\\overline{player}, \\sigma_{player})$$\n$$ \\overline{player} \\sim Normal(0,10)$$\n$$ \\overline{player}  =  \\overline{player} - mean(\\overline{player})$$\n$$ \\sigma_{player} \\sim HalfCauchy(5)$$\n\n$$ stadium_j \\sim Normal(\\overline{stadium}, \\sigma_{stadium})$$\n$$ \\overline{stadium} \\sim Normal(0,10)$$\n$$ \\overline{stadium}  =  \\overline{stadium} - mean(\\overline{stadium})$$\n$$ \\sigma_{stadium} \\sim HalfCauchy(5)$$\n\nFor $\\overline{player}$ and $\\overline{stadium}$ it's zero-sumed since in a sports game you would expect the average effect of player and stadium is zero.","metadata":{}},{"cell_type":"code","source":"with pm.Model() as model_player_stadium:\n\n    intercept = pm.Normal(\"intercept\", mu=0, sd=1)\n    d = pm.Normal(\"d\", mu=0, sd=10)\n    θ = pm.Normal(\"θ\", mu=0, sd=10)\n\n\n    sigma_player =pm.HalfCauchy(\"sigma_player\", 5)\n    player_bar = pm.Normal(\"player_bar\", mu=0, sd=10)\n    player_bar = player_bar - aet.mean(player_bar)\n    \n    player = pm.Normal(\"player\", mu=player_bar, sd=sigma_player, shape=len(np.unique(player_idxs)))\n    \n    sigma_stadium =pm.HalfCauchy(\"sigma_stadium\", 5)\n    stadium_bar = pm.Normal(\"stadium_bar\", mu=0, sd=10)\n    stadium_bar = stadium_bar - aet.mean(stadium_bar)\n    \n    stadium = pm.Normal(\"stadium\", mu=stadium_bar, sd=sigma_stadium, shape=len(np.unique(stadium_idxs)))\n    \n    z = intercept + pm.math.dot(fg_dist, d) + pm.math.dot(angle, θ)  + player[player_idxs] + stadium[stadium_idxs]\n    p = pm.Deterministic(\"p\", pm.math.invlogit(z))\n\n    y_obs = pm.Binomial(\"y_obs\", n=n, p=p, observed=y)","metadata":{"execution":{"iopub.status.busy":"2021-12-28T15:24:44.809956Z","iopub.execute_input":"2021-12-28T15:24:44.810745Z","iopub.status.idle":"2021-12-28T15:24:44.905005Z","shell.execute_reply.started":"2021-12-28T15:24:44.810709Z","shell.execute_reply":"2021-12-28T15:24:44.904052Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with model_player_stadium:\n    logit_player_stadium_trace_jax = pm.sampling_jax.sample_numpyro_nuts(\n        2000, tune=2000, target_accept=.9,chains=2)","metadata":{"execution":{"iopub.status.busy":"2021-12-28T15:24:46.050247Z","iopub.execute_input":"2021-12-28T15:24:46.050529Z","iopub.status.idle":"2021-12-28T15:27:29.198436Z","shell.execute_reply.started":"2021-12-28T15:24:46.050498Z","shell.execute_reply":"2021-12-28T15:27:29.197687Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"az.plot_trace(logit_player_stadium_trace_jax, var_names=('d','θ','sigma_player','sigma_stadium'))","metadata":{"execution":{"iopub.status.busy":"2021-12-28T14:53:57.176544Z","iopub.status.idle":"2021-12-28T14:53:57.17689Z","shell.execute_reply.started":"2021-12-28T14:53:57.176702Z","shell.execute_reply":"2021-12-28T14:53:57.176727Z"},"_kg_hide-input":true,"_kg_hide-output":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kicker_df = pd.DataFrame({'Kicker':players,'coef':np.mean(logit_player_stadium_trace_jax.posterior.player,axis=(0,1)),\n                         'hdi_95':np.quantile(logit_player_stadium_trace_jax.posterior.player,axis=(0,1),q=0.95),\n                         'hdi_5':np.quantile(logit_player_stadium_trace_jax.posterior.player,axis=(0,1),q=0.05)}).sort_values(by='coef',ascending=False)","metadata":{"execution":{"iopub.status.busy":"2021-12-28T14:53:57.17823Z","iopub.status.idle":"2021-12-28T14:53:57.178542Z","shell.execute_reply.started":"2021-12-28T14:53:57.17838Z","shell.execute_reply":"2021-12-28T14:53:57.1784Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kicker_df","metadata":{"execution":{"iopub.status.busy":"2021-12-28T14:53:57.179825Z","iopub.status.idle":"2021-12-28T14:53:57.180136Z","shell.execute_reply.started":"2021-12-28T14:53:57.179975Z","shell.execute_reply":"2021-12-28T14:53:57.179996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Unsuprisingly we see Justin Tucker kicker for Baltimore Ravens on top since he is one of the best kicker in NFL history, and he's the only kicker where the credible interval do not cross zero within the span on tracking data era.","metadata":{}},{"cell_type":"code","source":"stadium_df = pd.DataFrame({'Stadium':stadiums,'coef':np.mean(logit_player_stadium_trace_jax.posterior.stadium,axis=(0,1)),\n                         'hdi_95':np.quantile(logit_player_stadium_trace_jax.posterior.stadium,axis=(0,1),q=0.95),\n                         'hdi_5':np.quantile(logit_player_stadium_trace_jax.posterior.stadium,axis=(0,1),q=0.05)}).sort_values(by='coef',ascending=False)","metadata":{"execution":{"iopub.status.busy":"2021-12-28T14:53:57.181592Z","iopub.status.idle":"2021-12-28T14:53:57.181961Z","shell.execute_reply.started":"2021-12-28T14:53:57.18176Z","shell.execute_reply":"2021-12-28T14:53:57.181797Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stadium_df","metadata":{"execution":{"iopub.status.busy":"2021-12-28T14:53:57.183077Z","iopub.status.idle":"2021-12-28T14:53:57.18339Z","shell.execute_reply.started":"2021-12-28T14:53:57.183224Z","shell.execute_reply":"2021-12-28T14:53:57.183245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We saw for example MetLife Stadium which host the New York Giants is easier to kick while Ford Field which host the Detroit Lions is harder to kick.","metadata":{}},{"cell_type":"markdown","source":"Next we would investigate how direction of stadium would affect field goal percentage. In a [blog post by NFL Operation team](https://operations.nfl.com/gameday/analytics/stats-articles/field-goal-success-probabilities-by-direction/), it studies six different statdium that has largest difference between stadium direction. The next is similar to previous one, but each stadium direction are treated as distinct stadium.","metadata":{}},{"cell_type":"code","source":"fg_dist = np.array(field_goal_ball_df[\"fg_dist\"])\nangle = np.array(field_goal_ball_df[\"angle\"])\npressure = np.array(field_goal_ball_df['Pressure'])\n\nwith pm.Model() as model_stadium_direction:\n\n    intercept = pm.Normal(\"intercept\", mu=0, sd=1)\n    d = pm.Normal(\"d\", mu=0, sd=10)\n    θ = pm.Normal(\"θ\", mu=0, sd=10)\n\n\n    sigma_player =pm.HalfCauchy(\"sigma_player\", 5)\n    player_bar = pm.Normal(\"player_bar\", mu=0, sd=10)\n    player_bar = player_bar - aet.mean(player_bar)\n    \n    player = pm.Normal(\"player\", mu=player_bar, sd=sigma_player, shape=len(np.unique(player_idxs)))\n    \n    sigma_stadium_direction =pm.HalfCauchy(\"sigma_stadium_direction\", 5)\n    stadium_direction_bar = pm.Normal(\"stadium_direction_bar\", mu=0, sd=100)\n    stadium_direction_bar = stadium_direction_bar - aet.mean(stadium_direction_bar)\n    \n    stadium_direction = pm.Normal(\"stadium_direction\", mu=stadium_direction_bar, sd=sigma_stadium_direction, shape=len(np.unique(stadium_direction_idxs)))\n    \n    z = intercept + pm.math.dot(fg_dist, d) + pm.math.dot(angle, θ) + player[player_idxs] + stadium_direction[stadium_direction_idxs]\n\n    p = pm.Deterministic(\"p\", pm.math.invlogit(z))\n \n    y_obs = pm.Binomial(\"y_obs\", n=n, p=p, observed=y)","metadata":{"execution":{"iopub.status.busy":"2021-12-28T13:38:23.715811Z","iopub.execute_input":"2021-12-28T13:38:23.716606Z","iopub.status.idle":"2021-12-28T13:38:23.812004Z","shell.execute_reply.started":"2021-12-28T13:38:23.716559Z","shell.execute_reply":"2021-12-28T13:38:23.811149Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with model_stadium_direction:\n    logit_stadium_direction_trace_jax = pm.sampling_jax.sample_numpyro_nuts(\n        2000, tune=2000, target_accept=.9,chains=2)","metadata":{"execution":{"iopub.status.busy":"2021-12-28T13:38:24.570718Z","iopub.execute_input":"2021-12-28T13:38:24.571167Z","iopub.status.idle":"2021-12-28T13:41:20.115727Z","shell.execute_reply.started":"2021-12-28T13:38:24.571128Z","shell.execute_reply":"2021-12-28T13:41:20.11489Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stadium_directions_df = pd.DataFrame({'Stadium_direction':stadiums_direction,'coef':np.mean(logit_stadium_direction_trace_jax.posterior.stadium_direction,axis=(0,1)),\n                         'hdi_95':np.quantile(logit_stadium_direction_trace_jax.posterior.stadium_direction,axis=(0,1),q=0.95),\n                         'hdi_5':np.quantile(logit_stadium_direction_trace_jax.posterior.stadium_direction,axis=(0,1),q=0.05)}).sort_values(by='coef',ascending=False)","metadata":{"execution":{"iopub.status.busy":"2021-12-28T13:41:20.119566Z","iopub.execute_input":"2021-12-28T13:41:20.120001Z","iopub.status.idle":"2021-12-28T13:41:20.144199Z","shell.execute_reply.started":"2021-12-28T13:41:20.119968Z","shell.execute_reply":"2021-12-28T13:41:20.143384Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stadium_directions_df","metadata":{"execution":{"iopub.status.busy":"2021-12-28T13:41:20.14518Z","iopub.execute_input":"2021-12-28T13:41:20.145582Z","iopub.status.idle":"2021-12-28T13:41:20.160455Z","shell.execute_reply.started":"2021-12-28T13:41:20.145547Z","shell.execute_reply":"2021-12-28T13:41:20.159763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Both direction in MetLife stadium continue to be top at converting field goal, and also interesting to see how each side of Gillette Stadium are on two different spectrum. To investigate the difference between two direction, one can just simply use the difference of two side, but also one can define a \"deviation\" term that specify how large the two direction differ, and pooled by mean of stadium. For each stadium:\n\n$$Deviation_j \\sim Normal(0,1)$$\n$$\\sigma_{Stadium direction} \\sim HalfCauchy(5)$$\n$$\\overline{Stadium Direction} = mean(Stadium Direction + Deviation, Stadium Direction - Deviation)$$\n$$Stadium Direction \\sim Normal(\\overline{Stadium Direction},\\sigma_{Stadium Direction})$$\n\n","metadata":{}},{"cell_type":"code","source":"with pm.Model() as model_stadium_pooled:\n\n    intercept = pm.Normal(\"intercept\", mu=0, sd=1)\n    d = pm.Normal(\"d\", mu=0, sd=10)\n    θ = pm.Normal(\"θ\", mu=0, sd=10)\n\n\n    sigma_player =pm.HalfCauchy(\"sigma_player\", 5)\n    player_bar = pm.Normal(\"player_bar\", mu=0, sd=10)\n    player_bar = player_bar - aet.mean(player_bar)\n    \n    player = pm.Normal(\"player\", mu=player_bar, sd=sigma_player, shape=len(np.unique(player_idxs)))\n    \n    \n    \n    \n    sigma_stadium =pm.HalfCauchy(\"sigma_stadium\", 5)\n    stadium_bar = pm.Normal(\"stadium_bar\", mu=0, sd=10)\n    stadium_bar = stadium_bar - aet.mean(stadium_bar)\n    stadium = pm.Normal(\"stadium\", mu=stadium_bar, sd=sigma_stadium, shape=len(np.unique(stadium_idxs)))\n    \n    deviation = pm.Normal(\"deviation\", mu=0, sd=1, shape=len(np.unique(stadium_idxs)))\n    stadium_direction_bar = aet.stack([stadium+deviation, stadium-deviation]).reshape((2,-1)).T.flatten()\n    sigma_stadium_direction =pm.HalfCauchy(\"sigma_stadium_direction\", 1)\n    \n    stadium_direction = pm.Normal(\"stadium_direction\", mu=stadium_direction_bar, sd=sigma_stadium_direction, shape=len(np.unique(stadium_direction_idxs)))\n    \n    z = intercept + pm.math.dot(fg_dist, d) + pm.math.dot(angle, θ)  + player[player_idxs] + stadium_direction[stadium_direction_idxs]\n    p = pm.Deterministic(\"p\", pm.math.invlogit(z))\n\n    y_obs = pm.Binomial(\"y_obs\", n=n, p=p, observed=y)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-05T12:50:07.318757Z","iopub.execute_input":"2022-01-05T12:50:07.319292Z","iopub.status.idle":"2022-01-05T12:50:12.01577Z","shell.execute_reply.started":"2022-01-05T12:50:07.319252Z","shell.execute_reply":"2022-01-05T12:50:12.013666Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with model_stadium_pooled:\n    logit_stadium_pooled_trace_jax = pm.sampling_jax.sample_numpyro_nuts(\n        2000, tune=2000, target_accept=.9,chains=2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-05T12:50:14.591869Z","iopub.execute_input":"2022-01-05T12:50:14.593137Z","iopub.status.idle":"2022-01-05T12:54:52.133289Z","shell.execute_reply.started":"2022-01-05T12:50:14.59309Z","shell.execute_reply":"2022-01-05T12:54:52.132039Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"az.style.use(\"arviz-darkgrid\")\naz.plot_trace(logit_stadium_pooled_trace_jax, var_names=('d','θ','sigma_player','sigma_stadium','deviation')) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-05T12:55:39.05956Z","iopub.execute_input":"2022-01-05T12:55:39.059871Z","iopub.status.idle":"2022-01-05T12:55:45.508996Z","shell.execute_reply.started":"2022-01-05T12:55:39.059841Z","shell.execute_reply":"2022-01-05T12:55:45.508235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stadium_directions_pooled_df = pd.DataFrame({'Stadium':stadiums,'coef':abs(np.mean(logit_stadium_pooled_trace_jax.posterior.deviation,axis=(0,1)))}).sort_values(by='coef',ascending=False)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-01-05T12:55:51.969402Z","iopub.execute_input":"2022-01-05T12:55:51.969874Z","iopub.status.idle":"2022-01-05T12:55:51.981746Z","shell.execute_reply.started":"2022-01-05T12:55:51.96983Z","shell.execute_reply":"2022-01-05T12:55:51.980182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stadium_directions_pooled_df","metadata":{"execution":{"iopub.status.busy":"2022-01-05T12:55:58.965523Z","iopub.execute_input":"2022-01-05T12:55:58.966413Z","iopub.status.idle":"2022-01-05T12:55:58.983048Z","shell.execute_reply.started":"2022-01-05T12:55:58.966366Z","shell.execute_reply":"2022-01-05T12:55:58.982071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From the model, Gillette Stadium has the highest deviation two direction of stadium","metadata":{}},{"cell_type":"markdown","source":"Lastly we compare the perforamnce of model","metadata":{}},{"cell_type":"code","source":"compare_dict = {\"Logistic Model: Player + Stadium\": logit_player_stadium_trace_jax, \n                \"Logistic Model: Player + Stadium direction\": logit_stadium_direction_trace_jax,\n                \"Logistic Model: Player + Stadium direction pooled\": logit_stadium_pooled_trace_jax\n\n               }\ndf_compare = az.compare(compare_dict, ic=\"waic\")\ndf_compare","metadata":{"execution":{"iopub.status.busy":"2021-12-28T13:55:27.049506Z","iopub.execute_input":"2021-12-28T13:55:27.049817Z","iopub.status.idle":"2021-12-28T13:55:28.742721Z","shell.execute_reply.started":"2021-12-28T13:55:27.049789Z","shell.execute_reply":"2021-12-28T13:55:28.741208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_, ax = plt.subplots(1, 1, figsize=(10, 5))\naz.plot_compare(df_compare, ax=ax);","metadata":{"execution":{"iopub.status.busy":"2021-12-28T13:55:29.847197Z","iopub.execute_input":"2021-12-28T13:55:29.847586Z","iopub.status.idle":"2021-12-28T13:55:30.136945Z","shell.execute_reply.started":"2021-12-28T13:55:29.847546Z","shell.execute_reply":"2021-12-28T13:55:30.13596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Model with stadium effect only perform the best means stadium direction is probably a small second order effect only, also since field goal distance and angle are the most important feature so three different model perform similarly as same set of features are used.","metadata":{}}]}