{"cells":[{"metadata":{},"cell_type":"markdown","source":"## 1. Efficient Net download\nInstall efficientnet downloaded from github. If using Kaggle, please switch on internet."},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"# output is hidden in view version\n!pip3 install git+https://github.com/qubvel/efficientnet","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2. Import packages\n<span style = \"color:red\"> Could we remove the imports that aren't being used? </span>"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# Regular Imports\nimport numpy as np\nimport pandas as pd\nimport os\nimport random\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm \nfrom PIL import Image\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import GroupKFold, KFold, StratifiedKFold\nimport seaborn as sns\nimport math\nimport cv2\nimport pydicom\nimport os\nimport glob\nimport pickle as pkl\nimport matplotlib.image as mpimg\nfrom tabulate import tabulate\nimport missingno as msno \nfrom IPython.display import display_html\nfrom PIL import Image\nimport gc\nfrom skimage.transform import resize\nimport copy\nimport re\nfrom scipy.stats import pearsonr\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n%matplotlib inline\n\n# Segmentation\nimport glob\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nimport scipy.ndimage\nfrom skimage import morphology\nfrom skimage import measure\nfrom skimage.transform import resize\nfrom sklearn.cluster import KMeans\nfrom plotly import __version__\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nfrom plotly.tools import FigureFactory as FF\nfrom plotly.graph_objs import *\ninit_notebook_mode(connected=True) \n\n# Model imports\nimport tensorflow as tf \nfrom tensorflow.keras.layers import (\n                                    Dense, Dropout, Activation, Flatten, Input, BatchNormalization, GlobalAveragePooling2D,\n                                    Add, Conv2D, AveragePooling2D, LeakyReLU, Concatenate , Lambda\n                                    )\nfrom tensorflow.keras import Model\nimport tensorflow.keras.models as M\nfrom tensorflow.keras.utils import Sequence\nimport tensorflow.keras.backend as K\n# import tensorflow.keras.applications as tfa\nimport tensorflow_addons as tfa\nimport efficientnet.tfkeras as efn\nfrom sklearn.model_selection import train_test_split, KFold\nimport seaborn as sns\nimport plotly.express as px\n\n\n\npd.set_option(\"display.max_columns\", 100)\ncustom_colors = ['#74a09e','#86c1b2','#98e2c6','#f3c969','#f2a553', '#d96548', '#c14953']\nsns.palplot(sns.color_palette(custom_colors))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Tuned Hyperparameters\nWe define the number of -\n* **epochs** : the number of times we loop through the dataset. \n* **batch_size**: how many training examples to feed into network before updating the weights and internal nodes\n* **LR**: learning rate\n* **MODEL_CLASS**:  allow us to define which efficientnet model we want to use\n* **SAVE_BEST**: defined as true, to only save the model due to 'early stopping' implemented"},{"metadata":{"trusted":true},"cell_type":"code","source":"EPOCHS = 40\nBATCH_SIZE = 8\nNFOLD = 5\nLR = 0.003\nSAVE_BEST = True\nMODEL_CLASS = 'b1'\npath = '../input/osic-pulmonary-fibrosis-progression'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4. Train-Test split\nWe create a test dataset with 20% of unique patients and their associated data; 80% for training. \nDuplicates are dropped as some patients have more than one recorded FVC value for each week. "},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data = pd.read_csv(f'{path}/train.csv') \nall_data.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\nunique_patients = all_data.Patient.unique()\ntrain_ids, test_ids = train_test_split(unique_patients, test_size=0.2, random_state=42)\ntrain = all_data[all_data['Patient'].isin(train_ids)]\ntest = all_data[all_data['Patient'].isin(test_ids)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.groupby('Sex').agg('count')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.scatter(train, x=\"Sex\", y =\"FVC\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = px.scatter(test, x=\"Sex\", y =\"FVC\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# BAD_ID = ['ID00011637202177653955184', 'ID00052637202186188008618']\n# a = test[\"Patient\"].values\n# for b in BAD_ID:\n#     if b in a:\n#         print(True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 5. Data Preprocessing\nThe original features of the data included Age, Smoking_Status, Sex, Weeks, Percent. However, we would like to transform the data, and include altered features. \n\nThe features that will be fed into the model eventually include:\n* baseline Age (numerical)\n* baseline Percent (numerical)\n* Gender (encoded)\n* Smoking status (encoded)\n\nAlong with that, the numerical data is normalised and categorical variables are one-hot encoded."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_baseline_week(df):    \n    # make a copy to not change original df    \n    _df = df.copy()\n    # ensure all Weeks values are INT and not accidentaly saved as string\n    _df['Weeks'] = _df['Weeks'].astype(int)\n    _df['min_week'] = _df['Weeks']\n    _df[\"min_week\"] = _df.groupby('Patient')['Weeks'].transform('min')\n    _df['baselined_week'] = _df['Weeks'] - _df['min_week']\n    \n    return pd.DataFrame(_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = get_baseline_week(train)\n#train = get_baseline_FVC_new(train)\ntest = get_baseline_week(test)\n#test = get_baseline_FVC_new(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# define which attributes shall not be transformed, are numeric or categorical\nno_transform_attribs = ['Patient', 'Weeks', 'min_week', 'FVC']\nnum_attribs = ['Percent', 'Age']\ncat_attribs = ['Sex', 'SmokingStatus']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import preprocessing\n\nnum_attribs_scld = [s + '_scld' for s in num_attribs]\n\nmin_max_scaler = preprocessing.MinMaxScaler()\ntrain[num_attribs_scld] = min_max_scaler.fit_transform(train[num_attribs])\ntest[num_attribs_scld] = min_max_scaler.transform(test[num_attribs])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Encoding our categorical variables. \n* Gender: Male: 0, Female: 1\n* Smoking_Status: Never Smoked: [0,0], Ex-Smoker: [1,1], Currently Smokes: [0,1]\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_tab_scaled(df): #getting scaled variables\n    vector = [df.Percent_scld.values[0]] # only the first percent value is retained, so this forms the \"base_percent\"\n    vector.extend([df.Age_scld.values[0]]) # only the first age value is retained, so this forms the \"base age\"\n    if df.Sex.values[0].lower() == 'male':\n        vector.append(0)\n    else:\n        vector.append(1)\n    if df.SmokingStatus.values[0] == 'Never smoked':\n        vector.extend([0,0])\n    elif df.SmokingStatus.values[0] == 'Ex-smoker':\n        vector.extend([1,1])\n    elif df.SmokingStatus.values[0] == 'Currently smokes':\n        vector.extend([0,1])\n    else:\n        vector.extend([1,0])\n    \n    return np.array(vector)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We need to create the labels i.e y values for our dataset. We aim to predict the rate of FVC values deterioration, which would be the gradient of a linear regression model fit to every available patient. "},{"metadata":{"trusted":true},"cell_type":"code","source":"A = {} # Gradient of the linear regression of FVC against weeks \nTAB = {} # Initialize tabular data for each patient\nP = [] # Patient IDs\n\n# for all 140 train patients we compute the gradient \nfor i, p in tqdm(enumerate(train.Patient.unique())):\n    sub = train.loc[train.Patient == p, :] \n    fvc = sub.FVC.values    \n    weeks = sub.Weeks.values \n    c = np.vstack([weeks, np.ones(len(weeks))]).T\n    a, b = np.linalg.lstsq(c, fvc)[0]\n    \n    A[p] = a\n    TAB[p] = get_tab_scaled(sub)\n    # TAB[p] = get_tab_unscaled(sub)\n    P.append(p)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_grad = train.copy()\ntrain_grad['gradient'] = train_grad['Patient'].map(A)\ntrain_grad.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_grad.loc[train_grad.gradient == train_grad.gradient.max()]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.lineplot(data = train_grad.loc[train_grad.Patient == 'ID00197637202246865691526'].FVC)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set_style(\"whitegrid\") \nsns.boxplot(x = 'SmokingStatus', y = 'gradient', data = train_grad.drop_duplicates(subset=['Patient'])) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.boxplot(x = 'Sex', y = 'gradient', data = train_grad.drop_duplicates(subset=['Patient'])) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(train_grad['Age'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_grad['Age bin'] =  ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For each sliced ct scan, we need to resize into 512x512 pixel dimensions. "},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_img(path):\n    d = pydicom.dcmread(path)\n    return cv2.resize((d.pixel_array - d.RescaleIntercept) / (d.RescaleSlope * 1000), (512, 512))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are patients (labelled with BAD_IDs ( <span style = \"color:red\"> why are they bad_IDs maybe we need to explain that ? </span> )) that need to be removed from the training data. This class also transforms the dataframes into the required numpy arrays so that our data can be input into the neural network. "},{"metadata":{"trusted":true},"cell_type":"code","source":"class IGenerator(Sequence):\n    BAD_ID = ['ID00011637202177653955184', 'ID00052637202186188008618']\n    def __init__(self, keys, a, tab, batch_size=BATCH_SIZE):\n        self.keys = [k for k in keys if k not in self.BAD_ID] # all the patients in the keys  \n        self.a = a # gradients found by fitting a linear regression model\n        self.tab = tab # tabular data for patient specified in keys\n        self.batch_size = batch_size\n        \n        self.train_data = {} # initialise a dictionary to contain all the images pertaining to one key: patient id \n        for p in train.Patient.values:\n            self.train_data[p] = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/')\n    \n    def __len__(self):\n        return 1000 # 1000 batches per epoch\n    \n    def __getitem__(self, idx):\n        x = []\n        a, tab = [], [] \n        keys = np.random.choice(self.keys, size = self.batch_size) # randomly chosen n patients for one batch  \n        for k in keys:\n            try:\n                i = np.random.choice(self.train_data[k], size=1)[0]\n                img = get_img(f'../input/osic-pulmonary-fibrosis-progression/train/{k}/{i}')\n                x.append(img)\n                a.append(self.a[k])\n                tab.append(self.tab[k])\n            except:\n                print(k, i)\n       \n        x,a,tab = np.array(x), np.array(a), np.array(tab)\n        x = np.expand_dims(x, axis=-1)\n        return [x] , a","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 7. Model Definition"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_efficientnet(model, shape):\n    '''\n    From https://github.com/qubvel/efficientnet\n    EfficientNet is a CNN architecture achieving state of the art accuracy.\n    b0 is the simplest model, b7 is the most complex.\n    '''\n    models_dict = {\n        'b0': efn.EfficientNetB0(input_shape=shape,weights=None,include_top=False),\n        'b1': efn.EfficientNetB1(input_shape=shape,weights=None,include_top=False), # We use a b1 efficientnet\n        'b2': efn.EfficientNetB2(input_shape=shape,weights=None,include_top=False),\n        'b3': efn.EfficientNetB3(input_shape=shape,weights=None,include_top=False),\n        'b4': efn.EfficientNetB4(input_shape=shape,weights=None,include_top=False),\n        'b5': efn.EfficientNetB5(input_shape=shape,weights=None,include_top=False),\n        'b6': efn.EfficientNetB6(input_shape=shape,weights=None,include_top=False),\n        'b7': efn.EfficientNetB7(input_shape=shape,weights=None,include_top=False)\n    }\n    return models_dict[model]\n\ndef build_model(shape=(512, 512, 1), model_class=None):\n    inp = Input(shape=shape) # 512 x 512 input shape\n    base = get_efficientnet(model_class, shape) # A b1 pre-trained efficientnet is used\n    x = base(inp)\n    x = GlobalAveragePooling2D()(x)\n    \n    inp2 = Input(shape=(4,))\n    out_tab = tf.keras.layers.GaussianNoise(0.2)(inp2) # add some noise to our data\n    \n    x2 = Concatenate()([x, out_tab]) \n    \n    #### added ###\n    den_1 = Dense(20)(x) # linear layer with 50 output nodes\n    den_1 = LeakyReLU(alpha=0.3)(den_1) # default alpha is 0.3\n    den_2 = Dense(500)(den_1)\n    den_2 = LeakyReLU(alpha=0.3)(den_2)\n    bn_1 = BatchNormalization()(den_2)\n    den_3 = Dense(100)(bn_1)\n    \n    x3 = Dropout(0.5)(den_3) # move the dropout layer to before the linear layers \n    \n    # the linear layers here are akin to the final dense layer usually used in a \"before concat\" model\n    out_1 = Dense(1)(x3)\n    out_2 = Dense(1, activation='relu')(x3)\n    y = out_1 + tf.keras.backend.cumsum(out_2, axis=1)\n    \n    model = Model([inp], y)\n    return model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 8. Model Training\n\nWe split the dataset into 5 folds for cross validation purposes. There is early stopping employed in the model and the loss function that we decided upon is RMSE = Root Mean Squared Error. It is a commonly used loss function in regression problems.  <span style = \"color:red\"> do we need to change mse in the code to rmse ? </span>"},{"metadata":{"trusted":true},"cell_type":"code","source":"kf = KFold(n_splits=NFOLD, random_state=42,shuffle=False)\nP = np.array(P)\nsubs = []\nfolds_history = []\nfor fold, (tr_idx, val_idx) in enumerate(kf.split(P)):\n    print('#####################')\n    print('####### Fold %i ######'%fold)\n    print('#####################')\n    print('Training...')\n    \n    er = tf.keras.callbacks.EarlyStopping(\n        monitor=\"val_loss\",\n        min_delta=1e-3,\n        patience=8,\n        verbose=1,\n        mode=\"auto\",\n        baseline=None,\n        restore_best_weights=True,\n    )\n\n    cpt = tf.keras.callbacks.ModelCheckpoint(\n        filepath='fold-%i.h5'%fold,\n        monitor='val_loss', \n        verbose=1, \n        save_best_only=SAVE_BEST,\n        mode='auto'\n    )\n\n    rlp = tf.keras.callbacks.ReduceLROnPlateau(\n        monitor='val_loss', \n        factor=0.5,\n        patience=5, \n        verbose=1, \n        min_lr=1e-8\n    )\n    model = build_model(model_class=MODEL_CLASS)\n    model.compile(optimizer=tf.keras.optimizers.Adam(learning_rate=LR), loss=\"mse\") \n    history = model.fit_generator(IGenerator(keys=P[tr_idx], \n                                   a = A, \n                                   tab = TAB), \n                        steps_per_epoch = 32,\n                        validation_data=IGenerator(keys=P[val_idx], \n                                   a = A, \n                                   tab = TAB),\n                        validation_steps = 16, \n                        callbacks = [cpt, rlp], \n                        epochs=EPOCHS)\n    folds_history.append(history.history)\n    print('Training done!')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 9. Validation Loss examination"},{"metadata":{"trusted":true},"cell_type":"code","source":"min_array = []\nfor i in range(5):\n    min_array.append(min(folds_history[i]['val_loss']))\n    print(i, min(folds_history[i]['val_loss']))\nprint(min(min_array))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We pick the best model (weights) based on cross validation score.\nif SAVE_BEST:\n    mean_val_loss = np.mean([np.min(h['val_loss']) for h in folds_history])\nelse:\n    mean_val_loss = np.mean([h['val_loss'][-1] for h in folds_history])\nprint('Our mean CV MAE is: ' + str(mean_val_loss))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"min_fold finds out which fold gives the least validation loss"},{"metadata":{"trusted":true},"cell_type":"code","source":"min_fold = np.argmin([np.min(h['val_loss']) for h in folds_history])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"min_fold = 1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"^print what is the fold that gives min loss, keep a record below "},{"metadata":{"trusted":true},"cell_type":"code","source":"################################\n# min_fold = 1 #change accordingly\n################################\n#uncomment for using notebook","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"1. download the weights from the output and save to zip file and upload them \n2. Restart kernel at this point for fitting weights from best fold to model in model building"},{"metadata":{},"cell_type":"markdown","source":"building model for prediction"},{"metadata":{},"cell_type":"markdown","source":"## 10. Model for Test Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_pred_model(shape=(512, 512, 1), model_class=None, fold=None):\n    inp = Input(shape=shape) # 512 x 512 input shape\n    base = get_efficientnet(model_class, shape) # A b1 pre-trained efficientnet is used\n    x = base(inp)\n    x = GlobalAveragePooling2D()(x)\n    \n    inp2 = Input(shape=(5,))\n    out_tab = tf.keras.layers.GaussianNoise(0.2)(inp2) # add some noise to our data\n    \n    x2 = Concatenate()([x, out_tab]) \n    \n    #### added ###\n    den_1 = Dense(20)(x2) # linear layer with 50 output nodes\n    den_1 = LeakyReLU(alpha=0.3)(den_1) # default alpha is 0.3\n    den_2 = Dense(500)(den_1)\n    den_2 = LeakyReLU(alpha=0.3)(den_2)\n    bn_1 = BatchNormalization()(den_2)\n    den_3 = Dense(100)(bn_1)\n    \n    x3 = Dropout(0.5)(den_3) # move the dropout layer to before the linear layers \n    \n    # the linear layers here are akin to the final dense layer usually used in a \"before concat\" model\n    out_1 = Dense(1)(x3)\n    out_2 = Dense(1, activation='relu')(x3)\n    y = out_1 + tf.keras.backend.cumsum(out_2, axis=1)\n    \n    model = Model([inp, inp2], y)\n    \n    # Take from kaggle  working output\n    weights = [w for w in os.listdir('../input/best-model-fold/') if str(fold) in w][0] #use dir of './' if using training above, else can just use my trained weights in input\n    model.load_weights('../input/best-model-fold/' + weights) #use dir of './' if using training above, else can just use my trained weights in input\n    \n    #take weights from uploaded weights\n    #uncomment when using notebook, instead of when committing notebook\n    #weights = [w for w in os.listdir('../input/scaled-weights-genderagepercentsmoker/') if str(fold) in w][0] #use dir of './' if using training above, else can just use my trained weights in input\n    #model.load_weights('../input/scaled-weights-genderagepercentsmoker/' + weights) #use dir of './' if using training above, else can just use my trained weights in input\n    return model\n#models = [build_pred_model(shape=(512, 512, 1), model_class='b1', fold=min_fold)]\n#uncomment the above to build model from the weights trained above, else can use the below code for building model\nmodels = [build_pred_model(shape=(512, 512, 1), model_class='b1', fold=min_fold)]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":false,"trusted":true},"cell_type":"code","source":"subs = []\nq = 0.5\nweeks = list(range(-12, 134))\npredictions = []\nfor model in models:\n    metric = []\n    \n    A_test, B_test, P_test,W, FVC= {}, {}, {},{},{} \n    STD, WEEK = {}, {} \n    for p in test.Patient.unique():\n        x = [] \n        tab = [] \n        ldir = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/')\n        for i in ldir:\n            if int(i[:-4]) / len(ldir) < 0.8 and int(i[:-4]) / len(ldir) > 0.15: # what is this? only certain slices are being extracted for patients?\n                x.append(get_img(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/{i}')) \n                tab.append(get_tab_scaled(test.loc[test.Patient == p, :])) \n        if len(x) <= 1:\n            continue\n        print(len(test[test.Patient == p])) # number of fvc values to predict for every week\n        tab = np.array(tab) \n        print(\"number of patient data is: \", len(tab))\n        x = np.expand_dims(x, axis=-1) \n        _a = model.predict([x, tab]) # Predict from all image data and tabular data.\n        a = np.quantile(_a, q) # Get the value at the 50th percentile\n        \n        A_test[p] = a\n        print(\"number of patient FVC data is: \", len(test.FVC.values[test.Patient == p])) # number of times patient takes\n        #######################################################################################################################################\n        B_test[p] = test.FVC.values[test.Patient == p] - a*test.Weeks.values[test.Patient == p] #to find the y intercept   #\n        #######################################################################################################################################\n        print(\"B_test[p] is: \", B_test[p])\n        #B_test is actually to find the intercept\n        \n        #P_test[p] = test.Percent_scld.values[test.Patient == p] \n        #print(\"P_test[p] is: \", P_test[p])\n        #WEEK[p] = test.baselined_week_scld.values[test.Patient == p]\n        #print(\"WEEK[p] is: \", WEEK[p])\n    \n    for p in test.Patient.unique():\n        for w in weeks:\n            patient_prediction = {}\n            fvc = A_test[p] * w + B_test[p]  #y = mx + c, A_test[p] = gradient of progression for patient p, w = week_num, B_test[p] is the calculated y_intercept\n            print(\"fvc is: \", fvc)\n            patient_prediction = {\n                'Week': w,\n                'Patient': p,\n                'FVC': np.sum(fvc)/len(fvc)\n            }\n            predictions.append(patient_prediction)\nprint('done!')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# GET RESULTS INTO A DATAFRAME \npredictions_df = pd.DataFrame(predictions)\n# predictions_df.to_csv(\"test_predictions_GAPS_scaled.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"pred_patients are unique patientID of patients in the predictions "},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_patients  = predictions_df.Patient.unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions_df.rename(columns = {'Week': 'Weeks'}, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Meso Analysis to check for model biases"},{"metadata":{"trusted":true},"cell_type":"code","source":"test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"merged_data = pd.merge(predictions_df, test, on = ['Patient', 'Weeks'], how = 'inner')\nmerged_data = merged_data.rename(columns={\"FVC_x\": \"Predicted FVC\", \"FVC_y\": \"True FVC\"})\nmerged_data[\"Absolute Error\"] = abs(merged_data[\"Predicted FVC\"] - merged_data[\"True FVC\"])\nmerged_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def calculate_rmse(df):\n    return np.sqrt((1/len(df)) * sum(np.square(df[\"True FVC\"] - df[\"Predicted FVC\"])))  #change this part for non_scaled","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# RMSE by sex\nmales = merged_data[merged_data['Sex'] == \"Male\"]\nfemales = merged_data[merged_data['Sex'] == \"Female\"]\n\nmales_rmse = calculate_rmse(males)\nfemales_rmse = calculate_rmse(females)\nprint(f\"There are {males.shape[0]} males of avg rmse {males_rmse} and {females.shape[0]} of avg rmse {females_rmse}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# PLot of absolute errors\nfig = px.scatter(merged_data, x=\"Sex\", y =\"Absolute Error\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# PLot of absolute errors\nfig = px.scatter(merged_data, x=\"SmokingStatus\", y =\"Absolute Error\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.groupby(\"SmokingStatus\").agg(\"count\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Rmse by smoking status\nmerged_data.groupby(\"SmokingStatus\").agg(\"count\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"smoke_current = merged_data[merged_data[\"SmokingStatus\"] == \"Currently smokes\"]\nsmoke_ex = merged_data[merged_data[\"SmokingStatus\"] == \"Ex-smoker\"]\nsmoke_never = merged_data[merged_data[\"SmokingStatus\"] == \"Never smoked\"]\n\nsmoke_current_rmse = calculate_rmse(smoke_current)\nsmoke_ex_rmse = calculate_rmse(smoke_ex)\nsmoke_never_rmse = calculate_rmse(smoke_never)\nprint(f\"Current {smoke_current.shape[0]}: {smoke_current_rmse}, ex {smoke_ex.shape[0]}: {smoke_ex_rmse}, never {smoke_never.shape[0]}: {smoke_never_rmse}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# code from https://www.kaggle.com/piantic/osic-pulmonary-fibrosis-progression-basic-eda\nplt.figure(figsize=(16, 6))\nsns.kdeplot(merged_data.loc[merged_data['SmokingStatus'] == 'Ex-smoker', 'Age'], label = 'Ex-smoker',shade=True)\nsns.kdeplot(merged_data.loc[merged_data['SmokingStatus'] == 'Never smoked', 'Age'], label = 'Never smoked',shade=True)\nsns.kdeplot(merged_data.loc[merged_data['SmokingStatus'] == 'Currently smokes', 'Age'], label = 'Currently smokes', shade=True)\n\n# Labeling of plot\nplt.xlabel('Age (years)'); plt.ylabel('Density'); plt.title('Distribution of Ages');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"merged_data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"age_range1 = merged_data[merged_data[\"Age\"] <= 65]\nage_range2 = merged_data[(merged_data[\"Age\"] > 65) & (merged_data['Age'] <= 71)]\nage_range3 = merged_data[merged_data[\"Age\"] > 71]\n                          \na1_rmse = calculate_rmse(age_range1)\na2_rmse = calculate_rmse(age_range2)\na3_rmse = calculate_rmse(age_range3)\nprint(f\"{age_range1.shape[0]} {a1_rmse}, {age_range2.shape[0]} {a2_rmse}, {age_range3.shape[0]} {a3_rmse}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get train sizes\nprint(train[train['Age'] <= 65].shape[0])\nprint(train[(train[\"Age\"] > 65) & (train['Age'] <= 71)].shape[0])\nprint(train[train['Age'] > 71].shape[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# PLot of absolute errors\nfig = px.scatter(merged_data, x=\"Age\", y =\"Absolute Error\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rmse = calculate_rmse(merged_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# 144.413 (for pyotrch part 2 rmse)\n# around 142 for this version\nrmse ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#rmse for unscaled notebook with percent included = 143.2077201851251","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_test = merged_data[\"True FVC\"]\ny_pred = merged_data[\"Predicted FVC\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import mean_squared_error\nfrom math import sqrt\nrms = sqrt(mean_squared_error(y_test, y_pred))\nrms","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt # Impot the relevant module\n\n#fig, ax = plt.subplots() # Create the figure and axes object\n\nax = y_test.plot()\ny_pred.plot(ax=ax)\nplt.suptitle('CNN (efficientnet) + MLP best model predictions vs label')\nplt.legend(loc=\"upper left\")\nplt.xlabel('Patient Data Points')\nplt.ylabel('FVC score')\nax.set_xticklabels([])\nplt.savefig('foo.png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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":4}