{"cells":[{"metadata":{},"cell_type":"markdown","source":"![image](https://github.com/Lexie88rus/NFL-1st-and-Future/raw/master/assets/image.jpg)"},{"metadata":{},"cell_type":"markdown","source":"[Photo by Dave Adamson on Unsplash](https://unsplash.com/photos/-nATH0CrkMU)"},{"metadata":{},"cell_type":"markdown","source":"# NFL Injury Analysis"},{"metadata":{},"cell_type":"markdown","source":"It is incredible how we can apply data science to different areas of our life. Even in sports, data science can help to make the game safer for the athletes.\n\nThe main goal of this analysis is to find out factors, which lead to injury. In particular, find out if the effects that synthetic turf versus natural turf can have on player movements and the factors that may contribute to lower extremity injuries."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport itertools\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n        \n# image utils\nimport PIL\nfrom PIL import Image\n\n# import plotting\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport scipy.stats as ss\n\n# import machine learning\nfrom sklearn.model_selection import StratifiedKFold\nfrom imblearn.over_sampling import RandomOverSampler, SMOTE\nimport xgboost as xgb\nfrom sklearn.metrics import accuracy_score, confusion_matrix, cohen_kappa_score","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load Data"},{"metadata":{},"cell_type":"markdown","source":"Load the data as pandas dataframes:"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"play_df = pd.read_csv('../input/nfl-playing-surface-analytics/PlayList.csv')\nplayer_df = pd.read_csv('../input/nfl-playing-surface-analytics/PlayerTrackData.csv')\ninjury_df = pd.read_csv('../input/nfl-playing-surface-analytics/InjuryRecord.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data Exploration"},{"metadata":{},"cell_type":"markdown","source":"First of all, I would like to perform the initial cleaning and make visualizations to get familiar with the data and generate the ideas for further analysis."},{"metadata":{},"cell_type":"markdown","source":"### Play Dataset Exploration"},{"metadata":{},"cell_type":"markdown","source":"`1` General information:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"unique_players = play_df.PlayerKey.nunique()\nunique_games = play_df.GameID.nunique()\nunique_plays = play_df.PlayKey.nunique()\n\nprint('There are {} players in the dataset.'.format(unique_players))\nprint('There are {} games in the dataset.'.format(unique_games))\nprint('There are {} plays in the dataset.'.format(unique_plays))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`2` Game exploration:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# create a dataframe with game-level information\ngame_df = play_df[['GameID', 'StadiumType', 'FieldType', 'Weather', 'Temperature']].drop_duplicates().reset_index().drop(columns=['index'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# https://stackoverflow.com/questions/28931224/adding-value-labels-on-a-matplotlib-bar-chart\ndef add_value_labels(ax, spacing=5, decimals = 0):\n    \"\"\"Add labels to the end of each bar in a bar chart.\n\n    Arguments:\n        ax (matplotlib.axes.Axes): The matplotlib object containing the axes\n            of the plot to annotate.\n        spacing (int): The distance between the labels and the bars.\n    \"\"\"\n\n    # For each bar: Place a label\n    for rect in ax.patches:\n        # Get X and Y placement of label from rect.\n        y_value = rect.get_height()\n        x_value = rect.get_x() + rect.get_width() / 2\n\n        # Number of points between bar and label. Change to your liking.\n        space = spacing\n        # Vertical alignment for positive values\n        va = 'bottom'\n\n        # If value of bar is negative: Place label below bar\n        if y_value < 0:\n            # Invert space to place label below\n            space *= -1\n            # Vertically align label at top\n            va = 'top'\n\n        # Use Y value as label and format number with one decimal place\n        format_str = \"{:.\" + str(decimals) + \"f}\"\n        label = format_str.format(y_value)\n\n        # Create annotation\n        ax.annotate(\n            label,                      # Use `label` as label\n            (x_value, y_value),         # Place label at end of the bar\n            xytext=(0, space),          # Vertically shift label by `space`\n            textcoords=\"offset points\", # Interpret `xytext` as offset in points\n            ha='center',                # Horizontally center label\n            va=va)                      # Vertically align label differently for\n                                        # positive and negative values.\n\n\ndef visualize_game_features(game_df, rotation = 90, add_labels = False, figsize=(10,10)):\n    plt.style.use('ggplot')\n    fig = plt.figure(figsize=figsize)\n    grid = plt.GridSpec(4, 3, hspace=0.2, wspace=0.2)\n    stadium_ax = fig.add_subplot(grid[0, :2])\n    fieldtype_ax = fig.add_subplot(grid[0, 2])\n    weather_ax = fig.add_subplot(grid[1, 0:])\n    temperature_ax = fig.add_subplot(grid[2, 0:])\n    temperature_box_ax = fig.add_subplot(grid[3, 0:])\n\n    stadium_ax.bar(game_df.StadiumType.value_counts().keys(), game_df.StadiumType.value_counts().values, color='#00c2c7')\n    stadium_ax.set_title('StadiumType')\n    stadium_ax.set_xticklabels(game_df.StadiumType.value_counts().keys(), rotation=rotation)\n    \n    if add_labels:\n        add_value_labels(stadium_ax, spacing=5)\n\n    fieldtype_ax.bar(game_df.FieldType.value_counts().keys(), game_df.FieldType.value_counts().values, color=['#00c2c7', '#ff9e15'])\n    fieldtype_ax.set_title('FieldType')\n    fieldtype_ax.set_xticklabels(game_df.FieldType.value_counts().keys(), rotation=0)\n    \n    if add_labels:\n        add_value_labels(fieldtype_ax, spacing=5)\n\n    weather_ax.bar(game_df.Weather.value_counts().keys(), game_df.Weather.value_counts().values, color='#00c2c7')\n    weather_ax.set_title('Weather')\n    weather_ax.set_xticklabels(game_df.Weather.value_counts().keys(), rotation=rotation)\n    \n    if add_labels:\n        add_value_labels(weather_ax, spacing=5)\n        \n    temperature_ax.hist(game_df.Temperature.astype(int).values, bins=30, range=(0,90))\n    temperature_ax.set_xlim(0,110)\n    temperature_ax.set_xticks(range(0,110,10))\n    temperature_ax.set_xticklabels(range(0,110,10))\n    temperature_ax.set_title('Temperature')\n    \n    temperature_box_ax.boxplot(game_df.Temperature.astype(int).values, vert=False)\n    temperature_box_ax.set_xlim(0,110)\n    temperature_box_ax.set_xticks(range(0,110,10))\n    temperature_box_ax.set_xticklabels(range(0,110,10))\n    temperature_box_ax.set_yticklabels(['Temperature'])\n\n    plt.suptitle('Game-Level Exploration', fontsize=16)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def clean_weather(row):\n    cloudy = ['Cloudy 50% change of rain', 'Hazy', 'Cloudy.', 'Overcast', 'Mostly Cloudy',\n          'Cloudy, fog started developing in 2nd quarter', 'Partly Cloudy',\n          'Mostly cloudy', 'Rain Chance 40%',' Partly cloudy', 'Party Cloudy',\n          'Rain likely, temps in low 40s', 'Partly Clouidy', 'Cloudy, 50% change of rain','Mostly Coudy', '10% Chance of Rain',\n          'Cloudy, chance of rain', '30% Chance of Rain', 'Cloudy, light snow accumulating 1-3\"',\n          'cloudy', 'Coudy', 'Cloudy with periods of rain, thunder possible. Winds shifting to WNW, 10-20 mph.',\n         'Cloudy fog started developing in 2nd quarter', 'Cloudy light snow accumulating 1-3\"',\n         'Cloudywith periods of rain, thunder possible. Winds shifting to WNW, 10-20 mph.',\n         'Cloudy 50% change of rain', 'Cloudy and cold',\n       'Cloudy and Cool', 'Partly cloudy']\n    \n    clear = ['Clear, Windy',' Clear to Cloudy', 'Clear, highs to upper 80s',\n             'Clear and clear','Partly sunny',\n             'Clear, Windy', 'Clear skies', 'Sunny', 'Partly Sunny', 'Mostly Sunny', 'Clear Skies',\n             'Sunny Skies', 'Partly clear', 'Fair', 'Sunny, highs to upper 80s', 'Sun & clouds', 'Mostly sunny','Sunny, Windy',\n             'Mostly Sunny Skies', 'Clear and Sunny', 'Clear and sunny','Clear to Partly Cloudy', 'Clear Skies',\n            'Clear and cold', 'Clear and warm', 'Clear and Cool', 'Sunny and cold', 'Sunny and warm', 'Sunny and clear']\n    \n    rainy = ['Rainy', 'Scattered Showers', 'Showers', 'Cloudy Rain', 'Light Rain', 'Rain shower', 'Rain likely, temps in low 40s.', 'Cloudy, Rain']\n    \n    snow = ['Heavy lake effect snow']\n    \n    indoor = ['Controlled Climate', 'Indoors', 'N/A Indoor', 'N/A (Indoors)']\n        \n    if row.Weather in cloudy:\n        return 'Cloudy'\n    \n    if row.Weather in indoor:\n        return 'Indoor'\n    \n    if row.Weather in clear:\n        return 'Clear'\n    \n    if row.Weather in rainy:\n        return 'Rain'\n    \n    if row.Weather in snow:\n        return 'Snow'\n      \n    if row.Weather in ['Cloudy.', 'Heat Index 95', 'Cold']:\n        return np.nan\n    \n    return row.Weather\n\ndef clean_stadiumtype(row):\n    if row.StadiumType in ['Bowl', 'Heinz Field', 'Cloudy']:\n        return np.nan\n    else:\n        return row.StadiumType\n\ndef clean_play_df(play_df):\n    play_df_cleaned = play_df.copy()\n    \n    # clean StadiumType\n    play_df_cleaned['StadiumType'] = play_df_cleaned['StadiumType'].str.replace(r'Oudoor|Outdoors|Ourdoor|Outddors|Outdor|Outside', 'Outdoor')\n    play_df_cleaned['StadiumType'] = play_df_cleaned['StadiumType'].str.replace(r'Indoors|Indoor, Roof Closed|Indoor, Open Roof', 'Indoor')\n    play_df_cleaned['StadiumType'] = play_df_cleaned['StadiumType'].str.replace(r'Closed Dome|Domed, closed|Domed, Open|Domed, open|Dome, closed|Domed', 'Dome')\n    play_df_cleaned['StadiumType'] = play_df_cleaned['StadiumType'].str.replace(r'Retr. Roof-Closed|Outdoor Retr Roof-Open|Retr. Roof - Closed|Retr. Roof-Open|Retr. Roof - Open|Retr. Roof Closed', 'Retractable Roof')\n    play_df_cleaned['StadiumType'] = play_df_cleaned['StadiumType'].str.replace('Open', 'Outdoor')\n    play_df_cleaned['StadiumType'] = play_df_cleaned.apply(lambda row: clean_stadiumtype(row), axis=1)\n    \n    # clean Weather\n    play_df_cleaned['Weather'] = play_df_cleaned.apply(lambda row: clean_weather(row), axis=1)\n    \n    return play_df_cleaned","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"play_df_cleaned = clean_play_df(play_df)\ngame_df_cleaned = play_df_cleaned[['GameID', 'StadiumType', 'FieldType', 'Weather', 'Temperature']].drop_duplicates().reset_index().drop(columns=['index'])\nvisualize_game_features(game_df_cleaned, rotation=0, add_labels = True, figsize=(12,16))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The diagrams above show:\n\n* __The majority of the games are played outdoors.__ That's why the weather should be important.\n* __Synthetic turf is very common, but there are around 30% more games with natural turf.__\n* __The temperature and other weather conditions vary greately__. It will be interesting to look at how the temperature affects the injuries."},{"metadata":{},"cell_type":"markdown","source":"`3` Player level exploration:"},{"metadata":{"trusted":true},"cell_type":"code","source":"player_data_df = play_df_cleaned[['PlayerKey', 'RosterPosition', 'PlayerGamePlay', 'Position', 'PositionGroup']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def visualize_player_features(player_df, figsize=(25,20), add_labels=False):\n    plt.style.use('ggplot')\n    fig = plt.figure(figsize=figsize)\n    \n    grid = plt.GridSpec(3, 4, hspace=0.2, wspace=0.2)\n    \n    plays_ax = fig.add_subplot(grid[0, 0:2])\n    max_rolling_plays_ax = fig.add_subplot(grid[1, 0:2])\n    \n    rosterposition_ax = fig.add_subplot(grid[0, 2:])\n    positiongroup_ax = fig.add_subplot(grid[1, 2:])\n    position_ax = fig.add_subplot(grid[2, 0:])\n    \n    plays_ax.hist(player_df.groupby(by=['PlayerKey']).count()['RosterPosition'].values, bins=20, color='#00c2c7')\n    plays_ax.set_title('Number of plays per player')\n    \n    max_rolling_plays_ax.hist(player_df.groupby(by=['PlayerKey']).PlayerGamePlay.max().values, bins=20, color='#00c2c7')\n    max_rolling_plays_ax.set_title('Maximum number of rolling plays per player')\n    \n    rosterposition_ax.bar(player_df.RosterPosition.value_counts().keys().values, player_df.RosterPosition.value_counts().values)\n    rosterposition_ax.set_xticklabels(player_df.RosterPosition.value_counts().keys().values, rotation=20)\n    rosterposition_ax.set_title('Roster Position')\n    if add_labels:\n        add_value_labels(rosterposition_ax, spacing=5)\n    \n    position_ax.bar(player_df.Position.value_counts().keys().values, player_df.Position.value_counts().values, color='#ff9e15')\n    position_ax.set_title('Position')\n    if add_labels:\n        add_value_labels(position_ax, spacing=5)\n        \n    positiongroup_ax.bar(player_df.PositionGroup.value_counts().keys().values, player_df.PositionGroup.value_counts().values)\n    positiongroup_ax.set_title('Position Group')\n    if add_labels:\n        add_value_labels(positiongroup_ax, spacing=5)\n    \n    plt.suptitle('Player-Level Exploration', fontsize=16)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_player_features(player_data_df, add_labels=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The diagrams above show that the circumstances around the players vary greatly: some players do much more games than the others, there is a lot of different playing position. All of these factors can contribute to injuries. "},{"metadata":{},"cell_type":"markdown","source":"`4` Play-level exploration:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def visualize_play(play_df_cleaned):\n    plt.style.use('ggplot')\n    fig, ax = plt.subplots(1,1,figsize=(15,5))\n    \n    plt.bar(play_df_cleaned.PlayType.value_counts().keys().values, play_df_cleaned.PlayType.value_counts().values)\n    plt.xticks(range(len(play_df_cleaned.PlayType.value_counts().keys().values)), play_df_cleaned.PlayType.value_counts().keys().values, rotation=20)\n    add_value_labels(ax, spacing=5)\n    plt.title('Play-Level Exploration: PlayType', fontsize=16)\n    \n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_play(play_df_cleaned)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To my mind, it will be exciting to see how play types affect the injury. I guess that some play types can cause much more dangerous situations than others!"},{"metadata":{},"cell_type":"markdown","source":"### Player Dataset Visualization and Exploration"},{"metadata":{},"cell_type":"markdown","source":"`1` Visualize players positions and paths:\n\nTo my mind, this is the most exciting part: we can take actual data and track the player's position!"},{"metadata":{},"cell_type":"markdown","source":"Let's start with visualizing just one line of the dataset:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def get_position(player_df, playkey, time):\n    '''\n    Helper function to get the player position at a certain time\n    '''\n    row = player_df[(player_df['PlayKey']==playkey) & (player_df['time']==time)]\n    event = row['event'].values[0]\n    x = row['x'].values[0]\n    y = row['y'].values[0]\n    direction = row['dir'].values[0]\n    distance = row['dis'].values[0]\n    orientation = row['o'].values[0]\n    speed = row['s'].values[0]\n    \n    return event, x, y, direction, distance, orientation, speed\n    \ndef visualize_player_position(player_df, playkey, time, figsize=(24,10)):\n    '''\n    Helper function to visualize a single player position.\n    '''\n    \n    event, x, y, direction, distance, orientation, speed = get_position(player_df, playkey, time)\n    \n    background = Image.open('../input/nfl-utils/nfl_coordinates.png')\n    background = background.resize((1200,533))\n    background = background.transpose(PIL.Image.FLIP_TOP_BOTTOM)\n    \n    # show background\n    fig = plt.figure(figsize=figsize)\n    plt.imshow(np.array(background).transpose(0,1,2), origin='lower')\n    \n    # show orientation\n    dx = 5\n    dy = dx*np.tan(np.radians(90 + orientation))\n    plt.arrow(x*10,y*10,dx,dy, color = '#767676', width=5)\n    plt.plot(x*10,y*10,color = '#767676', label='orientation')\n    \n    # show speed\n    dx = speed*20\n    dy = dx*np.tan(np.radians(90 + direction))\n    plt.arrow(x*10,y*10,dx,dy, color = '#004c97', width=5)\n    plt.plot(x*10,y*10, color = '#004c97', label='speed')\n    \n    # show current player position\n    plt.scatter(x*10,y*10,s=200,color='#e01e5a',marker='x')\n    plt.annotate('({x:.1f},{y:.1f})'.format(x=x,y=y), (x*10, y*10), xytext=(x*10,y*10-30), color='#e01e5a')\n    \n    plt.xticks(range(0,1200,100), range(0,120,10))\n    plt.yticks(range(0,533,100), range(0,53,10))\n    \n    plt.title('{play}:{time} {event}'.format(play=playkey, time=time, event = event))\n    \n    plt.legend()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_player_position(player_df, '26624-1-1', 0.0, figsize=(10,5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_player_position(player_df, '47888-13-55', 35.9, figsize=(10,5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can also plot the entire path of the player during the play and see how he moved:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def visualize_player_track(player_df, playkey, figsize=(24,10), max_positions = 5):\n    '''\n    Helper function to visualize the whole path of the player\n    '''\n    # get timestamps per player\n    timestamps = player_df[player_df['PlayKey'] == playkey].time.unique()\n    \n    background = Image.open('../input/nfl-utils/nfl_coordinates.png')\n    background = background.resize((1200,533))\n    background = background.transpose(PIL.Image.FLIP_TOP_BOTTOM)\n    \n    # show background\n    fig = plt.figure(figsize=figsize)\n    plt.imshow(np.array(background).transpose(0,1,2), origin='lower')\n    \n    # iterate through timestamps\n    positions_x, positions_y = [], []\n    for i in range(0, len(timestamps), len(timestamps) // max_positions):\n        time = timestamps[i]\n        # get position and orientation\n        event, x, y, direction, distance, orientation, speed = get_position(player_df, playkey, time)\n        positions_x.append(x*10)\n        positions_y.append(y*10)\n        \n        # plot orientation and speed for the lasi point\n        if (len(timestamps) - i < len(timestamps) // max_positions):\n            # show orientation\n            dx = 5\n            dy = dx*np.tan(np.radians(90 + orientation))\n            plt.arrow(x*10,y*10,dx,dy, color = '#767676', width=5)\n            plt.plot(x*10,y*10,color = '#767676', label='orientation')\n\n            # show speed\n            dx = speed*20\n            dy = dx*np.tan(np.radians(90 + direction))\n            plt.arrow(x*10,y*10,dx,dy, color = '#004c97', width=5)\n            plt.plot(x*10,y*10, color = '#004c97', label='speed')\n            \n            # add x-marker and annotation only for the last position\n            plt.scatter(x*10,y*10,s=200,color='#e01e5a',marker='x')\n            plt.annotate('({x:.1f},{y:.1f})'.format(x=x,y=y), (x*10, y*10), xytext=(x*10,y*10-30), color='#e01e5a')\n        \n    # show player's path\n    plt.scatter(positions_x,positions_y,s=50,color='#e01e5a',marker='o')\n    plt.plot(positions_x,positions_y,color='#e01e5a', label='player path', linestyle='--')\n    \n    # creat custom tick labels\n    plt.xticks(range(0,1200,100), range(0,120,10))\n    plt.yticks(range(0,533,100), range(0,53,10))\n    \n    # add title\n    plt.title('{play}'.format(play=playkey))\n    \n    plt.legend()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_player_track(player_df, '26624-1-1', max_positions = 5, figsize=(10,5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_player_track(player_df, '47888-13-55', max_positions = 5, figsize=(10,5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`2` Create a heatmap of the field"},{"metadata":{},"cell_type":"markdown","source":"Let's look at the heatmap of the field:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def visualize_field_heatmap(player_df, xbins=13, ybins=6, annotate = False):\n    # create a grid\n    x = np.linspace(0, 120, xbins)\n    y = np.linspace(0, 53, ybins)\n    \n    # initialize heatmap\n    hmap = np.zeros((xbins,ybins))\n    \n    for i in range(xbins-1):\n        for j in range(ybins-1):\n            hmap[i,j] = len(player_df[(player_df.x >= x[i]) & (player_df.x <= x[i+1]) & (player_df.y >= y[j]) & (player_df.y <= y[j+1])])\n            \n    fig = plt.figure(figsize=(10,5))\n    ax = sns.heatmap(np.transpose(hmap), annot=annotate, fmt = '.0f')\n    plt.title('Field Heatmap \\n the most visited areas of the field are highlighted')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_field_heatmap(player_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that some areas of the field are more 'busy' than the others."},{"metadata":{},"cell_type":"markdown","source":"### Injuries Dataset EDA"},{"metadata":{},"cell_type":"markdown","source":"Let's look quickly at the injuries:"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('There are {} injury records in total.'.format(len(injury_df)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('{} unique players injured'.format(len(injury_df.PlayerKey.unique())))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This means that there are players injured twice!"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('There are {} PlayKey values missing.'.format(len(injury_df) - injury_df.PlayKey.count()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For around 30% of the injuries, we don't know the play, when the injury took place. It is quite a significant number."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def visualize_injury(injury_df):\n    injury_df_cleaned = injury_df.copy()\n    injury_df_cleaned.DM_M1 = injury_df_cleaned.DM_M1 - injury_df_cleaned.DM_M7\n    injury_df_cleaned.DM_M7 = injury_df_cleaned.DM_M7 - injury_df_cleaned.DM_M28\n    injury_df_cleaned.DM_M28 = injury_df_cleaned.DM_M28 - injury_df_cleaned.DM_M42\n    \n    fig, axs = plt.subplots(1,3, figsize=(15,5))\n    \n    axs[0].bar(injury_df_cleaned.BodyPart.value_counts().keys().values, injury_df_cleaned.BodyPart.value_counts().values, color='#00c2c7')\n    axs[0].set_title('Body Part')\n    add_value_labels(axs[0], spacing=5)\n    \n    axs[1].bar(injury_df_cleaned.Surface.value_counts().keys().values, injury_df_cleaned.Surface.value_counts().values, color='#ff9e15')\n    axs[1].set_title('Surface')\n    add_value_labels(axs[1], spacing=5)\n    \n    M1 = injury_df_cleaned.DM_M1.sum()\n    M7 = injury_df_cleaned.DM_M7.sum()\n    M28 = injury_df_cleaned.DM_M28.sum()\n    M42 = injury_df_cleaned.DM_M42.sum()\n    \n    axs[2].bar(['1-7', '7-28', '28-42', '>=42'], [M1, M7, M28, M42])\n    axs[2].set_title('Missed Days')\n    add_value_labels(axs[2], spacing=5)\n    \n    plt.suptitle('Injury', fontsize=16)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_injury(injury_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"All in all, we have just **105** injury records! So few!\n\n* We can see that the __knees and ankles are the most frequently injured body parts.__\n* We have almost __equal number of observations for synthetic and natural surface.__"},{"metadata":{},"cell_type":"markdown","source":"## Injury Analysis"},{"metadata":{},"cell_type":"markdown","source":"`1` Merge injury data with game data (drop the play-level information), so we can explore the dependencies between the injuries and game-level factors."},{"metadata":{"trusted":true},"cell_type":"code","source":"# joined cleaned games dataset and injury dataset\ngame_injury_df = injury_df.set_index('GameID').join(game_df_cleaned.set_index('GameID'), how = 'outer')\n\n# fill null values for the injury columns with zeros\ngame_injury_df['DM_M1'] = game_injury_df['DM_M1'].fillna(0).astype(int)\ngame_injury_df['DM_M7'] = game_injury_df['DM_M7'].fillna(0).astype(int)\ngame_injury_df['DM_M28'] = game_injury_df['DM_M28'].fillna(0).astype(int)\ngame_injury_df['DM_M42'] = game_injury_df['DM_M42'].fillna(0).astype(int)\n\ngame_injury_df.DM_M1 = game_injury_df.DM_M1 - game_injury_df.DM_M7\ngame_injury_df.DM_M7 = game_injury_df.DM_M7 - game_injury_df.DM_M28\ngame_injury_df.DM_M28 = game_injury_df.DM_M28 - game_injury_df.DM_M42\n\n# introduce a column with a flag indicating an injury\ngame_injury_df['Injury'] = game_injury_df['DM_M1'] + game_injury_df['DM_M7'] + game_injury_df['DM_M28'] + game_injury_df['DM_M42']\n\n# drop duplicated surface column\ngame_injury_df = game_injury_df.drop(columns=['Surface'])\n\n# drop play-level features just for now\ngame_injury_df = game_injury_df.drop(columns=['PlayerKey', 'PlayKey'])\n\n# create dummy variables\ngame_injury_df_dummies = pd.get_dummies(game_injury_df, dummy_na = True, drop_first = True).drop(columns=['FieldType_nan'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`2` Now let's visualize the correlations:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"corr_df = game_injury_df_dummies[['Temperature', 'StadiumType_Indoor', 'StadiumType_Outdoor', 'StadiumType_Retractable Roof',\n                                 'FieldType_Synthetic', 'Weather_Cloudy', 'Weather_Rain', 'Weather_Snow', 'Injury']].corr()\n\nfig = plt.figure(figsize=(10,7))\nsns.heatmap(corr_df, annot=True, cmap=sns.diverging_palette(220, 20, as_cmap=True))\nplt.title('Correlation Heatmap')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Unfortunately, there is __no correlation between our features and the injury.__"},{"metadata":{},"cell_type":"markdown","source":"Let's look at [Cramer's V](https://en.wikipedia.org/wiki/Cram%C3%A9r%27s_V):"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Source:\n# https://stackoverflow.com/questions/46498455/categorical-features-correlation/46498792#46498792\ndef cramers_v(confusion_matrix):\n    \"\"\" calculate Cramers V statistic for categorial-categorial association.\n        uses correction from Bergsma and Wicher,\n        Journal of the Korean Statistical Society 42 (2013): 323-328\n    \"\"\"\n    chi2 = ss.chi2_contingency(confusion_matrix)[0]\n    n = confusion_matrix.sum()\n    phi2 = chi2 / n\n    r, k = confusion_matrix.shape\n    phi2corr = max(0, phi2 - ((k-1)*(r-1))/(n-1))\n    rcorr = r - ((r-1)**2)/(n-1)\n    kcorr = k - ((k-1)**2)/(n-1)\n    return np.sqrt(phi2corr / min((kcorr-1), (rcorr-1)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Source:\n# https://stackoverflow.com/questions/51859894/how-to-plot-a-cramer-s-v-heatmap-for-categorical-features\ndef plot_cramers_v_heatmap(df, cols):\n    corrM = np.zeros((len(cols),len(cols)))\n\n    for col1, col2 in itertools.combinations(cols, 2):\n        idx1, idx2 = cols.index(col1), cols.index(col2)\n        corrM[idx1, idx2] = cramers_v(pd.crosstab(df[col1], df[col2]).as_matrix())\n        corrM[idx2, idx1] = corrM[idx1, idx2]\n\n    corr = pd.DataFrame(corrM, index=cols, columns=cols)\n    fig, ax = plt.subplots(figsize=(7, 6))\n    ax = sns.heatmap(corr, annot=True, cmap=sns.diverging_palette(220, 20, as_cmap=True), ax=ax); ax.set_title(\"Cramer V Correlation between Variables\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = [\"StadiumType\", \"FieldType\", \"Weather\", \"Temperature\", \"Injury\"]\nplot_cramers_v_heatmap(game_injury_df, cols)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This table looks way more exciting!\n* We can see that there are correlations between Weather and Temperature (obviously, snow and rain means lower temperature than sunny weather).\n* There is also some correlation between Stadium Type and Temperature.\n\nBut unfortunately, still no correlation with the injury."},{"metadata":{},"cell_type":"markdown","source":"Let's try [Theil’s U](https://docs.oracle.com/cd/E40248_01/epm.1112/cb_statistical/frameset.htm?ch07s02s03s04.html) to find some insights:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Source:\n# https://github.com/shakedzy/dython/blob/master/dython/nominal.py\n\nfrom collections import Counter\nimport math\n\ndef conditional_entropy(x, y, nan_strategy='replace', nan_replace_value=0):\n    \"\"\"\n    Calculates the conditional entropy of x given y: S(x|y)\n    Wikipedia: https://en.wikipedia.org/wiki/Conditional_entropy\n    **Returns:** float\n    Parameters\n    ----------\n    x : list / NumPy ndarray / Pandas Series\n        A sequence of measurements\n    y : list / NumPy ndarray / Pandas Series\n        A sequence of measurements\n    nan_strategy : string, default = 'replace'\n        How to handle missing values: can be either 'drop' to remove samples with missing values, or 'replace'\n        to replace all missing values with the nan_replace_value. Missing values are None and np.nan.\n    nan_replace_value : any, default = 0.0\n        The value used to replace missing values with. Only applicable when nan_strategy is set to 'replace'.\n    \"\"\"\n    y_counter = Counter(y)\n    xy_counter = Counter(list(zip(x,y)))\n    total_occurrences = sum(y_counter.values())\n    entropy = 0.0\n    for xy in xy_counter.keys():\n        p_xy = xy_counter[xy] / total_occurrences\n        p_y = y_counter[xy[1]] / total_occurrences\n        entropy += p_xy * math.log(p_y/p_xy)\n    return entropy\n\ndef theils_u(x, y):\n    s_xy = conditional_entropy(x,y)\n    x_counter = Counter(x)\n    total_occurrences = sum(x_counter.values())\n    p_x = list(map(lambda n: n/total_occurrences, x_counter.values()))\n    s_x = ss.entropy(p_x)\n    if s_x == 0:\n        return 1\n    else:\n        return (s_x - s_xy) / s_x","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def plot_theils_u_heatmap(df, cols):\n    corrM = np.zeros((len(cols),len(cols)))\n\n    for col1, col2 in itertools.combinations(cols, 2):\n        idx1, idx2 = cols.index(col1), cols.index(col2)\n        corrM[idx1, idx2] = theils_u(df[col1].values, df[col2].values)\n        corrM[idx2, idx1] = corrM[idx1, idx2]\n\n    corr = pd.DataFrame(corrM, index=cols, columns=cols)\n    fig, ax = plt.subplots(figsize=(7, 6))\n    ax = sns.heatmap(corr, annot=True, cmap=sns.diverging_palette(220, 20, as_cmap=True), ax=ax); ax.set_title(\"Theil's U Correlation between Variables\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = [\"StadiumType\", \"FieldType\", \"Weather\", \"Temperature\", \"Injury\"]\nplot_theils_u_heatmap(game_injury_df, cols)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Again, no insights on  the injury."},{"metadata":{},"cell_type":"markdown","source":"`3` Explore the game features with respect to the injury:"},{"metadata":{},"cell_type":"markdown","source":"Temperature:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# get the temperature values for games \n# be sure to only take the values where the game was played outdoors\nnon_injury_temp = game_injury_df[(game_injury_df.StadiumType == 'Outdoor') & (game_injury_df.Temperature >= 0) & (game_injury_df.Injury == 0)].Temperature.values\ninjury_temp = game_injury_df[(game_injury_df.StadiumType == 'Outdoor') & (game_injury_df.Temperature >= 0) & (game_injury_df.Injury == 1)].Temperature.values","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10,5))\nplt.boxplot([non_injury_temp, injury_temp], vert = False)\nplt.title('Temperature Distribution (Box Plot)')\nplt.yticks([1,2], ['No Injury', 'Injury'])\nplt.xlim(0,100)\nplt.xlabel('Temperature')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(1,1, figsize=(10,5))\nsns.kdeplot(non_injury_temp, label = 'No Injury')\nsns.kdeplot(injury_temp, label = 'Injury')\nplt.title('Temperature Distribution')\nplt.xlabel('Temperature')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Check the Kolmogorov-Smirnov test to check if the temperature samples come from one distribution:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare maximum play speed samples\nD, pvalue = ss.ks_2samp(non_injury_temp, injury_temp)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The p-value is high. We can't say that temperature values come from different distributions."},{"metadata":{},"cell_type":"markdown","source":"`4` Hypothesis testing"},{"metadata":{},"cell_type":"markdown","source":"Let's try a simple hypothesis test to check if the number of injuries on synthetic turf is significantly higher. We see that basically, the probability of injury is higher when playing on synthetic turf. But is this statistically significant?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compute the injury probability for synthetic and natural turf\np_injury = game_injury_df[['FieldType', 'Injury']].groupby('FieldType').mean()['Injury']\np_injury","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The null hypothesis is that the probability of the injury on natural turf is less or equal to the probability of the injury on synthetic turf."},{"metadata":{"trusted":true},"cell_type":"code","source":"# get number of trials and overall injury rate under null\nn_natural = game_injury_df[['FieldType', 'Injury']].groupby('FieldType').size()[0]\nn_synthetic = game_injury_df[['FieldType', 'Injury']].groupby('FieldType').size()[1]\n\np_null = game_injury_df[['FieldType', 'Injury']][game_injury_df.FieldType == 'Natural'].mean()[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# compute standard error\nse_p = np.sqrt(p_null * (1-p_null) * (1/n_natural + 1/n_synthetic))\n\n#  compute z-score and p-value\nz = (p_injury[1] - p_injury[0]) / se_p\n\nprint('The z-score is: {}'.format(z))\nprint('The p-value is: {}'.format(1-ss.norm.cdf(z)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Under 95% confidence, we reject the null hypothesis if the p-value is less than 0.05. Our computed p-value is approximately 0.002, so we can reject the null hypothesis that the probability of injury on natural turf is less or equal than the injury probability on the synthetic turf. We can conclude that the probability of the injury on synthetic turf is higher than the probability of the injury on natural turf (with 99% confidence)."},{"metadata":{},"cell_type":"markdown","source":"`5` Play-level features:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# joined cleaned games dataset and injury dataset\nplay_injury_df = injury_df.dropna(subset=['PlayKey']).set_index('PlayKey').join(play_df_cleaned.set_index('PlayKey'), how = 'outer', lsuffix='_left', rsuffix='_right')\n\n# fill null values for the injury columns with zeros\nplay_injury_df['DM_M1'] = play_injury_df['DM_M1'].fillna(0).astype(int)\nplay_injury_df['DM_M7'] = play_injury_df['DM_M7'].fillna(0).astype(int)\nplay_injury_df['DM_M28'] = play_injury_df['DM_M28'].fillna(0).astype(int)\nplay_injury_df['DM_M42'] = play_injury_df['DM_M42'].fillna(0).astype(int)\n\n# introduce a column with a flag indicating an injury\nplay_injury_df.DM_M1 = play_injury_df.DM_M1 - play_injury_df.DM_M7\nplay_injury_df.DM_M7 = play_injury_df.DM_M7 - play_injury_df.DM_M28\nplay_injury_df.DM_M28 = play_injury_df.DM_M28 - play_injury_df.DM_M42\n\nplay_injury_df['Injury'] = play_injury_df['DM_M1'] + play_injury_df['DM_M7'] + play_injury_df['DM_M28']+ play_injury_df['DM_M42']\n\n# drop duplicated surface column\nplay_injury_df = play_injury_df.drop(columns=['Surface'])\n\n# create dummy variables\nplay_injury_df_dummies = pd.get_dummies(play_injury_df, columns = ['PlayType', 'PositionGroup'], dummy_na = True, drop_first = True)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"corr_df = play_injury_df_dummies[['PlayType_Pass', 'PlayType_Kickoff', 'PlayType_Punt', 'PlayType_Rush',\n                                  'PositionGroup_QB', 'PositionGroup_DL', 'PositionGroup_LB', 'PositionGroup_OL',\n                                  'PositionGroup_RB', 'PositionGroup_SPEC', 'PositionGroup_TE', 'PositionGroup_WR',\n                                  'Injury']].corr()\n\nfig = plt.figure(figsize=(15,7))\nsns.heatmap(corr_df, annot=True, cmap=sns.diverging_palette(220, 20, as_cmap=True))\nplt.title('Correlation Heatmap')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Cramer's V:"},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = [\"RosterPosition\", \"StadiumType\", \"FieldType\", \"Temperature\", \"Weather\", 'PlayType', 'PlayerGamePlay', 'Position', 'PositionGroup', 'Injury']\nplot_cramers_v_heatmap(play_injury_df, cols)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Again, this chart shows us some obvious correlations, but still, nothing is correlated with the injury."},{"metadata":{},"cell_type":"markdown","source":"Let's tru Theil's U:"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_theils_u_heatmap(play_injury_df, cols)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Unfortunately, there are no insights on injury."},{"metadata":{},"cell_type":"markdown","source":"Visualize the heatmap of the field with the plays with the injury:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# find the keys of the plays related to injuries\nplay_injuries = play_injury_df.reset_index().dropna()[['PlayKey']]\n# merge the playkeys associated with injuries with players' positions\nplayer_injuries = player_df.merge(play_injuries, on='PlayKey', how='inner') # use inner join!\n# visualize the heatmap\nvisualize_field_heatmap(player_injuries, annotate = True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It looks like the heatmap of injuries is somewhat different from the general heatmap of the field! Probably, there are more dangerous areas on the field, where players are more likely to get injured!\n\nWe can use this insight for the feature engineering for the injury prediction model!"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def visualize_field_kde(player_df):\n    # plot kde\n    fig = plt.figure(figsize=(10,5))\n    ax = sns.jointplot(x=\"x\", y=\"y\", data=player_df, kind=\"kde\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is the KDE plot for the injury play locations on the field:"},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_field_kde(player_injuries)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Game Timeline and Injury Analysis"},{"metadata":{},"cell_type":"markdown","source":"Let's plot a timeline of player's games and highlight the days with injuries:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def player_games_timeline(player_key, play_df, injury_df):\n    '''\n    Function to plot the player's timeline\n    '''\n    player_games = play_df[play_df.PlayerKey == player_key][['GameID', 'PlayKey', 'PlayerDay', 'PlayerGame']]\n    \n    # plot timeline for the player\n    plt.figure(figsize=(20,5))\n    plt.title('Player Games Timeline \\n PlayerKey: ' + str(player_key))  \n    plt.plot(player_games.PlayerDay.unique(), np.zeros(len(player_games.PlayerDay.unique())), color='#00c2c7')\n    plt.scatter(player_games.PlayerDay.unique(), np.zeros(len(player_games.PlayerDay.unique())), s=100, color='#00c2c7', label='games')\n    \n    # add games with injury\n    injured_players = injury_df.PlayerKey.unique()\n    if player_key in injured_players:\n        injury_games = injury_df[injury_df.PlayerKey == player_key].GameID.values\n        injury_days = player_games[player_games.GameID.isin(injury_games)].PlayerDay.unique()\n        \n        plt.scatter(injury_days, np.zeros(len(injury_days)), s=100, color='#e01e5a', label='injury')\n    \n    plt.legend()\n    plt.xlabel('days')\n    plt.yticks([])\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"player_games_timeline(26624, play_df, injury_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"player_games_timeline(33337, play_df, injury_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"player_games_timeline(43540, play_df, injury_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"player_games_timeline(39873, play_df, injury_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I picked some random injured players from the dataset. It seems that injuries happen at the beginning of the season (for each player)."},{"metadata":{},"cell_type":"markdown","source":"So is the speed of the player, his position on the field and the number of the game (`PlayerGame`) correlated with the injury?"},{"metadata":{"trusted":true},"cell_type":"code","source":"play_injury = play_injury_df[['PlayerGame', 'PlayerGamePlay', 'Injury']]\ncorrs = play_injury.corr()\n\nfig = plt.figure(figsize=(7,5))\nsns.heatmap(corrs, annot=True, cmap=sns.diverging_palette(220, 20, as_cmap=True))\nplt.title('Correlation Heatmap')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The diagram above shows that there is no correlation between the injury and the number of games/plays per game played."},{"metadata":{},"cell_type":"markdown","source":"## Data Engineering"},{"metadata":{},"cell_type":"markdown","source":"My next step is to try to extract features and try to build a model to predict the injury. Hopefully, I will be able to obtain the feature importances out of this model and find out the injury drivers."},{"metadata":{"trusted":true},"cell_type":"code","source":"# reset index \nfeatures_df = play_injury_df.copy().reset_index()\n# drop the irrelevant columns\nfeatures_df = features_df.drop(columns=['PlayerKey_left','GameID_left','BodyPart', 'PlayKey', 'PlayerKey_right', 'GameID_right', 'DM_M1', 'DM_M7', 'DM_M28', 'DM_M42'])\n# convert into dummies\nfeatures_df = pd.get_dummies(features_df, dummy_na = False, drop_first = True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Machine Learning"},{"metadata":{},"cell_type":"markdown","source":"We have the prepared data, now we can move to the machine learning part."},{"metadata":{},"cell_type":"markdown","source":"`1` Split into features and targets:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# split into X and y\ny = features_df['Injury']\nX = features_df.drop(columns=['Injury'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`2` Use stratified split to split into train and test set so that class labels will be evenly distributed among the train and the test set:"},{"metadata":{"trusted":true},"cell_type":"code","source":"skf = StratifiedKFold(n_splits=2)\n\nfor train_index, test_index in skf.split(X, y):\n    X_train, X_test = X.values[train_index, :], X.values[test_index, :]\n    y_train, y_test = y[train_index], y[test_index]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`3` Resampling the dataset:\n\nThe number of plays with injury is so tiny that we can't just use the data as-is. Our model will only always predict no injury. That is why we need to resample our dataset. Resampling will help us to generate more examples of cases from the rare class (the injury), so we can train our model! I will use [imbalanced-learn](https://imbalanced-learn.readthedocs.io/en/stable/index.html) library."},{"metadata":{"trusted":true},"cell_type":"code","source":"res = RandomOverSampler(random_state=0)\nX_resampled, y_resampled = res.fit_resample(X_train, y_train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`4` Training the model:"},{"metadata":{"trusted":true},"cell_type":"code","source":"model = xgb.XGBClassifier(max_depth=3,\n                      learning_rate=0.1,\n                      n_estimators=100,\n                      objective='binary:logistic',\n                      booster='gbtree',\n                      tree_method='auto',\n                      n_jobs=50,\n                      gamma=0,\n                      min_child_weight=1,\n                      max_delta_step=0,\n                      subsample=1,\n                      colsample_bytree=1,\n                      colsample_bylevel=1,\n                      colsample_bynode=1,\n                      reg_alpha=0,\n                      reg_lambda=1,\n                      scale_pos_weight=1,\n                      base_score=0.5,\n                      random_state=42)\nmodel.fit(X_resampled, y_resampled)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`5` Evaluate the model:"},{"metadata":{"trusted":true},"cell_type":"code","source":"y_pred = model.predict(X_test)\n\naccuracy = accuracy_score(y_test, y_pred)\nconf_matrix = confusion_matrix(y_test, y_pred)\ncohen_kappa = cohen_kappa_score(y_test, y_pred)\n\nprint('Accuracy: {}'.format(accuracy))\nprint('Cohen kappa: {}'.format(cohen_kappa))\nprint('Confusion Matrix: \\n {}'.format(conf_matrix))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The model's performance is quite poor, but we can still look at the feature importances:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# extract the feature importances from XGBoost model\nfeature_importances = model.feature_importances_\nfeature_importances = pd.DataFrame(feature_importances, index = X.columns).reset_index().rename(columns={'index':'feature', 0:'importance'}).sort_values(by=['importance'], ascending=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot as a bar chart\nplt.figure(figsize=(20,7))\nplt.bar(range(len(feature_importances)), feature_importances.importance.values)\nplt.xticks(range(len(feature_importances)), feature_importances.feature.values, rotation=90)\n\nplt.title('Feature importances')\nplt.xlabel('features')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looking at the chart above, we can say:\n* Play type is among features impacting injury;\n* Field type feature is also among the important features in our injury predicting model;\n* Player game and plyer day are also among the important features impacting the injury."},{"metadata":{},"cell_type":"markdown","source":"## Motion Data Analysis"},{"metadata":{},"cell_type":"markdown","source":"One of the most exciting things in this competition is the analysis of player movements. We can extract a lot of different features for each play, for example:\n* the average speed per play (might be calculated out of the coordinates or just the average per play);\n* the sharpest turn per play;\n* the average angle between the direction and the orientation;\n* the maximum angle between the direction and the orientation."},{"metadata":{},"cell_type":"markdown","source":"The dataset is too big to process on Kaggle, so I will include my preprocessing code here (see `create_motion_data_df`) and upload the result in the dataset."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def create_motion_data_df(injury_df, play_df, player_df):\n    '''\n    Function to create aggregated values for motion data.\n    '''\n    # create a field with the angle between direction and orientation\n    player_df['angle'] = player_df['o'] - player_df['dir']\n    \n    # greate grouped dataframes\n    grouped_max = player_df[['PlayKey','time', 'dir', 'dis', 'o', 's', 'angle']].groupby(by=['PlayKey']).max()\n    grouped_average = player_df[['PlayKey','time', 'dir', 'dis', 'o', 's', 'angle']].groupby(by=['PlayKey']).mean()\n    \n    # merge with play data\n    play_df = play_df.merge(grouped_max.reset_index(), on=['PlayKey'])\n    play_df = play_df.merge(grouped_average.reset_index(), on=['PlayKey'], suffixes=('_max', '_avg'))\n    \n    # merge with injury data\n    injury_df = injury_df.drop(columns=['PlayerKey', 'GameID', 'BodyPart', 'Surface']).merge(play_df, on=['PlayKey'], how='outer').fillna(0)\n    \n    return injury_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load the preprocessed data\nmotion_df = pd.read_csv('../input/nfl-utils/motion_data.csv').drop(columns=['Unnamed: 0'])\nmotion_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# introduce a column with a flag indicating an injury\nmotion_df.DM_M1 = motion_df.DM_M1 - motion_df.DM_M7\nmotion_df.DM_M7 = motion_df.DM_M7 - motion_df.DM_M28\nmotion_df.DM_M28 = motion_df.DM_M28 - motion_df.DM_M42\nmotion_df['Injury'] = motion_df['DM_M1'] + motion_df['DM_M7'] + motion_df['DM_M28']+ motion_df['DM_M42']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Explore the correlations:"},{"metadata":{"trusted":true},"cell_type":"code","source":"corr_df = motion_df[['RosterPosition', 'PlayerDay', 'PlayerGame', 'StadiumType', 'FieldType',\n       'Temperature', 'Weather', 'PlayType', 'PlayerGamePlay', 'Position',\n       'PositionGroup', 'time_max', 'dir_max', 'dis_max', 'o_max', 's_max',\n       'angle_max', 'time_avg', 'dir_avg', 'dis_avg', 'o_avg', 's_avg',\n       'angle_avg', 'Injury']].corr()\n\nfig = plt.figure(figsize=(20,20))\nsns.heatmap(corr_df, annot=True, cmap=sns.diverging_palette(220, 20, as_cmap=True))\nplt.title('Correlation Heatmap')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the distances and speed per play have a negative correlation with PlayerGamePlay. Does it mean that players get tired and move slower?"},{"metadata":{},"cell_type":"markdown","source":"### Explore Motion for Plays with Injury vs Plays with No Injury"},{"metadata":{},"cell_type":"markdown","source":"Let's look at the distribution of motion features for plays with injury:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(3,2,figsize=(15,10))\n\naxs[0,0].hist(motion_df.s_max, bins = 20)\naxs[0,0].set_title('All: Maximum Speed Distribution')\n\naxs[1,0].hist(motion_df[motion_df.Injury == 1].s_max, bins = 20)\naxs[1,0].set_title('Injury: Maximum Speed Distribution')\n\naxs[0,1].hist(motion_df.s_avg, bins = 20)\naxs[0,1].set_title('All: Average Speed Distribution')\n\naxs[1,1].hist(motion_df[motion_df.Injury == 1].s_avg, bins = 20)\naxs[1,1].set_title('Injury: Average Speed Distribution')\n\nsns.kdeplot(motion_df.s_max, ax = axs[2,0], label='all')\nsns.kdeplot(motion_df[motion_df.Injury == 1].s_max,ax = axs[2,0], label='injury')\naxs[2,0].set_title('Maximum Speed Distribution')\n\nsns.kdeplot(motion_df.s_avg, ax = axs[2,1], label='all')\nsns.kdeplot(motion_df[motion_df.Injury == 1].s_avg,ax = axs[2,1], label='injury')\naxs[2,1].set_title('Average Speed Distribution')\n\nplt.suptitle('Speed Distribution', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems that speed for plays with injuries may come from a different distribution. Let's try [Kolmogorov–Smirnov test](https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test) to check if those samples come from one distribution:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare maximum play speed samples\nD, pvalue = ss.ks_2samp(motion_df.s_max, motion_df[motion_df.Injury == 1].s_max)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare average play speed samples\nD, pvalue = ss.ks_2samp(motion_df.s_avg, motion_df[motion_df.Injury == 1].s_avg)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The p-value is very small, so we can reject the null hypothesis that both samples are drawn from one distribution. It means that most certainly, we do have some difference between speed for normal plays and plays with the injury."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(3,2,figsize=(15,10))\n\naxs[0,0].hist(motion_df.dis_max, bins = 20)\naxs[0,0].set_title('All: Maximum Distance Distribution')\n\naxs[1,0].hist(motion_df[motion_df.Injury == 1].dis_max, bins = 20)\naxs[1,0].set_title('Injury: Maximum Distance Distribution')\n\naxs[0,1].hist(motion_df.dis_avg, bins = 20)\naxs[0,1].set_title('All: Average Distance Distribution')\n\naxs[1,1].hist(motion_df[motion_df.Injury == 1].dis_avg, bins = 20)\naxs[1,1].set_title('Injury: Average Distance Distribution')\n\nsns.kdeplot(motion_df.dis_max, ax = axs[2,0], label='all')\nsns.kdeplot(motion_df[motion_df.Injury == 1].dis_max,ax = axs[2,0], label='injury')\naxs[2,0].set_title('Maximum Distance Distribution')\n\nsns.kdeplot(motion_df.dis_avg, ax = axs[2,1], label='all')\nsns.kdeplot(motion_df[motion_df.Injury == 1].dis_avg,ax = axs[2,1], label='injury')\naxs[2,1].set_title('Average Distance Distribution')\n\nplt.suptitle('Distance Distribution', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Apply Kolmogorov-Smirnov test:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare maximum play speed samples\nD, pvalue = ss.ks_2samp(motion_df.dis_max, motion_df[motion_df.Injury == 1].dis_max)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare maximum play speed samples\nD, pvalue = ss.ks_2samp(motion_df.dis_avg, motion_df[motion_df.Injury == 1].dis_avg)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, axs = plt.subplots(3,2,figsize=(15,10))\n\naxs[0,0].hist(motion_df.angle_max, bins = 20)\naxs[0,0].set_title('All: Maximum Angle Distribution')\n\naxs[1,0].hist(motion_df[motion_df.Injury == 1].angle_max, bins = 20)\naxs[1,0].set_title('Injury: Maximum Angle Distribution')\n\naxs[0,1].hist(motion_df.angle_avg, bins = 20)\naxs[0,1].set_title('All: Average Angle Distribution')\n\naxs[1,1].hist(motion_df[motion_df.Injury == 1].angle_avg, bins = 20)\naxs[1,1].set_title('Injury: Average Angle Distribution')\n\nsns.kdeplot(motion_df.angle_max, ax = axs[2,0], label='all')\nsns.kdeplot(motion_df[motion_df.Injury == 1].angle_max,ax = axs[2,0], label='injury')\naxs[2,0].set_title('Maximum Angle Distribution')\n\nsns.kdeplot(motion_df.angle_avg, ax = axs[2,1], label='all')\nsns.kdeplot(motion_df[motion_df.Injury == 1].angle_avg,ax = axs[2,1], label='injury')\naxs[2,1].set_title('Average Angle Distribution')\n\nplt.suptitle('Angle between Direction and Orientation Distribution', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Kolmogorov-Smirnov test:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare maximum play speed samples\nD, pvalue = ss.ks_2samp(motion_df.angle_max, motion_df[motion_df.Injury == 1].angle_max)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare maximum play speed samples\nD, pvalue = ss.ks_2samp(motion_df.angle_avg, motion_df[motion_df.Injury == 1].angle_avg)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see a little different result here: average angle distribution for all plays and plays with injury are close. We can see that the p-value for these distributions is higher, and the KDE plots look somewhat similar."},{"metadata":{},"cell_type":"markdown","source":"In general, the motion features for plays related to injury look different: their distributions have two spikes. One of the peaks is around zero value. Probably, the injury made the player stop for a while. That is why we see speed and distance values concentrated around zero."},{"metadata":{},"cell_type":"markdown","source":"### Explore Motion for Plays on Synthetic vs Natural Turf"},{"metadata":{},"cell_type":"markdown","source":"Explore the speed:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, axs = plt.subplots(3,2,figsize=(15,10))\n\naxs[0,0].hist(motion_df[motion_df.FieldType == 'Natural'].s_max, bins = 20)\naxs[0,0].set_title('Natural Turf: Maximum Speed Distribution')\n\naxs[1,0].hist(motion_df[motion_df.FieldType == 'Synthetic'].s_max, bins = 20)\naxs[1,0].set_title('Synthetic Turf: Maximum Speed Distribution')\n\naxs[0,1].hist(motion_df[motion_df.FieldType =='Natural'].s_avg, bins = 20)\naxs[0,1].set_title('Natural Turf: Average Speed Distribution')\n\naxs[1,1].hist(motion_df[motion_df.FieldType == 'Synthetic'].s_avg, bins = 20)\naxs[1,1].set_title('Synthetic Turf: Average Speed Distribution')\n\nsns.kdeplot(motion_df[motion_df.FieldType == 'Natural'].s_max, ax = axs[2,0], label='natural')\nsns.kdeplot(motion_df[motion_df.FieldType == 'Synthetic'].s_max,ax = axs[2,0], label='synthetic')\naxs[2,0].set_title('Maximum Speed Distribution')\n\nsns.kdeplot(motion_df[motion_df.FieldType == 'Natural'].s_avg, ax = axs[2,1], label='natural')\nsns.kdeplot(motion_df[motion_df.FieldType == 'Synthetic'].s_avg,ax = axs[2,1], label='synthetic')\naxs[2,1].set_title('Average Speed Distribution')\n\nplt.suptitle('Speed Distribution', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the speed distribution for the natural and synthetic field type are almost identical!"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, axs = plt.subplots(3,2,figsize=(15,10))\n\naxs[0,0].hist(motion_df[motion_df.FieldType == 'Natural'].dis_max, bins = 20)\naxs[0,0].set_title('Natural Turf: Maximum Distance Distribution')\n\naxs[1,0].hist(motion_df[motion_df.FieldType == 'Synthetic'].dis_max, bins = 20)\naxs[1,0].set_title('Synthetic Turf: Maximum Distance Distribution')\n\naxs[0,1].hist(motion_df[motion_df.FieldType =='Natural'].dis_avg, bins = 20)\naxs[0,1].set_title('Natural Turf: Average Distance Distribution')\n\naxs[1,1].hist(motion_df[motion_df.FieldType == 'Synthetic'].dis_avg, bins = 20)\naxs[1,1].set_title('Synthetic Turf: Average Distance Distribution')\n\nsns.kdeplot(motion_df[motion_df.FieldType == 'Natural'].dis_max, ax = axs[2,0], label='natural')\nsns.kdeplot(motion_df[motion_df.FieldType == 'Synthetic'].dis_max,ax = axs[2,0], label='synthetic')\naxs[2,0].set_title('Maximum Distance Distribution')\n\nsns.kdeplot(motion_df[motion_df.FieldType == 'Natural'].dis_avg, ax = axs[2,1], label='natural')\nsns.kdeplot(motion_df[motion_df.FieldType == 'Synthetic'].dis_avg,ax = axs[2,1], label='synthetic')\naxs[2,1].set_title('Average Distance Distribution')\n\nplt.suptitle('Distance Distribution', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The same thing is for the distance."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, axs = plt.subplots(3,2,figsize=(15,10))\n\naxs[0,0].hist(motion_df[motion_df.FieldType == 'Natural'].angle_max, bins = 20)\naxs[0,0].set_title('Natural Turf: Maximum Angle Distribution')\n\naxs[1,0].hist(motion_df[motion_df.FieldType == 'Synthetic'].angle_max, bins = 20)\naxs[1,0].set_title('Synthetic Turf: Maximum Angle Distribution')\n\naxs[0,1].hist(motion_df[motion_df.FieldType =='Natural'].angle_avg, bins = 20)\naxs[0,1].set_title('Natural Turf: Average Angle Distribution')\n\naxs[1,1].hist(motion_df[motion_df.FieldType == 'Synthetic'].angle_avg, bins = 20)\naxs[1,1].set_title('Synthetic Turf: Average Angle Distribution')\n\nsns.kdeplot(motion_df[motion_df.FieldType == 'Natural'].angle_max, ax = axs[2,0], label='natural')\nsns.kdeplot(motion_df[motion_df.FieldType == 'Synthetic'].angle_max,ax = axs[2,0], label='synthetic')\naxs[2,0].set_title('Maximum Angle Distribution')\n\nsns.kdeplot(motion_df[motion_df.FieldType == 'Natural'].angle_avg, ax = axs[2,1], label='natural')\nsns.kdeplot(motion_df[motion_df.FieldType == 'Synthetic'].angle_avg,ax = axs[2,1], label='synthetic')\naxs[2,1].set_title('Average Angle Distribution')\n\nplt.suptitle('Angle between Direction and Orientation Distribution', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We get identical distributions for the angle too.\n\nIt means that the type of turf does not affect the main motion features."},{"metadata":{},"cell_type":"markdown","source":"We didn't see any difference between the motion features on different playing surfaces in general, but maybe some differences are depending on the weather conditions or the player position?"},{"metadata":{},"cell_type":"markdown","source":"`1` Movement under different weather conditions:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# clean the weather and stadium type values in the motion dataset\nmotion_df_cleaned = clean_play_df(motion_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def compare_distributions_weather(df, field, title=''):\n    '''\n    Helper function to compare the distribution of the specified field\n    '''\n    fig, axs = plt.subplots(4,1,figsize=(10,12))\n    sns.kdeplot(df[(df.Weather == 'Cloudy') & (df.FieldType == 'Natural') & (df.StadiumType == 'Outdoor')][field].values, ax = axs[0], label='Natural')\n    sns.kdeplot(df[(df.Weather == 'Cloudy') & (df.FieldType == 'Synthetic') & (df.StadiumType == 'Outdoor')][field].values, ax = axs[0], label='Synthetic')\n    axs[0].set_title('Cloudy Weather')\n    \n    sns.kdeplot(df[(df.Weather == 'Clear') & (df.FieldType == 'Natural') & (df.StadiumType == 'Outdoor')][field].values, ax = axs[1], label='Natural')\n    sns.kdeplot(df[(df.Weather == 'Clear') & (df.FieldType == 'Synthetic') & (df.StadiumType == 'Outdoor')][field].values, ax = axs[1], label='Synthetic')\n    axs[1].set_title('Clear Weather')\n    \n    sns.kdeplot(df[(df.Weather == 'Rain') & (df.FieldType == 'Natural') & (df.StadiumType == 'Outdoor')][field].values, ax = axs[2], label='Natural')\n    sns.kdeplot(df[(df.Weather == 'Rain') & (df.FieldType == 'Synthetic') & (df.StadiumType == 'Outdoor')][field].values, ax = axs[2], label='Synthetic')\n    axs[2].set_title('Rainy Weather')\n    \n    sns.kdeplot(df[(df.Weather == 'Snow') & (df.FieldType == 'Natural') & (df.StadiumType == 'Outdoor')][field].values, ax = axs[3], label='Natural')\n    sns.kdeplot(df[(df.Weather == 'Snow') & (df.FieldType == 'Synthetic') & (df.StadiumType == 'Outdoor')][field].values, ax = axs[3], label='Synthetic')\n    axs[3].set_title('Snowy Weather')\n    \n    plt.suptitle(title, fontsize=16)\n    \n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_distributions_weather(motion_df_cleaned, 's_max', 'Maximum Speed Distribution')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_distributions_weather(motion_df_cleaned, 's_avg', 'Average Speed Distribution')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"__Observations:__\n* We have a slight difference for maximum speed distributions for rainy and snowy weather (when playing on the outdoor stadium, of course).  We see that the maximum speed on synthetic turf is a bit lower.\n* We can see that the average speed for the rainy and snowy weather is lower on the synthetic turf."},{"metadata":{},"cell_type":"markdown","source":"Kolmogorov-Smirnov test:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare the samples for rainy weather\nsynth_rain = motion_df_cleaned[(motion_df_cleaned.Weather == 'Rain') & (motion_df_cleaned.FieldType == 'Synthetic') & (motion_df_cleaned.StadiumType == 'Outdoor')]['s_avg'].values\nnat_rain = motion_df_cleaned[(motion_df_cleaned.Weather == 'Rain') & (motion_df_cleaned.FieldType == 'Natural') & (motion_df_cleaned.StadiumType == 'Outdoor')]['s_avg'].values\n\nD, pvalue = ss.ks_2samp(synth_rain, nat_rain)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare the samples for snowy weather\nsynth_snow = motion_df_cleaned[(motion_df_cleaned.Weather == 'Snow') & (motion_df_cleaned.FieldType == 'Synthetic') & (motion_df_cleaned.StadiumType == 'Outdoor')]['s_avg'].values\nnat_snow = motion_df_cleaned[(motion_df_cleaned.Weather == 'Snow') & (motion_df_cleaned.FieldType == 'Natural') & (motion_df_cleaned.StadiumType == 'Outdoor')]['s_avg'].values\n\nD, pvalue = ss.ks_2samp(synth_snow, nat_snow)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The p-value is very small for both snowy and rainy weather. So with a certain level of confidence, we can conclude that maximum speed values for snowy and rainy weather come from different distributions. __Probably, the synthetic turf under severe weather conditions makes the players slow down.__"},{"metadata":{},"cell_type":"markdown","source":"`2` Movement depending on a player position:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def compare_distributions_position(df, field, title=''):\n    '''\n    Helper function to compare the distribution of a specific variable for different player positions\n    '''\n    # get the positions\n    positions = df.Position.unique()\n    \n    # creat subplots\n    fig, axs = plt.subplots(len(positions),1,figsize=(7,25))\n    \n    for i, pos in enumerate(positions):\n        sns.kdeplot(df[(df.Position == pos) & (df.FieldType == 'Natural')][field].values, ax = axs[i], label='Natural')\n        sns.kdeplot(df[(df.Position == pos) & (df.FieldType == 'Synthetic')][field].values, ax = axs[i], label='Synthetic')\n        axs[i].set_title(pos)\n    \n    plt.suptitle(title, fontsize=16)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_distributions_position(motion_df_cleaned, 's_max', title='Average Speed Distribution Comparison')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_distributions_position(motion_df_cleaned, 's_avg', title='Maximum Speed Distribution Comparison')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"__Observations:__\nThere are some slight differences in speed distribution for DB (defensive back) and HB (halfback) player positions. Surprisingly, in both cases, the speed is more distributed towards higher values.\n\nLet's conduct a test:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare the samples for DB\nsynth_DB = motion_df_cleaned[(motion_df_cleaned.Position == 'DB') & (motion_df_cleaned.FieldType == 'Synthetic')]['s_avg'].values\nnat_DB = motion_df_cleaned[(motion_df_cleaned.Position == 'DB') & (motion_df_cleaned.FieldType == 'Natural')]['s_avg'].values\n\nD, pvalue = ss.ks_2samp(synth_DB, nat_DB)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# compare the samples for HB\nsynth_HB = motion_df_cleaned[(motion_df_cleaned.Position == 'HB') & (motion_df_cleaned.FieldType == 'Synthetic')]['s_avg'].values\nnat_HB = motion_df_cleaned[(motion_df_cleaned.Position == 'HB') & (motion_df_cleaned.FieldType == 'Natural')]['s_avg'].values\n\nD, pvalue = ss.ks_2samp(synth_HB, nat_HB)\nD, pvalue","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The Kolmogorov-Smirnov test shows that the distribution of the average speed for the DB (defensive back) position is different for synthetic and natural turf. Probably, the players move slightly faster on the synthetic turf."},{"metadata":{},"cell_type":"markdown","source":"## Synthetic vs Natural Turf\n\nHypothesis testing proved that there is a higher risk of injury for the synthetic turf. But why? Let's try to look at different factors, such as weather, temperature, movement features.\n\n`1` Synthetic and natural turf and the weather conditions:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# group by FieldType and StaduimType\ngrouped_df = game_injury_df.groupby(by=['FieldType', 'StadiumType']).count().reset_index()\n\n# get the number of games on natural turf for each stadium type\nnatural_gr = grouped_df[grouped_df.FieldType == 'Natural'][['StadiumType','Injury']]\nnatural_gr = natural_gr.append({'StadiumType':'Indoor',  'Injury' : 0}, ignore_index = True).sort_values(by=['StadiumType'])\n\n# get the number of games on synthetic turf for each stadium type\nsynthetic_gr = grouped_df[grouped_df.FieldType == 'Synthetic'][['StadiumType','Injury']].sort_values(by=['StadiumType'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# visualize the stadium type for natural and synthetic turf\nn_groups = len(natural_gr)\nfig, ax = plt.subplots(figsize=(7,7))\nindex = np.arange(n_groups)\nbar_width = 0.35\nopacity = 0.8\n\nrects1 = plt.bar(index, natural_gr.Injury.values, bar_width,\nalpha=opacity,\ncolor='b',\nlabel='Natural Turf')\n\nrects2 = plt.bar(index + bar_width, synthetic_gr.Injury.values, bar_width,\nalpha=opacity,\ncolor='g',\nlabel='Synthetic Turf')\n\nplt.xticks(index + bar_width/2, natural_gr.StadiumType.values)\nplt.xlabel('Stadium Type')\nplt.ylabel('Number of Games')\nplt.legend()\nplt.title('Stadium Type for Synthetic/Natural turf', fontsize=16)\nadd_value_labels(ax)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Most of the games are played outdoors. We see that there are a lot of outdoor games played on synthetic turf.\n\nCan extreme weather conditions increase the risk of injury?"},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_injury_nat = game_injury_df[(game_injury_df['StadiumType'] == 'Outdoor') & (game_injury_df['FieldType'] == 'Natural')]\\\n.groupby(by=['Weather']).sum().reset_index().sort_values(by=['Weather'])\n\nweather_injury_synth = game_injury_df[(game_injury_df['StadiumType'] == 'Outdoor') & (game_injury_df['FieldType'] == 'Synthetic')]\\\n.groupby(by=['Weather']).sum().reset_index().sort_values(by=['Weather'])","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# visualize the staduin type for natural and synthetic turf\nn_groups = len(weather_injury_synth)\nfig, ax = plt.subplots(figsize=(7,7))\nindex = np.arange(n_groups)\nbar_width = 0.35\nopacity = 0.8\n\nrects1 = plt.bar(index, weather_injury_nat.Injury.values, bar_width,\nalpha=opacity,\ncolor='b',\nlabel='Natural Turf')\n\nrects2 = plt.bar(index + bar_width, weather_injury_synth.Injury.values, bar_width,\nalpha=opacity,\ncolor='g',\nlabel='Synthetic Turf')\n\nplt.xticks(index + bar_width/2, weather_injury_synth.Weather.values)\nplt.xlabel('Weather')\nplt.ylabel('Number of Injuries')\nplt.legend()\nadd_value_labels(ax)\nplt.title('Number of Injuries for Different Weather Conditions', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_injury_nat = game_injury_df[(game_injury_df['StadiumType'] == 'Outdoor') & (game_injury_df['FieldType'] == 'Natural')]\\\n.groupby(by=['Weather']).mean().reset_index().sort_values(by=['Weather'])\n\nweather_injury_synth = game_injury_df[(game_injury_df['StadiumType'] == 'Outdoor') & (game_injury_df['FieldType'] == 'Synthetic')]\\\n.groupby(by=['Weather']).mean().reset_index().sort_values(by=['Weather'])","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# visualize the staduin type for natural and synthetic turf\nn_groups = len(weather_injury_synth)\nfig, ax = plt.subplots(figsize=(7,7))\nindex = np.arange(n_groups)\nbar_width = 0.35\nopacity = 0.8\n\nrects1 = plt.bar(index, weather_injury_nat.Injury.values, bar_width,\nalpha=opacity,\ncolor='b',\nlabel='Natural Turf')\n\nrects2 = plt.bar(index + bar_width, weather_injury_synth.Injury.values, bar_width,\nalpha=opacity,\ncolor='g',\nlabel='Synthetic Turf')\n\nplt.xticks(index + bar_width/2, weather_injury_synth.Weather.values)\nplt.xlabel('Weather')\nplt.ylabel('Probability of Injury')\nplt.legend()\nplt.title('Probability of Injury for Different Weather Conditions', fontsize=16)\nadd_value_labels(ax, decimals = 3)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The highest probability of injury is:\n* for rainy weather on natural turf;\n* sunny weather on synthetic turf.\n\nSo there are two questions here:\n1. Could the natural turf be more slippery when it is raining? Do we have a higher chance of injury on rainy days for natural turf?\n2. Is synthetic turf more dangerous on sunny days? Does the high temperature also raise chances of injury on synthetic turf?"},{"metadata":{},"cell_type":"markdown","source":"Do some hypothesis testing to prove that the risk of injury on natural turf on rainy days is higher. The null hypothesis states that the probability of injury on natural turf on rainy days is lower or equal than the injury probability when those conditions are not met."},{"metadata":{"trusted":true},"cell_type":"code","source":"# get number of trials and overall injury rate under null\nhyp_df = game_injury_df[['StadiumType', 'FieldType', 'Weather', 'Injury']][game_injury_df['StadiumType'] == 'Outdoor']\n\nn_other = len(hyp_df[(hyp_df.FieldType != 'Natural') & (hyp_df.Weather != 'Rain')])\n\nn_natural_rain = len(hyp_df[(hyp_df.FieldType == 'Natural') & (hyp_df.Weather == 'Rain')])\n\np_null = len(hyp_df[(hyp_df.FieldType != 'Natural') & (hyp_df.Weather != 'Rain')]) / len(hyp_df)\n\np_injury_1 = hyp_df[(hyp_df.FieldType == 'Natural') & (hyp_df.Weather == 'Rain')].Injury.mean()\np_injury_0 = hyp_df[(hyp_df.FieldType != 'Natural') & (hyp_df.Weather != 'Rain')].Injury.mean()\n\n# compute standard error\nse_p = np.sqrt(p_null * (1-p_null) * (1/n_other + 1/n_natural_rain))\n\n#  compute z-score and p-value\nz = (p_injury_1 - p_injury_0) / se_p\n\nprint('The z-score is: {}'.format(z))\nprint('The p-value is: {}'.format(1-ss.norm.cdf(z)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"P-value is too high, and we can't reject the null hypothesis. There is no statistical significance that natural turf is riskier on rainy weather.\n\nNow let's try to look at the synthetic turf at high temperatures and clear weather."},{"metadata":{"trusted":true},"cell_type":"code","source":"# get the distribution for temperature for plays associated with the injury on synthetic turf\ntemp = game_injury_df[(game_injury_df.StadiumType == 'Outdoor') & (game_injury_df.FieldType == 'Synthetic')& (game_injury_df.Injury == 1)]\\\n.Temperature.values","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# plot the temperature distribution for injury on synthetic turf\nplt.figure(figsize=(10,7))\nplt.title('Injury vs Temperature for Synthetic Turf', fontsize=16)\nplt.xlabel('temperature')\nplt.ylabel('number of injuries')\nplt.hist(temp)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I can't see any significant spike for the number of injuries for high temperatures. I don't think that extremely high temperatures can cause a higher risk of injury. On the contrary, we see that there are more injuries at lower temperatures.\n\nLet's try hypothesis testing to prove that there are higher chances of injury on sunny days and synthetic turf:"},{"metadata":{"trusted":true},"cell_type":"code","source":"n_other = len(hyp_df[(hyp_df.FieldType != 'Synthetic') & (hyp_df.Weather != 'Clear')])\nn_synth_clear = len(hyp_df[(hyp_df.FieldType == 'Synthetic') & (hyp_df.Weather == 'Clear')])\n\np_null = len(hyp_df[(hyp_df.FieldType != 'Synthetic') & (hyp_df.Weather != 'Clear')]) / len(hyp_df)\n\np_injury_1 = hyp_df[(hyp_df.FieldType == 'Synthetic') & (hyp_df.Weather == 'Clear')].Injury.mean()\np_injury_0 = hyp_df[(hyp_df.FieldType != 'Synthetic') & (hyp_df.Weather != 'Clear')].Injury.mean()\n\n# compute standard error\nse_p = np.sqrt(p_null * (1-p_null) * (1/n_other + 1/n_synth_clear))\n\n#  compute z-score and p-value\nz = (p_injury_1 - p_injury_0) / se_p\n\nprint('The z-score is: {}'.format(z))\nprint('The p-value is: {}'.format(1-ss.norm.cdf(z)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The p-value is too high. We can't reject the null hypothesis.\n\nMy analysis shows that weather conditions interacting with field types don't lead to a higher risk of injury. "},{"metadata":{},"cell_type":"markdown","source":"`2` Synthetic and natural turf and player positions:"},{"metadata":{},"cell_type":"markdown","source":"Let's plot the probability of injury on synthetic turf for different player positions:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# get the probabilities of injury for different player positions\n\n# natural turf\npos_injury_nat = play_injury_df[(play_injury_df['FieldType'] == 'Natural')]\\\n.groupby(by=['Position']).mean().reset_index().sort_values(by=['Position'])\n\n# synthetic turf\npos_injury_synth = play_injury_df[(play_injury_df['FieldType'] == 'Synthetic')]\\\n.groupby(by=['Position']).mean().reset_index().sort_values(by=['Position'])","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"n_groups = len(pos_injury_nat)\nfig, ax = plt.subplots(figsize=(15,7))\nindex = np.arange(n_groups)\nbar_width = 0.35\nopacity = 0.8\n\nrects1 = plt.bar(index, pos_injury_nat.Injury.values*100, bar_width,\nalpha=opacity,\ncolor='b',\nlabel='Natural Turf')\n\nrects2 = plt.bar(index + bar_width, pos_injury_synth.Injury.values*100, bar_width,\nalpha=opacity,\ncolor='g',\nlabel='Synthetic Turf')\n\nplt.xticks(index + bar_width/2, pos_injury_synth.Position.values, rotation=90)\nplt.xlabel('Position')\nplt.ylabel('Probability of Injury, %')\nplt.legend()\nplt.title('Probability of Injury for Different Player Positions, %', fontsize=16)\nadd_value_labels(ax, decimals = 1)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It looks like players in DB (defensive back) position have higher chances of injury on synthetic turf!\n\nLet's do some hypothesis testing to prove that there are higher chances of injury for DB players on synthetic turf:"},{"metadata":{"trusted":true},"cell_type":"code","source":"hyp_df = play_injury_df[['FieldType', 'Position', 'Injury']]\n\nn_other = len(hyp_df[(hyp_df.FieldType != 'Synthetic') & (hyp_df.Position != 'DB')])\nn_synth_db = len(hyp_df[(hyp_df.FieldType == 'Synthetic') & (hyp_df.Position == 'DB')])\n\np_null = len(hyp_df[(hyp_df.FieldType != 'Synthetic') & (hyp_df.Position != 'DB')]) / len(hyp_df)\n\np_injury_1 = hyp_df[(hyp_df.FieldType == 'Synthetic') & (hyp_df.Position == 'DB')].Injury.mean()\np_injury_0 = hyp_df[(hyp_df.FieldType != 'Synthetic') & (hyp_df.Position != 'DB')].Injury.mean()\n\n# compute standard error\nse_p = np.sqrt(p_null * (1-p_null) * (1/n_other + 1/n_synth_db))\n\n#  compute z-score and p-value\nz = (p_injury_1 - p_injury_0) / se_p\n\nprint('The z-score is: {}'.format(z))\nprint('The p-value is: {}'.format(1-ss.norm.cdf(z)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Unfortunately, we can't reject the null hypothesis.\n\nBut let's look at the motion features for the DB position and compare them for different field types."},{"metadata":{"trusted":true},"cell_type":"code","source":"db_motion = motion_df_cleaned[motion_df_cleaned.Position == 'DB']","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def compare_distributions(db_motion, field1, field2, title=''):\n    fig, axs = plt.subplots(3,2,figsize=(15,10), sharex=True)\n\n    axs[0,0].hist(db_motion[db_motion.FieldType == 'Natural'][field1], bins = 20)\n    axs[0,0].set_title('Natural: {}'.format(field1))\n\n    axs[1,0].hist(db_motion[db_motion.FieldType == 'Synthetic'][field1], bins = 20)\n    axs[1,0].set_title('Synthetic: {}'.format(field1))\n\n    axs[0,1].hist(db_motion[db_motion.FieldType == 'Natural'][field2], bins = 20)\n    axs[0,1].set_title('Natural: {}'.format(field2))\n\n    axs[1,1].hist(db_motion[db_motion.FieldType == 'Synthetic'][field2], bins = 20)\n    axs[1,1].set_title('Synthetic: {}'.format(field2))\n\n    sns.kdeplot(db_motion[db_motion.FieldType == 'Natural'][field1], ax = axs[2,0], label='natural')\n    sns.kdeplot(db_motion[db_motion.FieldType == 'Synthetic'][field1],ax = axs[2,0], label='synthetic')\n    axs[2,0].set_title('{} Distribution'.format(field1))\n\n    sns.kdeplot(db_motion[db_motion.FieldType == 'Natural'][field2], ax = axs[2,1], label='natural')\n    sns.kdeplot(db_motion[db_motion.FieldType == 'Synthetic'][field2],ax = axs[2,1], label='synthetic')\n    axs[2,1].set_title('{} Distribution'.format(field2))\n\n    plt.suptitle(title, fontsize=16)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_distributions(db_motion, 'dis_max', 'dis_avg', title='Distance Distribution for DB')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_distributions(db_motion, 's_max', 's_avg', title='Distance Distribution for DB')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_distributions(db_motion, 'angle_max', 'angle_avg', title='Distance Distribution for DB')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"When we look at player movement in the defensive back position we see:\n* The distributions for distances traveled are almost the same for natural and synthetic turf;\n* The angle and speed distributions for the synthetic and natural turf fields are slightly different. A higher speed spike for synthetic turf is especially interesting. "},{"metadata":{},"cell_type":"markdown","source":"Let's look at the injury on synthetic turf and DB player position:"},{"metadata":{"trusted":true},"cell_type":"code","source":"db_motion[db_motion.Injury == 1][['RosterPosition', 'StadiumType', 'FieldType','PlayType', 'Position',\n       'PositionGroup', 'time_max', 'dir_max', 'dis_max', 'o_max', 's_max',\n       'angle_max', 'time_avg', 'dir_avg', 'dis_avg', 'o_avg', 's_avg',\n       'angle_avg', 'Injury']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have only one injury record for DB position on synthetic turf, and we see that the maximum speed is `10.11`, which is quite high! And as we noticed, players in the DB position tend to have higher speed on synthetic turf. __Possibly, higher speeds on synthetic turf could lead to an elevated risk of injury.__\n\nLet's look at players in the ILB (Inside linebackers) position:"},{"metadata":{"trusted":true},"cell_type":"code","source":"db_motion = motion_df_cleaned[motion_df_cleaned.Position == 'ILB']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_distributions(db_motion, 's_max', 's_avg', title='Distance Distribution for ILB')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The distributions for speed for ILB position are almost identical. But we still see a little spike for high speeds for synthetic turf.\n\nLet's look closely at the injury records:"},{"metadata":{"trusted":true},"cell_type":"code","source":"db_motion[db_motion.Injury == 1][['RosterPosition', 'StadiumType', 'FieldType','PlayType', 'Position',\n       'PositionGroup', 'time_max', 'dir_max', 'dis_max', 'o_max', 's_max',\n       'angle_max', 'time_avg', 'dir_avg', 'dis_avg', 'o_avg', 's_avg',\n       'angle_avg', 'Injury']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we can see that only one of three injury records corresponds to the high speed of `9.20`."},{"metadata":{},"cell_type":"markdown","source":"## Conclusions:\n* There is a higher risk of injury when playing on the synthetic turf\n* Under severe weather conditions (rain or snow) players on the synthetic turf tend to have a lower speed than the players on natural turf\n* Probably, the synthetic turf can't sustain the high speed. It makes the players slow down when the weather is bad. Probably, moving too fast on synthetic turf could even raise the risk of injury when playing in a certain position"},{"metadata":{},"cell_type":"markdown","source":"## References and Credits\n1. [Great Medium article](https://towardsdatascience.com/the-search-for-categorical-correlation-a1cf7f1888c9) on correlation for categorical variables. I highly encourage to read it!\n2. [dython library](https://github.com/shakedzy/dython/blob/master/dython/nominal.py): The implementation of Theil's U"},{"metadata":{},"cell_type":"markdown","source":"__I will further update this notebook, so, please, feel free to ask questions and add comments!__"}],"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":4,"nbformat_minor":1}