{"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":"The focus on the 2021 - 2022 Big Data Bowl is special teams. This demo will provide sample code to help with the process of getting started. This demo includes sample code to do the following:\n\n[Speed on Kickoff Plays](#whatis) - explains the metric of interest and how it can be used to analyze members of kickoff coverage unit.\n\n[Read Data](#ReadData) - will read the tracking / non-tracking data.\n\n[Clean Data](#CleanData) - will filter for plays of interest and filter out players.\n\n[Mixed-Effects Model](#ModelAll) - will create a mixed effects model for using player and surface.\n\n[Analyze Stability](#ModelAcross) - will analyze the stability of the modelling technique.","metadata":{}},{"cell_type":"markdown","source":"<a id=\"whatis\"></a>\n# Speed on Kickoff Plays\n\nWhen NFL kickoff specialists kick the ball to the opposing team at the beginning of a half or after a score, the other members of the kickoff coverage unit typically will sprint down the field to tackle the returner. In order to limit the return yards, the kickoff team wants to get down the field and meet ther returner as fast as possible. Thus, player max speed reached during the initial sprint (before they are blocked / obstructed by the return team) can be a valuable metric to analyze those on the kickoff team. However, taking the raw average of max speed can be decieving as other factors, such as field surface, can have an effect on player speeds as well. This demo will show how to model player speeds on kickoff coverage unit while accounting for surface.\n\n\n<img src=\"https://operations.nfl.com/images/infocus-rules/evolution-of-kickoff-poster.jpg\" style=\"width:500px;\">","metadata":{}},{"cell_type":"markdown","source":"<a id=\"ReadData\"></a>\n# Read Data\n\nBefore we can do anything, we must read in the data.","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2021-10-22T15:53:48.866964Z","iopub.execute_input":"2021-10-22T15:53:48.867793Z","iopub.status.idle":"2021-10-22T15:53:48.892236Z","shell.execute_reply.started":"2021-10-22T15:53:48.867685Z","shell.execute_reply":"2021-10-22T15:53:48.891612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"##reading in non-tracking data\n\n#includes info on games\ndf_games = pd.read_csv(\"../input/nfl-big-data-bowl-2022/games.csv\")\n\n#includes play-by-play info on specific plays\ndf_plays = pd.read_csv(\"../input/nfl-big-data-bowl-2022/plays.csv\")\n\n#includes background info for players\ndf_players = pd.read_csv(\"../input/nfl-big-data-bowl-2022/players.csv\")\n\n#Reading tracking data (2020 only for this case)\ndf_tracking = pd.read_csv(\"../input/nfl-big-data-bowl-2022/tracking2020.csv\")\n\n#includes scouting info on specific plays\ndf_PFFScouting = pd.read_csv(\"../input/nfl-big-data-bowl-2022/PFFScoutingData.csv\")\n\n#loading data from Lee Sharpe's public GitHub repository. It includes info on field surface.\ndf_leeSharpeGames = pd.read_csv(\"https://raw.githubusercontent.com/nflverse/nfldata/master/data/games.csv\")","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:53:48.893491Z","iopub.execute_input":"2021-10-22T15:53:48.894048Z","iopub.status.idle":"2021-10-22T15:54:23.565184Z","shell.execute_reply.started":"2021-10-22T15:53:48.894017Z","shell.execute_reply":"2021-10-22T15:54:23.564355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"CleanData\"></a>\n# Clean Data\n\nTo start the cleaning process, we must first use PFF scouting data to identify special teams safeties who are not actively advancing downfield. To be able to merge this information, we will use tracking data to get a map of a player's jersey number to their `nflId` in each game.","metadata":{}},{"cell_type":"code","source":"#using df_tracking to merge to jersey numbers\n\n#selecting variables of interest & dropping duplicates - jersey # is constant throughout game\ndf_jerseyMap = df_tracking.drop_duplicates(subset = [\"gameId\", \"team\", \"jerseyNumber\", \"nflId\"])\n\n#joining to games\ndf_jerseyMap = pd.merge(df_jerseyMap, df_games, left_on=['gameId'], right_on =['gameId'])\n\n#getting name of team\nconditions = [\n    (df_jerseyMap['team'] == \"home\"),\n    (df_jerseyMap['team'] != \"home\"),\n]\n\nvalues = [df_jerseyMap['homeTeamAbbr'], df_jerseyMap['visitorTeamAbbr']]\n\n#adjusting jersey number so that it includes 0 when < 10\ndf_jerseyMap['team'] = np.select(conditions, values)\n\ndf_jerseyMap['jerseyNumber'] = df_jerseyMap['jerseyNumber'].astype(str)\n\ndf_jerseyMap.loc[df_jerseyMap['jerseyNumber'].map(len) < 4, 'jerseyNumber'] = \"0\"+df_jerseyMap.loc[df_jerseyMap['jerseyNumber'].map(len) < 4, 'jerseyNumber'].str[:2]\n\ndf_jerseyMap['jerseyNumber'] = df_jerseyMap['jerseyNumber'].str[:2]\n\n#getting team and jersey\ndf_jerseyMap['teamJersey'] = df_jerseyMap['team'] + ' ' + df_jerseyMap['jerseyNumber'].str[:2]\n\n#map to merge nflId to teamJersey\ndf_jerseyMap = df_jerseyMap[['gameId', 'nflId', 'teamJersey']]\n\ndf_jerseyMap = df_jerseyMap.sort_values(['gameId', 'nflId', 'teamJersey'])\n\ndf_jerseyMap.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:54:23.566569Z","iopub.execute_input":"2021-10-22T15:54:23.566836Z","iopub.status.idle":"2021-10-22T15:54:25.203771Z","shell.execute_reply.started":"2021-10-22T15:54:23.566791Z","shell.execute_reply":"2021-10-22T15:54:25.203167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#dataframe will include gameId, playId and nflId for each special teams safety\ndf_PFF_specialTeamSafeties = df_PFFScouting.copy()\n\n#splitting into a column for each special teams safety\ndf_PFF_specialTeamSafeties[['teamJersey1', 'teamJersey2', 'teamJersey3', 'teamJersey4', 'teamJersey5', 'teamJersey6']] = df_PFF_specialTeamSafeties['specialTeamsSafeties'].str.split('; ',expand=True)\n\n#selecting jersey numbers for each team\ndf_PFF_specialTeamSafeties = df_PFF_specialTeamSafeties[['gameId', 'playId', 'teamJersey1', 'teamJersey2', 'teamJersey3', 'teamJersey4', 'teamJersey5', 'teamJersey6']]\n\n#gathering data\ndf_PFF_specialTeamSafeties = pd.melt(df_PFF_specialTeamSafeties, id_vars =['gameId', 'playId'], value_vars =['teamJersey1', 'teamJersey2', 'teamJersey3', 'teamJersey4', 'teamJersey5', 'teamJersey6'],\n               value_name = 'teamJersey')\n\n#dropping NA rows\ndf_PFF_specialTeamSafeties.dropna()\n\n#joining to jersey map\ndf_PFF_specialTeamSafeties = pd.merge(df_PFF_specialTeamSafeties, df_jerseyMap, on = ['gameId', 'teamJersey'])\n\n#selecting variables of interest\ndf_PFF_specialTeamSafeties = df_PFF_specialTeamSafeties[['gameId', 'playId', 'nflId']]\n\ndf_PFF_specialTeamSafeties = df_PFF_specialTeamSafeties.sort_values(['gameId', 'playId', 'nflId'])\n\ndf_PFF_specialTeamSafeties.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:54:25.205396Z","iopub.execute_input":"2021-10-22T15:54:25.205787Z","iopub.status.idle":"2021-10-22T15:54:25.317247Z","shell.execute_reply.started":"2021-10-22T15:54:25.205760Z","shell.execute_reply":"2021-10-22T15:54:25.316531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we will create a dataframe to only include deep kickoffs so that we can remove all unnecessary rows in the tracking data.","metadata":{}},{"cell_type":"code","source":"#creating data frame that will only include deep kickoffs\ndf_deepKickoffs = df_plays.copy()\n\n#joining the scouting data\ndf_deepKickoffs = pd.merge(df_deepKickoffs, df_PFFScouting, on = ['gameId', 'playId'])\n\n#filtering for kickoff plays only & deep kickoffs only\ndf_deepKickoffs = df_deepKickoffs[(df_deepKickoffs['specialTeamsPlayType'] == 'Kickoff') & (df_deepKickoffs['kickType'] == 'D')]\n\n#selecting variables of interest\ndf_deepKickoffs = df_deepKickoffs[['gameId', 'playId', 'kickerId', 'possessionTeam']]\n\ndf_deepKickoffs.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:54:25.318382Z","iopub.execute_input":"2021-10-22T15:54:25.318638Z","iopub.status.idle":"2021-10-22T15:54:25.386789Z","shell.execute_reply.started":"2021-10-22T15:54:25.318612Z","shell.execute_reply":"2021-10-22T15:54:25.385748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, we will use the data frames we created to filter the tracking data. We will use `df_PFF_specialTeamSafeties` to filter out players who were special teams safeties and use `df_deepKickoffs` to remove plays where there was not a deep kickoff. For each player on each play in the the tracking data, we will filter for the first 40 frames which approximately corresponds to the initial sprint portion of the play. Over that interval, we calculate the maximum speed reached for each player in the play.","metadata":{}},{"cell_type":"code","source":"df_maxSpeeds = df_tracking.copy()\n\n#joining games\ndf_maxSpeeds = pd.merge(df_maxSpeeds, df_games, on = 'gameId')\n\n#using a join to remove special teams safeties\ndf_maxSpeeds = pd.merge(left = df_maxSpeeds, right = df_PFF_specialTeamSafeties, how='left', indicator=True, on = ['gameId', 'playId', 'nflId'])\n\ndf_maxSpeeds = df_maxSpeeds.loc[df_maxSpeeds._merge == 'left_only', :].drop(columns = '_merge')\n\n#joining deep kickoffs\ndf_maxSpeeds = pd.merge(df_maxSpeeds, df_deepKickoffs, on = ['gameId', 'playId'])\n\n#removing the kicker from the tracking data\ndf_maxSpeeds = df_maxSpeeds[(df_maxSpeeds['kickerId'] != df_maxSpeeds['nflId']) &\n                            \n                            #player is on home team and kicking team is home\n                            ((df_maxSpeeds['team']=='home') & (df_maxSpeeds['possessionTeam'] == df_maxSpeeds['homeTeamAbbr']) |\n                             \n                             #or player is on away team and kicking team is away\n                            (df_maxSpeeds['team']=='away') & (df_maxSpeeds['possessionTeam'] == df_maxSpeeds['visitorTeamAbbr']))]\n\n#select variables of interest\ndf_maxSpeeds = df_maxSpeeds[['gameId', 'playId', 'frameId', 'nflId', 'event', 's']]\n\n#arranging by gameId, playId, frameId and nflId\ndf_maxSpeeds = df_maxSpeeds.sort_values(['gameId', 'playId', 'frameId'])\n\n#grouping by gameId, playId and nflId & filtering for frames after kickoff\ndf_maxSpeeds = df_maxSpeeds.loc[df_maxSpeeds.groupby(['gameId', 'playId']).event.transform(lambda z: np.cumsum(z.isin(['kickoff', 'free_kick'])) >= 1)]\n\n#grouping by gameId, playId and nflId & filtering for first 40 observations\ndf_maxSpeeds = df_maxSpeeds.groupby(['gameId', 'playId', 'nflId']).head(40).reset_index()\n\n#calculating max speed for given play / player\ndf_maxSpeeds = df_maxSpeeds.groupby(['gameId', 'playId', 'nflId']).s.apply(lambda z: z.max()).reset_index()\n\n#renaming speed column as maxSpeed\ndf_maxSpeeds = df_maxSpeeds.rename(columns={\"s\" : \"maxSpeed\"})","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:54:25.387937Z","iopub.execute_input":"2021-10-22T15:54:25.388152Z","iopub.status.idle":"2021-10-22T15:55:20.340217Z","shell.execute_reply.started":"2021-10-22T15:54:25.388126Z","shell.execute_reply":"2021-10-22T15:55:20.339610Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_maxSpeeds.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:20.341621Z","iopub.execute_input":"2021-10-22T15:55:20.343654Z","iopub.status.idle":"2021-10-22T15:55:20.355668Z","shell.execute_reply.started":"2021-10-22T15:55:20.343615Z","shell.execute_reply":"2021-10-22T15:55:20.354711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#will have max speeds for each player with surface info\ndf_maxSpeeds2 = df_maxSpeeds.copy()\n\n#merging to Lee Sharpe's data\ndf_maxSpeeds2 = pd.merge(df_maxSpeeds, df_leeSharpeGames, left_on = ['gameId'], right_on = ['old_game_id'])\n\n#selecting variables of interest\ndf_maxSpeeds2 = df_maxSpeeds2[['gameId', 'playId', 'nflId', 'week', 'surface', 'maxSpeed']]\n\n#striping the surface column to remove extra spaces at the end\ndf_maxSpeeds2['surface'] = df_maxSpeeds2['surface'].transform(lambda x : x.str.strip())\n\ndf_maxSpeeds2['surface'].unique()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:20.356879Z","iopub.execute_input":"2021-10-22T15:55:20.357132Z","iopub.status.idle":"2021-10-22T15:55:20.426468Z","shell.execute_reply.started":"2021-10-22T15:55:20.357095Z","shell.execute_reply":"2021-10-22T15:55:20.425528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"ModelAll\"></a>\n# Mixed Effects Model\n\nTo account for the lack of independence among observations from the same player, we fit a mixed effects model of player speed that includes a random intercept for each player (called a random effect) and a fixed effect surface type. The random effect represents how that player, after some shrinkage towards the overall league mean, compares to the league average. Random effects models allow us to avoid overfitting (e.g., a fixed effect for each player, and the shrinkage towards league average pulls some players more than others (those with fewer observations get pulled further towards the league average).","metadata":{}},{"cell_type":"code","source":"import statsmodels.api as sm\nimport statsmodels.formula.api as smf\n\n#fitting model\nmodel = smf.mixedlm(formula= \"maxSpeed ~ C(surface, Treatment(reference='grass'))\", data = df_maxSpeeds2, groups = 'nflId')\n\nmodelF = model.fit(method=[\"lbfgs\"])","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:20.428003Z","iopub.execute_input":"2021-10-22T15:55:20.428289Z","iopub.status.idle":"2021-10-22T15:55:24.335776Z","shell.execute_reply.started":"2021-10-22T15:55:20.428252Z","shell.execute_reply":"2021-10-22T15:55:24.334674Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's view a summary of the model:","metadata":{}},{"cell_type":"code","source":"print(modelF.summary())","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:24.339187Z","iopub.execute_input":"2021-10-22T15:55:24.339829Z","iopub.status.idle":"2021-10-22T15:55:24.678706Z","shell.execute_reply.started":"2021-10-22T15:55:24.339787Z","shell.execute_reply":"2021-10-22T15:55:24.677639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you can see in the \"Random effects\" section of the summay, the `nflId` effect has a standard deviation of ~0.5. This would imply that the model estimated the player population as a distribution with mean 0 and standard deviation of ~0.5.","metadata":{}},{"cell_type":"code","source":"modelF.params","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:24.679759Z","iopub.execute_input":"2021-10-22T15:55:24.679969Z","iopub.status.idle":"2021-10-22T15:55:24.687116Z","shell.execute_reply.started":"2021-10-22T15:55:24.679945Z","shell.execute_reply":"2021-10-22T15:55:24.686230Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Moreover, each of the surface effects appear to be signifcant. Matrix Turf stands out as the biggest effect. The intercept corresponds to an overall average.\n\nLet's analyze the player effects:","metadata":{}},{"cell_type":"code","source":"#saving results of model in data frame. extracting ids:\ndf_modelResults = pd.DataFrame.from_dict(modelF.random_effects, orient = 'index').reset_index()\n\n#renaming the columns appropriately\ndf_modelResults = df_modelResults.rename(columns={\"index\" : \"nflId\", \"nflId\" : \"effect\"})\n\n#joining by nflId\ndf_modelResults = pd.merge(df_modelResults, df_players, on = 'nflId')\n\n#selecting variables of interest\ndf_modelResults = df_modelResults[['nflId', 'displayName', 'effect']]\n\ndf_modelResults.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:24.688286Z","iopub.execute_input":"2021-10-22T15:55:24.688538Z","iopub.status.idle":"2021-10-22T15:55:25.001948Z","shell.execute_reply.started":"2021-10-22T15:55:24.688490Z","shell.execute_reply":"2021-10-22T15:55:25.001074Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our results give us a player effect for each player which we can use to rank them. First, we plot the top 25 players with the highest player effect:","metadata":{}},{"cell_type":"code","source":"#ordering for top effects\ndf_modelResults = df_modelResults.sort_values(by = 'effect', ascending = False)\n\n#filtering for top 25\ndf_visual1 = df_modelResults.head(25)","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:25.003291Z","iopub.execute_input":"2021-10-22T15:55:25.003802Z","iopub.status.idle":"2021-10-22T15:55:25.012302Z","shell.execute_reply.started":"2021-10-22T15:55:25.003758Z","shell.execute_reply":"2021-10-22T15:55:25.011408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#plotting\nimport matplotlib.pyplot as plt\n\nfig = plt.figure(figsize=(15,10))\nplt.rc('grid', linestyle=':', color='lightgray', linewidth=0.5)\nplt.grid(True, zorder = 0)\n\n#reordering to show largest effect at the top\ndf_visual1 = df_visual1.sort_values(by = 'effect')\n\n#creating bar graph\nplt.barh(list(df_visual1['displayName']), list(df_visual1['effect']), color = \"lightblue\")\nplt.xticks([0.0, 0.5, 1, 1.5])\n\n#setting labels\nplt.title(\"Top 25 Players in Max Speed Player Effect\", fontsize = 16)\nplt.xlabel(\"Max Speed Player Effect\", fontsize = 14)\nplt.ylabel(\"Player Name\", fontsize = 14)\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:25.013796Z","iopub.execute_input":"2021-10-22T15:55:25.015239Z","iopub.status.idle":"2021-10-22T15:55:25.375961Z","shell.execute_reply.started":"2021-10-22T15:55:25.015190Z","shell.execute_reply":"2021-10-22T15:55:25.375162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we analyze the density of player effects:","metadata":{}},{"cell_type":"code","source":"import seaborn as sns\n\nfig = plt.figure(figsize=(15,10))\n\n#plotting density plot\nsns.kdeplot(df_modelResults['effect'], fill = True)\n\n#setting theme\nplt.rc('grid', linestyle=':', color='lightgray', linewidth=0.5)\nplt.grid(True, zorder = 0)\n\n#setting labels\nplt.title(\"Density of Max Speed Player Effects\", fontsize = 16)\nplt.xlabel(\"Max Speed Player Effects\", fontsize = 14)\nplt.xticks([-3, -2, -1, 0, 1])\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:25.377058Z","iopub.execute_input":"2021-10-22T15:55:25.377275Z","iopub.status.idle":"2021-10-22T15:55:25.684470Z","shell.execute_reply.started":"2021-10-22T15:55:25.377252Z","shell.execute_reply":"2021-10-22T15:55:25.683643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The density curve shows that the population of player effects is somewhat different from a standard normal distribution. While most of the players are between -1 and 1, there are some outliers on the left side of curve. Although, we removed special teams safeties and kickers in the data cleaning process, some of these players likely have a non-traditional role on the kickoff team.","metadata":{}},{"cell_type":"markdown","source":"<a id=\"ModelAcross\"></a>\n# Analyze Stability\n\nTo ensure that our results are meaningful, it is important to check the stability of this analysis across weeks. If each player effect is mostly independent of the player's previous weeks' player effect, it might suggest that the analysis cannot be used to properly measure performance. Here, we will make a model using data from weeks 1-8 only and compare the results to a model using weeks 9-17 only.","metadata":{}},{"cell_type":"code","source":"#dataframe for week 1-8\ndf_wk1to8 = df_maxSpeeds2.loc[df_maxSpeeds2['week'] <= 8].copy()\n\n#fitting model filtering for week 1-8\nmodel1 = smf.mixedlm(formula= \"maxSpeed ~ C(surface, Treatment(reference='grass'))\", data = df_wk1to8, groups = 'nflId')\nmodel_wk1to8 = model1.fit(method=[\"lbfgs\"])\n\n#saving results of effects\ndf_model_wk1to8Results = pd.DataFrame.from_dict(model_wk1to8.random_effects, orient = 'index').reset_index()\ndf_model_wk1to8Results = df_model_wk1to8Results.rename(columns={\"index\" : \"nflId\", \"nflId\" : \"effect\"})\ndf_model_wk1to8Results = pd.merge(df_model_wk1to8Results, df_players, on = 'nflId')\ndf_model_wk1to8Results = df_model_wk1to8Results[['nflId', 'displayName', 'effect']]\n\ndf_model_wk1to8Results.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:25.685727Z","iopub.execute_input":"2021-10-22T15:55:25.686271Z","iopub.status.idle":"2021-10-22T15:55:27.732398Z","shell.execute_reply.started":"2021-10-22T15:55:25.686238Z","shell.execute_reply":"2021-10-22T15:55:27.731555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#dataframe for week 9-17\ndf_wk9to17 = df_maxSpeeds2.loc[df_maxSpeeds2['week'] >= 9].copy()\n\n#fitting model filtering for week 9-17\nmodel2 = smf.mixedlm(formula= \"maxSpeed ~ C(surface, Treatment(reference='grass'))\", data = df_wk9to17, groups = 'nflId')\nmodel_wk9to17 = model2.fit(method=[\"lbfgs\"])\n\n#saving results of effects\ndf_model_wk9to17Results = pd.DataFrame.from_dict(model_wk9to17.random_effects, orient = 'index').reset_index()\ndf_model_wk9to17Results = df_model_wk9to17Results.rename(columns={\"index\" : \"nflId\", \"nflId\" : \"effect\"})\ndf_model_wk9to17Results = pd.merge(df_model_wk9to17Results, df_players, on = 'nflId')\ndf_model_wk9to17Results = df_model_wk9to17Results[['nflId', 'displayName', 'effect']]\n\ndf_model_wk9to17Results.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:27.737161Z","iopub.execute_input":"2021-10-22T15:55:27.737673Z","iopub.status.idle":"2021-10-22T15:55:29.907801Z","shell.execute_reply.started":"2021-10-22T15:55:27.737628Z","shell.execute_reply":"2021-10-22T15:55:29.906918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#merging model results\ndf_model_ResultsCompare = pd.merge(df_model_wk9to17Results, df_model_wk1to8Results, on = ['nflId', 'displayName'], suffixes = ('_wk1to8', '_wk9to17'))\n\ndf_model_ResultsCompare.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T15:55:29.909151Z","iopub.execute_input":"2021-10-22T15:55:29.909660Z","iopub.status.idle":"2021-10-22T15:55:29.936001Z","shell.execute_reply.started":"2021-10-22T15:55:29.909617Z","shell.execute_reply":"2021-10-22T15:55:29.935090Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"An easy way to measure stability is to graph the values in a scatter plot. We also will use a line with a slope of 1 and intercept 0 to visualize the correlation. If the points align well with the line, we can feel good that the first sample is well correlated with the second sample.","metadata":{}},{"cell_type":"code","source":"#plotting scatter plot\nfig = plt.figure(figsize=(10,10))\nplt.rc('grid', linestyle=':', color='lightgray', linewidth=0.5)\nplt.grid(True, zorder = 0)\n\n#finding the outliers\ndf_model_ResultsCompare_highlight = df_model_ResultsCompare.loc[((abs(df_model_ResultsCompare['effect_wk1to8'] - df_model_ResultsCompare['effect_wk9to17']) > 1) |\n                                                                (df_model_ResultsCompare['effect_wk9to17'] == max(df_model_ResultsCompare['effect_wk9to17'])) |\n                                                                (df_model_ResultsCompare['effect_wk1to8'] == max(df_model_ResultsCompare['effect_wk1to8'])))]\n\n#adding points\nplt.plot(df_model_ResultsCompare['effect_wk1to8'], df_model_ResultsCompare['effect_wk9to17'], 'o', color = 'black')\n\n#labeling outliers\nfor x, y, l in zip(df_model_ResultsCompare_highlight['effect_wk1to8'],df_model_ResultsCompare_highlight['effect_wk9to17'], df_model_ResultsCompare_highlight['displayName']):\n    plt.text(x, y, l)\n\n#adding linear model trend line\nx = np.linspace(-3, 1)\ny = x\nplt.plot(x, y,  color='red')\nplt.gca().set_aspect('equal')\n\n#setting labels\nplt.title(\"Stability of Max Speed Player Effect from Weeks 1-8 to Weeks 9-17 \\n (Each dot represents a player)\", fontsize = 16)\nplt.xlabel(\"Weeks 1-8 Max Speed Player Effect\", fontsize = 14)\nplt.xticks([-3, -2, -1, 0, 1])\nplt.ylabel(\"Weeks 9-17 Max Speed Player Effect\", fontsize = 14)\nplt.yticks([-4, -2, 0])\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T16:11:22.859019Z","iopub.execute_input":"2021-10-22T16:11:22.859298Z","iopub.status.idle":"2021-10-22T16:11:23.030461Z","shell.execute_reply.started":"2021-10-22T16:11:22.859271Z","shell.execute_reply":"2021-10-22T16:11:23.029617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As suggested under the density plot, some of the players that are shown as outliers on the graph likely have a non-traditional role on the kickoff team. To analyze the trend more easily, lets remove the players who are outliers:","metadata":{}},{"cell_type":"code","source":"#plotting scatter plot\nfig = plt.figure(figsize=(10,10))\nplt.rc('grid', linestyle=':', color='lightgray', linewidth=0.5)\nplt.grid(True, zorder = 0)\n\n#filtering out outliers\ndf_model_ResultsCompareNoOutliers = df_model_ResultsCompare.loc[(df_model_ResultsCompare['effect_wk1to8'] >= -2) & \n                                        (df_model_ResultsCompare['effect_wk1to8'] <= 2) &\n                                        (df_model_ResultsCompare['effect_wk9to17'] >= -2) &\n                                        (df_model_ResultsCompare['effect_wk9to17'] <= 2)].copy()\n\n#finding the outliers\ndf_model_ResultsCompare_highlight = df_model_ResultsCompareNoOutliers.loc[((abs(df_model_ResultsCompareNoOutliers['effect_wk1to8'] - df_model_ResultsCompareNoOutliers['effect_wk9to17']) > 1) |\n                                                                (df_model_ResultsCompareNoOutliers['effect_wk9to17'] == max(df_model_ResultsCompareNoOutliers['effect_wk9to17'])) |\n                                                                (df_model_ResultsCompareNoOutliers['effect_wk1to8'] == max(df_model_ResultsCompareNoOutliers['effect_wk1to8'])))]\n\n#adding points\nplt.plot(df_model_ResultsCompareNoOutliers['effect_wk1to8'], df_model_ResultsCompareNoOutliers['effect_wk9to17'], 'o', color = 'black')\n\n#labeling outliers\nfor x, y, l in zip(df_model_ResultsCompare_highlight['effect_wk1to8'],df_model_ResultsCompare_highlight['effect_wk9to17'], df_model_ResultsCompare_highlight['displayName']):\n    plt.text(x, y, l)\n\n#adding linear model trend line\nx = np.linspace(-2, 1)\ny = x\nplt.plot(x, y,  color='red')\nplt.gca().set_aspect('equal')\n\n#setting labels\nplt.title(\"Stability of Max Speed Player Effect from Weeks 1-8 to Weeks 9-17 \\n (Each dot represents a player)\", fontsize = 16)\nplt.xlabel(\"Weeks 1-8 Max Speed Player Effect\", fontsize = 14)\nplt.xticks([-1, 0, 1])\nplt.ylabel(\"Weeks 9-17 Max Speed Player Effect\", fontsize = 14)\nplt.yticks([-1, 0, 1])\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-10-22T16:14:48.124823Z","iopub.execute_input":"2021-10-22T16:14:48.125175Z","iopub.status.idle":"2021-10-22T16:14:48.293987Z","shell.execute_reply.started":"2021-10-22T16:14:48.125135Z","shell.execute_reply":"2021-10-22T16:14:48.293119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Another way to measure stability is to calculate the $R^2$ of the player effects across weeks:","metadata":{}},{"cell_type":"code","source":"#calculating R^2 with all the data\ncorr_AllData = np.corrcoef(df_model_ResultsCompare['effect_wk1to8'], df_model_ResultsCompare['effect_wk9to17'])\ncorr_xy_AllData = corr_AllData[0,1]\nRsquared_AllData = corr_xy_AllData**2\n\n#printing R^2\nprint('R^2 Correlation from week 1-8 player effect to week 9-17 player effect:',Rsquared_AllData)\n\n#calculating R^2 after removing the outliers\ncorr_NoOutliers = np.corrcoef(df_model_ResultsCompareNoOutliers['effect_wk1to8'], df_model_ResultsCompareNoOutliers['effect_wk9to17'])\ncorr_xy_NoOutliers = corr_NoOutliers[0,1]\nRsquared_NoOutliers = corr_xy_NoOutliers**2\n\n#printing R^2\nprint('R^2 Correlation from week 1-8 player effect to week 9-17 player effect:',Rsquared_NoOutliers)","metadata":{"execution":{"iopub.status.busy":"2021-10-22T16:18:39.957258Z","iopub.execute_input":"2021-10-22T16:18:39.957582Z","iopub.status.idle":"2021-10-22T16:18:39.966804Z","shell.execute_reply.started":"2021-10-22T16:18:39.957546Z","shell.execute_reply":"2021-10-22T16:18:39.965903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Considering the modestly high $R^2$ value and the clear trend in the scatter plot, it is fair to say that this technique has reasonable stability.","metadata":{}}]}