<jupyter_start><jupyter_text>EDA + Preprocessing + Modeling Importing necessary packages<jupyter_code># Let's import some of the libraries for 
# data manipulation and visualization
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
# Loaing trai_data and demographic data
train = pd.read_csv('train.csv')
demo = pd.read_csv('demog.csv')
# Getting shape of data
print(train.shape,demo.shape)
train.head()
demo.head()<jupyter_output><empty_output><jupyter_text>Preparing combined Data<jupyter_code># First merge demo data with train data to have a complete data
all_data = train.merge(right=demo,on='HCP_ID',how='left')
all_data.shape<jupyter_output><empty_output><jupyter_text>EDA<jupyter_code># Let's undergo a univariate analysis of all of 
# our features to better understand the distribution 
# of our features
# ======================================================================================================
# ======================================================================================================
#**************************************************************************************************************************	
	
	
#  different types of utility functions used for various univariate analysis (monks advanced method for univariate analysis)
from matplotlib import pyplot as plt
import seaborn as sns


def plot(x=None, y=None, hue=None, data=None,name='',barplot=False,countplot=False,boxplot=True,factorplot=False):
    data=data.copy()
    plt.figure(figsize=(14,9))
    if type(data[x].values[1])==str:
        from sklearn.preprocessing import LabelEncoder
        le=LabelEncoder()
        data[x]=data[x].apply(lambda y:treat_cat(y))
        data[x]=le.fit_transform(data[x])
    else:
        me=np.mean(data[x])
        data[x]=data[x].apply(lambda y: treat_num(y,m=me))
    try:
        print(name,'\n')
        data[x].hist()
        plt.show()
    except:
        print('cannot plot histogram\n')
    
    try:
        print(name,'\n')
        sns.distplot(train[x])
        plt.show()
    except:
        print('cannot plot distplot\n')
        
    if barplot:
        try:
            print(name,'\n')
            sns.barplot(x=x,y=y,hue=hue,data=data)
            plt.xlabel(x)
            plt.ylabel(y)
            plt.show()
        except:
            print('cannot plot barplot\n')
    if countplot:        
        try:
            print(name,'\n')
            sns.countplot(data[x])
            plt.xlabel(x)
            plt.ylabel(y)
            plt.show()
        except:
            print('cannot plot countplot\n')
    if boxplot:        
        try:
            print(name,'\n')
            sns.boxplot(x=x,y=y,hue=hue,data=data)
            plt.xlabel(x)
            plt.ylabel(y)
            plt.show()
        except:
            print('cannot plot boxplot\n')  
    if factorplot:        
        try:
            print(name,'\n')
            data=pd.DataFrame(data)
            sns.factorplot(x=x,y=y,data=data,size=14)
            plt.xlabel(x)
            plt.ylabel(y)
            plt.show()
        except:
            print('cannot plot factorplot\n')  
            
def plotHeatMap(df, size=(16,9),an=True):
    plt.figure(figsize=size)
    cmap = plt.cm.viridis
    sns.heatmap(df.corr(),cmap=cmap,linewidths=0.2,vmax=1, linecolor="white", square=True,annot=an)
    plt.show()
    
def treat_cat(tmp):
    if tmp!=tmp:
        return 'unknown'
    else:
        return tmp
    
def max_min(data):
    data=pd.DataFrame(data)
    print('max = '+str(max(data.values.ravel()))+'  min = '+str(min(data.values.ravel())))

def data_type(data):
    data=pd.DataFrame(data)
    for i in data.values.ravel():
        if i==i:
            return type(i)

			
# ======================================================================================================
# ======================================================================================================
#**************************************************************************************************************************	

#	best method for univariate analysis			

def describe(data,x=None, y=None, with_plot=True):
    data = data.copy()
    null_values(data[x])
    #for uniqueness
    try:
        if len(set(data[x].values))==len(data[x].values):
            print('is_unique : Yes')
        else:
            print('is_unique : No')
    except:
        pass
    
    no_of_unique_values(data[x])
    if len(set(data[x].values.ravel()))<=10:
        print(set(data[x].values.ravel()))
        count(data=data[x])
        
    if data_type(data[x])==str:
        from sklearn.preprocessing import LabelEncoder
        le=LabelEncoder()
        data[x]=data[x].apply(lambda y:treat_cat(y))
        data[x]=le.fit_transform(data[x])
    else:
        me=np.mean(data[x])
        data[x]=data[x].apply(lambda y: treat_num(y,m=me))
    data[x]=pd.DataFrame(data[x])
    max_min(data[x])
    non_zero(data[x])
    print(data[x].describe(),'\n')
    if with_plot:
        plot(x=x,data=data,y=y)



# ======================================================================================================
# ======================================================================================================
#**************************************************************************************************************************	
		
def non_zero(data):
    data=pd.DataFrame(data)
    print('non zero values = '+str(len(data[data.values.ravel()!=0]))+
         '   zero values = '+str(len(data[data.values.ravel()==0])))

def drop(data=None,col=None):
    return data.drop([i for i in col],axis=1)
    
def null_values(data):
     print('null values = '+str(len(data[data!=data])))
        
def treat_num(x,m):
    if x!=x:
        return m
    else:
        return x

def no_of_unique_values(data):
    print('Number of unique values : ',len(set(data)))
    
def count(data):
    from collections import Counter
    print(Counter(data))
    
def consider_unknown(data,li):
    if data in li:
        return data
    else :
        return 'unknown'


for i in all_data.columns:
    print(i,' =================================================>')
    describe(data=all_data,x=i)<jupyter_output>HCP_ID  =================================================>
null values = 0
is_unique : Yes
Number of unique values :  37186
max = 46513  min = 9304
non zero values = 37186   zero values = 0
count    37186.000000
mean     27910.118297
std      10741.107350
min       9304.000000
25%      18612.250000
50%      27910.500000
75%      37210.750000
max      46513.000000
Name: HCP_ID, dtype: float64<jupyter_text>EDA Insights<jupyter_code># Some conclusions from the univariate analysis

# HCP_ID (Unique ID for each Doctor) : 
# It appears that HCP_ID is unique, there are no any errors as it should be unique as we have records for unique doctors

# RL(Rep_Live) : Its the contact with doctors via video conferencing
# There are 4849 missing values that we need to predict
# It appears that more than 75% of the doctors have affinity >0.51
# ie they are not having a good affinity via this channel
# Also from the histogram plot, we can easily observe that distribution is -ve skewed and we have a more
# population biased towards right, so it appears this channel is not so convincing to doctors

# P2P (Peer to Peer such as conferences etc)
# Wow, 29136 missing values, its huge, :O
# It looks interseting to further investiate it.
# btw, From the histogram plot, it appears that most of the population
# who used this channel are having good affinity as we have more than
# 75% of population belongs to affinty < ~ 0.19 

# OLV (Online video (such as online video, etc – Digital_Pull))
# This channel has also too many missing values : 17027
# From the histogram plot, it appears that most of the population
# who used this channel are having good affinity as we have more than
# 75% of population belongs to affinty < ~ 0.15 

# Remote representative (Medical rep via video conference – Rep_Remote)
# We have here 9648 missing values here
# distribution is +ve skewed ie we have more population which is
# biased towards left ie less affinity

# DRT Direct mail (Snail mail – Direct_Mail)
# Here we have 4492 missing values 
# and distribution seems highly +ve skewed
# we have 23789 doctors having very high affinty ie zero
# and more than 90% doctors have affinty <.2
# wow, a good channel

# DMS Direct/Digital message (such as mobile messages – Digital_Push)
# Here we have 4379 missing values
# and more than 90% of population are having affinity <0.2
# which makes it a good channel comparitvely


# OLA Online advertisement (Such as website banner – Digital_Pull)
# Here we are having 11170 missing values
# and its distribution seems +ve skewed,
# most of the populations (>80%) are having good affinity
# so a good channel

# DEM (Direct Email)
# it has 11140 missing values
# It is highly +ve skewed so most of the populations
# are having a good affinity <0.2
# Thus its another good channel


# Region : Doctor office location {Urban, Rural}
# So, we do not have any missing values here :)
# we are having rural and urban reasons
# with more pouplation from urban reasons (>80%)
# As, its a categorical variable, we will analyze it 
# thorugh crosstab for better understanding of the distrubution


# Values of the HCP to Mediplexis 
# (High (H), Medium (M), Low (L), Undefined(U))
# 'H': 23620, 'M': 12430, 'L': 999, 'U': 137
# We have here 137 doctors for which we do not have 
# any values, so they are treated as another category
# It seems, many doctors have high values for Mediplexis


# Doctor's specialty code
# This feature represents the doctors specialization
# in various fields. we have 16 categories in total
# and there are >75% of doctors having speciality
# in category 155443


# Gender
# here we have 28 missing values and we should replace 
# these missing values with some further investigation
# moreover, we have mostly males as doctors 


# Age 
# Here, we have ~25 % population having age > 110
# and we have maximum age as 121 which seems like an
# outlier as its too old age to do a proffession
# we can treat it in multiple ways, either by Winsorization
# or by binning
# We will try both of them and lets see what helps


# So that was little bit of univariate stats related to our features
# Now let's move on
plotHeatMap(df=all_data)
# It seems not so worthy to extract info from above heatmap, so let's move on
# Lets see the missing distribution
all_data.isna().sum().sort_values().plot(kind='bar')
plt.show()
# Speciality wise missing values
for i in list(set(all_data['Speciality_ID'])):
    print('Speciality_ID : ',i)
    all_data[all_data['Speciality_ID']==i].isna().sum().sort_values().plot(kind='bar')
    plt.show()
# Above, it seems that some Speciality_ID have only nan in some channels, but not in all channels, interesting :O<jupyter_output><empty_output><jupyter_text>Preprocessing<jupyter_code># So that was a little bit expolration
# but thinking practically,we can conclude the fact
# that all channels are dependent on each other
# as if a doctor is very much satisfied by one channel
# then he should most probably have low affinity towards 
# other channel.
# so we will create some features that will provide us
# total, max and avg affinity of a doctor from all 8 channels

# One Hot Encoding
from sklearn.preprocessing import OneHotEncoder
# Label encoding
from sklearn.preprocessing import LabelEncoder
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split


# In histplot for Age, distribution seems to have a
# good volume of all aged doctors
# So, lets create a feature on the basis of
# carier experience, It seems setting <50 as low
#  >50 and <70 as Medium
# >70 as High experience would be good as per distribution

def get_experience(x):
    if x<50:
        return 'low'
    elif x>50 and x<70:
        return 'medium'
    else:
        return 'high'



def get_train_test_valid_data(data=None,col = None, test_size=0.2, ignore_cols_to_drop = []):
    data = data.copy()
        
    # simply filling gender with modes as only 28 missing values
    data['gender'] = data['gender'].fillna('MALE')
    
    # fill all nan's with mean other than given col to show its 
    # mean contribution, Note that as we are only using aggregating features,
    # so, mean wouldn't be having that much impact of biasing and as we
    # proceed to ignore some cols during mean imputation because at each step, we will be having true values
    # for one new missing col each time, biasing wouldn't be a problem :)
    for i in [j for j in ['DMS','DRT','RL','RR','DEM','OLA','OLV','P2P'] if j != col and j not in ignore_cols_to_drop]:
        data[i] = data[i].fillna(np.mean(all_data[all_data[i]==all_data[i]][i]))
    
    # feature generation
    # lets create some features using aggregation but column wise :) as all channels are dependent on each other
    # init
    data['total_channel_impact'] = 0
    
    for i in ['DMS','DRT','RL','RR','DEM','OLA','OLV','P2P']:
        if i!=col:
            data['total_channel_impact'] += data[i]
    data['avg_channel_impact'] = data['total_channel_impact']/7
    data['max_channel_impact'] = data[[i for i in ['DMS','DRT','RL','RR','DEM','OLA','OLV','P2P'] if i != col]].max(axis=1)
    
    # getting experience feature
    data['experience'] = data['Age'].apply(lambda x : get_experience(x))
    
    # Now let's normalize the Age feature with std scaler
    # to have 0 mean and unit variance
    data['Age'] = StandardScaler().fit_transform(data['Age'].reshape(-1,1)) 
    
    
    # getting rid of unneccessary dependent features
    data.drop([j for j in ['DMS','DRT','RL','RR','DEM','OLA','OLV','P2P'] if j != col],1,inplace=True)
    
    
    # Now lets one hot encode the categorical values
    cols_for_ohe = ['Region','Value','Speciality_ID','gender','experience']
    
    #Preprocess before OHE
    for i in cols_for_ohe:
        le = LabelEncoder()
        le.fit(data[i])
        data[i] = le.transform(data[i])
        
    
    # print('Started one hot encoding')
    for i in cols_for_ohe:
        # print(i,'\n',set(data[i]))

        ohe=OneHotEncoder(sparse =False)
        ohe.fit(data[i].values.reshape(-1,1))
        tmp=ohe.transform(data[i].values.reshape(-1,1))
        tmp=pd.DataFrame(tmp,dtype=int)
        tmp.columns=[i+'_'+str(j) for j in range(tmp.shape[1])]
        data=data.drop([i],1)
        data=pd.concat([data,tmp],1)        
        # print(data.shape)
    # print('Completed One Hot Encoding')
    
    # getting test data ie all those which are nan's in given col
    test_data = data.loc[(data[col]!=data[col]),[i for i in data.columns]]
    
    # get remaining non nan's data in a given col 
    tmp_data = data.loc[(data[col]==data[col]),[i for i in data.columns]]
    
    # note that, below i am creating train and validation data with given test_size
    # this
    train_data,validation_data,y_tr,y_ts = train_test_split(tmp_data,tmp_data[col],test_size=test_size,random_state=41526)
    
    return train_data,test_data,validation_data
# testing our preprocessor for getting data
train,test,validaton = get_train_test_valid_data(data=all_data,col='RR')
test.head()

train.head()
test.shape,train.shape
train.isnull().sum()
test.isnull().sum()<jupyter_output><empty_output><jupyter_text>Modeling<jupyter_code># Now my modeling strategy is to use the missing values imputation in ascending order of the no of missing items
# column wise, like DMS have least missing values, so we will start with DMS and then ignore its imputation with 
# mean value during new aggregate features generation of next missing column 
import xgboost as xgb
from sklearn.metrics import mean_squared_error<jupyter_output><empty_output><jupyter_text>Refer : Missing Data Imputation Through MachineLearning Algorithms by <NAME>, <NAME>, and <NAME><jupyter_code># As per above paper, iterative approach of filling nan's seems to outperform the classic method of imputation.
# DMS imputation

train,test,validation = get_train_test_valid_data(data=all_data,col='DMS',ignore_cols_to_drop=[])

dtrain = xgb.DMatrix(data=train.drop(['HCP_ID','DMS'],1), label=train['DMS'])
dvalid = xgb.DMatrix(data=validation.drop(['HCP_ID','DMS'],1), label=validation['DMS'])
dtest = xgb.DMatrix(data=test.drop(['HCP_ID','DMS'],1))

params = {}
params['objective'] = 'reg:linear'
params['eta'] = 0.01
params['silent'] = True
params['max_depth'] = 6
params['subsample'] = 0.8
params['colsample_bytree'] = 0.9
params['nthread']=-1
params['eval_metric']='rmse'

watchlist = [(dtrain, 'train'),(dvalid, 'eval')]

model1_1 = xgb.train(params, dtrain, 20000, watchlist, maximize=False, verbose_eval=50, early_stopping_rounds=260)<jupyter_output>C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:65: FutureWarning: reshape is deprecated and will raise in a subsequent release. Please use .values.reshape(...) instead
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\utils\validation.py:475: DataConversionWarning: Data with input dtype int64 was converted to float64 by StandardScaler.
  warnings.warn(msg, DataConversionWarning)<jupyter_text>Ridge regression<jupyter_code>from sklearn.linear_model import Ridge
model=Ridge(alpha=0.025, fit_intercept=False, normalize=True, copy_X=True, 
          max_iter=800, tol=0.01, solver='auto', random_state=64683725)
model.fit(train.drop(['HCP_ID','DMS'],1),train['DMS'])
print(np.power(mean_squared_error(y_true=validation['DMS'],y_pred=model.predict(validation.drop(['HCP_ID','DMS'],1))),.5))<jupyter_output>0.171748018789<jupyter_text>ElasticNet<jupyter_code>from sklearn.linear_model import ElasticNet
model=ElasticNet(alpha=0.00001, fit_intercept=False, normalize=True, 
                   copy_X=True, max_iter=800, tol=0.001, warm_start=False, 
                   random_state=None, selection='cyclic')
model.fit(train.drop(['HCP_ID','DMS'],1),train['DMS'])
print(np.power(mean_squared_error(y_true=validation['DMS'],y_pred=model.predict(validation.drop(['HCP_ID','DMS'],1))),.5))<jupyter_output>0.171733182096<jupyter_text>RandomForest<jupyter_code>from sklearn.ensemble import RandomForestRegressor

model = RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
           max_leaf_nodes=None,
           min_impurity_decrease=0.0, min_impurity_split=None,
           min_samples_leaf=50, min_samples_split=2,
           min_weight_fraction_leaf=0.0, n_estimators=700, n_jobs=-1,
           oob_score=False, random_state=573869792, verbose=0, warm_start=False)

model.fit(train.drop(['HCP_ID','DMS'],1),train['DMS'])
print(np.power(mean_squared_error(y_true=validation['DMS'],y_pred=model.predict(validation.drop(['HCP_ID','DMS'],1))),.5))
#	Plotting feature importance through random forest

def plot_importance(model,train, top_n = None):
	# Classifier model must be pretrained
    plt.figure(figsize=(12,20))
    importance=model.feature_importances_
    features=np.array([str(i) for i in train.columns])
    indices = np.argsort(importance)
    indices = np.array([int(i)  for i in indices[(len(indices)-top_n):]])
    plt.title('Feature Importances')
    plt.barh(range(top_n), importance[indices], color='r', align='center')
    plt.yticks(range(len(indices)), features[indices])
    plt.xlabel('Relative Importance')
    plt.show()
plot_importance(model,train.drop(['HCP_ID','DMS'],1),top_n=25)
plot_importance(model,train.drop(['HCP_ID','DMS'],1),top_n=29)
# In above two feature importance graph, first one without aggreagate features and 
# other one with aggregate feature, shows that our aggragated feature is doing good in predictions
# As random forest seems to have comparatively low rmse than other, so we will use it for others too.
# Defining a imputer function

def get_missing_values(model = None, data = None, col = None, ignore_cols_to_drop = []):
    data = data.copy()
    
    train,test,validation = get_train_test_valid_data(data=all_data,col=col,ignore_cols_to_drop=ignore_cols_to_drop)
    model.fit(train.drop(['HCP_ID',col],1),train[col])
    print(col,'  - > ',np.power(mean_squared_error(y_true=validation[col],y_pred=model.predict(validation.drop(['HCP_ID',col],1))),.5))
    
    test[col] = model.predict(test.drop(['HCP_ID',col],1))
    return test[['HCP_ID',col]]


sub = pd.read_csv('submission.csv')
cols_to_ignore = []

#  Filling missing values in increasing order of no of missing items, columns wise
for col in ['DMS','DRT','RL','RR','DEM','OLA','OLV','P2P']:
    # initializing model
    # Choosing parameters such that there will be many estimators overfitting in 
    # different reason which will most likely capture most of the patterns and then
    # bagging all of them with RandomForest
    model = RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
               max_leaf_nodes=None,
               min_impurity_decrease=0.0, min_impurity_split=None,
               min_samples_leaf=50, min_samples_split=2,
               min_weight_fraction_leaf=0.0, n_estimators=1000, n_jobs=-1,
               oob_score=False, random_state=573869792, verbose=0, warm_start=False)

    result = get_missing_values(model=model,data=all_data,col=col,ignore_cols_to_drop=cols_to_ignore)
    sub.loc[sub['HCP_ID'].isin(result['HCP_ID'].values),col] = result[col].values
    
    # updating the data with new imputation
    cols_to_ignore += [col]
    all_data.loc[all_data['HCP_ID'].isin(result['HCP_ID'].values),col] = result[col].values
    
    
# score mean
np.mean([0.171414250306,0.0944364039064,0.15190042531,0.183543004329,0.0744971943994,0.0810242033781,0.184085600183,0.179187825425])<jupyter_output><empty_output><jupyter_text>Note : was trying to implement below the referred paper iterative approach, but unfortunately time runs out :(  Let's try iterative approach as discussed in reference paperdefine no of iterationno_of_iteration = 5cols_to_ignore = []for i in range(no_of_iteration):        sub = pd.read_csv('submission.csv')        all_data_copy = all_data.copy()                  Filling missing values in increasing order of no of missing items, columns wise    for col in ['DMS','DRT','RL','RR','DEM','OLA','OLV','P2P']:                if i>0:             reseting this column again and use transfer learning in iterative way             note that cols_to_ignore are still same            all_data[col] = all_data_copy[col]                 initializing model         Choosing parameters such that there will be many estimators overfitting in          different reason which will most likely capture most of the patterns and then         bagging all of them with RandomForest        model = RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,                   max_leaf_nodes=None,                   min_impurity_decrease=0.0, min_impurity_split=None,                   min_samples_leaf=50, min_samples_split=2,                   min_weight_fraction_leaf=0.0, n_estimators=700, n_jobs=-1,                   oob_score=False, random_state=573869792, verbose=0, warm_start=False)        result = get_missing_values(model=model,data=all_data,col=col,ignore_cols_to_drop=cols_to_ignore)        sub.loc[sub['HCP_ID'].isin(result['HCP_ID'].values),col] = result[col].values         updating the data with new imputation        cols_to_ignore += [col]        all_data.loc[all_data['HCP_ID'].isin(result['HCP_ID'].values),col] = result[col].values<jupyter_code># Sanity checkup
sub.isnull().sum()
# Thus all missing values are imputed
# Saving our submission file
sub.to_csv('Mayank_Kumar_06101996.csv',index=False)<jupyter_output><empty_output><jupyter_text>Error Analysis<jupyter_code># Let's compare the resulted distribution with previous
old_data = pd.read_csv('submission.csv')
new_data = pd.read_csv('Mayank_Kumar_06101996.csv')

for i in ['DMS','DRT','RL','RR','DEM','OLA','OLV','P2P']:
    print(i)
    describe(data=old_data,x=i)
    describe(data=new_data,x=i)<jupyter_output>DMS
null values = 4379
is_unique : No
Number of unique values :  10624
max = 1.0  min = 0.0
non zero values = 30178   zero values = 7008
count    37186.000000
mean         0.118864
std          0.163822
min          0.000000
25%          0.015625
50%          0.075188
75%          0.127780
max          1.000000
Name: DMS, dtype: float64<jupyter_text>Result<jupyter_code># From the old and new distributions , It appears that previous and new 
# mean and std are still stable comparatively, so our model is doing good job :)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1528166443606278, 'fraction_numerical': 0.02757837530924281, 'mean_word_length': 4.013826758845059, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 63, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '439707', 'license': 'MIT', 'n_tokens_mistral': 8245, 'n_tokens_neox': 7555, 'n_words': 2132, 'path': 'ZS-Young-Data-scientist-2018-Winner/ZS offline phase II/Mayank_Kumar_06101996.ipynb', 'repo_name': 'Lilasrinivasreddy/NLP_Speech_recog', 'size': '24691'}