<jupyter_start><jupyter_text>Problem Statement Customer Lifetime Value represents a customer’s value to a company over a period of time. It’s a competitive market for insurance companies in 2019, and insurance premium isn’t the only determining factor in a customer’s decisions. CLV is a customer-centric metric, and a powerful base to build upon to retain valuable customers, increase revenue from less valuable customers, and improve the customer experience overall. Auto Insurance company is facing issues in retaining its customers and wants to advertise promotional offers for its loyal customers. They are considering CLV as a parameter to classify loyal customers.  Import libraries<jupyter_code>import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import scipy.stats as stats
df=pd.read_csv("Data/AutoInsurance.csv")
df.head()
print("The number of rows:",df.shape[0])
print("The number of columns:",df.shape[1])
df.rename(columns={'Customer Lifetime Value':'CLV'},inplace=True)
df.info()
numerical_cols = df.select_dtypes(include=["int64","float64"])
numerical_cols.columns
numerical_cols = numerical_cols.drop(["Number of Policies","Number of Open Complaints"],axis=1)<jupyter_output><empty_output><jupyter_text>Summary Statistics of Numerical features<jupyter_code>numerical_cols.describe()<jupyter_output><empty_output><jupyter_text>Checking Null values<jupyter_code>df.isnull().sum()<jupyter_output><empty_output><jupyter_text>There are no null values in the data  EDA  Numerical Features  CLV - Target variable analysis<jupyter_code>sns.distplot(df["CLV"])
# plt.savefig("CLV.png")
plt.show()
sns.boxplot(df["CLV"])
plt.show()
df["CLV"].skew()<jupyter_output><empty_output><jupyter_text>From the distribution and box plot we can infer that CLV is highly right skewed. There are lot of outliers in our dependent variable but we cant treat them directly because they are influential points.  Income **Univariate Analysis**<jupyter_code>sns.distplot(df["Income"])
plt.show()
sns.boxplot(df["Income"])
plt.show()<jupyter_output><empty_output><jupyter_text>**Bivariate Analysis**<jupyter_code>sns.scatterplot(df["Income"],df["CLV"])
plt.show()<jupyter_output><empty_output><jupyter_text>There is no linear relationship between income and CLV  Monthly Premium Auto<jupyter_code>sns.distplot(df["Monthly Premium Auto"],color='c')
plt.show()
sns.boxplot(df["Monthly Premium Auto"],color='c')
plt.show()<jupyter_output><empty_output><jupyter_text>**Bivariate Analysis**<jupyter_code>sns.scatterplot(df["Monthly Premium Auto"],df["CLV"],color='g')
# plt.savefig("bi.png")
plt.show()<jupyter_output><empty_output><jupyter_text>There is linear relationship between income and CLV  Months Since Last Claim<jupyter_code>sns.distplot(df['Months Since Last Claim'])
plt.show()
sns.boxplot(df["Months Since Last Claim"])
plt.show()<jupyter_output><empty_output><jupyter_text>**Bivariate Analysis**<jupyter_code>sns.scatterplot(df["Months Since Last Claim"],df["CLV"],color='c')
plt.show()<jupyter_output><empty_output><jupyter_text>There is no linear relationship  Months Since Policy Inception<jupyter_code>sns.distplot(df['Months Since Policy Inception'])
plt.show()
sns.boxplot(df["Months Since Policy Inception"])
plt.show()<jupyter_output><empty_output><jupyter_text>Total claim amount<jupyter_code>sns.distplot(df['Total Claim Amount'])
plt.show()
sns.boxplot(df["Total Claim Amount"])
plt.show()
sns.scatterplot(df["Total Claim Amount"],df["CLV"],color='c')
plt.show()<jupyter_output><empty_output><jupyter_text>There is a linear relationship between CLV and Total Claim Amount The monthly premium auto and income feature has multiple peak's so to remove those peaks we can apply any of the power transformation (SQUARE / CUBE)<jupyter_code>sns.distplot(np.square(df["Monthly Premium Auto"]),color='r')
plt.show()
sns.distplot(np.square(df["Income"]),color='r')
plt.show()<jupyter_output><empty_output><jupyter_text>But as we can see that after the square transformation the data is getting heavily skewed, so we stick with the actual distribution again. We will just stick with the actual distribution of the feature in our model building.  Categorical variables<jupyter_code>cat_cols = df.select_dtypes(include="object")
no_col = df[["Number of Open Complaints","Number of Policies"]]
cat_cols = pd.concat([cat_cols,no_col],axis=1)
cat_cols.head()
cat_cols.drop("Effective To Date",axis=1,inplace=True)
cat_cols.columns
for i in cat_cols:
    print("Unique values in",str(i),"is",df[i].nunique())
    print(df[i].value_counts())
    print("-----------------------------------------")<jupyter_output>Unique values in Customer is 9134
AP52565    1
YU72630    1
EV66951    1
RF58958    1
KG83280    1
          ..
LY17372    1
DQ93036    1
PS57128    1
KE70943    1
BK21146    1
Name: Customer, Length: 9134, dtype: int64
-----------------------------------------
Unique values in State is 5
California    3150
Oregon        2601
Arizona       1703
Nevada         882
Washington     798
Name: State, dtype: int64
-----------------------------------------
Unique values in Response is 2
No     7826
Yes    1308
Name: Response, dtype: int64
-----------------------------------------
Unique values in Coverage is 3
Basic       5568
Extended    2742
Premium      824
Name: Coverage, dtype: int64
-----------------------------------------
Unique values in Education is 5
Bachelor                2748
College                 2681
High School or Below    2622
Master                   741
Doctor                   342
Name: Education, dtype: int64
-----------------------------------------
Unique values in Em[...]<jupyter_text>Customer column which is the id has unique values for all records so it can be removed for model building as it adds no significance  State<jupyter_code>sns.countplot(df["State"])
plt.show()<jupyter_output><empty_output><jupyter_text>Most number of people are residing in california  Response<jupyter_code>sns.countplot(df["Response"])
plt.show()<jupyter_output><empty_output><jupyter_text>Its seems that less than 1400 customers have responded to the marketing calls  Coverage<jupyter_code>sns.countplot(df["Coverage"])
plt.show()<jupyter_output><empty_output><jupyter_text>Education<jupyter_code>plt.figure(figsize=(10,5))
sns.countplot(df["Education"])
plt.show()<jupyter_output><empty_output><jupyter_text>Most of the customers have either a bachelor or college degree  EmploymentStatus<jupyter_code>plt.figure(figsize=(10,5))
sns.countplot(df["EmploymentStatus"])
# plt.savefig("empstatus.png")
plt.show()
df["EmploymentStatus"].value_counts(normalize=True)*100<jupyter_output><empty_output><jupyter_text>Around 62.38% of the customers are employed<jupyter_code>plt.figure(figsize=(10,5))
sns.boxplot(x=df['EmploymentStatus'],y=df['CLV'])
plt.title("CLV vs EmploymentStatus")
# plt.savefig("bi1.png")
plt.show()<jupyter_output><empty_output><jupyter_text>Gender<jupyter_code>sns.countplot(df["Gender"])
plt.show()<jupyter_output><empty_output><jupyter_text>Location Code<jupyter_code>sns.countplot(df["Location Code"])
plt.savefig("location.png")
plt.show()<jupyter_output><empty_output><jupyter_text>Marital Status<jupyter_code>sns.countplot(df["Marital Status"])
plt.show()

sns.barplot(df["Marital Status"],df["CLV"])
plt.show()<jupyter_output><empty_output><jupyter_text>Policy Type<jupyter_code>sns.countplot(df["Policy Type"])
plt.show()<jupyter_output><empty_output><jupyter_text>Policy<jupyter_code>plt.figure(figsize=(15,5))
sns.countplot(df["Policy"])
plt.show()<jupyter_output><empty_output><jupyter_text>Personal L3 policy subcategory has the most number of customers  Renew Offer type<jupyter_code>sns.countplot(df["Renew Offer Type"])
plt.show()<jupyter_output><empty_output><jupyter_text>Sales Channel<jupyter_code>sns.countplot(df["Sales Channel"])
plt.show()<jupyter_output><empty_output><jupyter_text>Vehicle Class<jupyter_code>plt.figure(figsize=(10,5))
sns.countplot(df["Vehicle Class"])
plt.show()<jupyter_output><empty_output><jupyter_text>Vehicle Size<jupyter_code>plt.figure(figsize=(10,5))
sns.countplot(df["Vehicle Size"])
plt.show()<jupyter_output><empty_output><jupyter_text>Number of Open Complaints<jupyter_code>plt.figure(figsize=(10,5))
sns.countplot(df["Number of Open Complaints"])
plt.show()
plt.figure(figsize=(10,5))
sns.countplot(df["Number of Policies"])
plt.show()
df['Effective To Date']=pd.to_datetime(df['Effective To Date'],infer_datetime_format=True)
df["Months"] = df["Effective To Date"].dt.month
df['Months'] = df['Months'].astype('object')<jupyter_output><empty_output><jupyter_text>Months<jupyter_code>plt.figure(figsize=(10,5))
sns.countplot(df["Months"])
plt.show()
cat_cols.columns
df.columns<jupyter_output><empty_output><jupyter_text>HeatMap<jupyter_code>heatmap = df[['CLV', 'Income', 'Monthly Premium Auto','Total Claim Amount']]
plt.figure(figsize=(7,4))
sns.heatmap(heatmap.corr(),annot=True)
# plt.savefig("Heatmap.png")
plt.show()<jupyter_output><empty_output><jupyter_text>We can clearly see in the heatmap, that customer lifetime value has a better correlation with monthly premium auto and acceptable correlation with total claim amount.  Statistical Significance - H0:The null hypothesis of the test is the data is normally distributed- HA:The data is not normally distributed<jupyter_code>stats.shapiro(df['CLV'])<jupyter_output>C:\Users\Nantha\Anaconda3\lib\site-packages\scipy\stats\morestats.py:1660: UserWarning: p-value may not be accurate for N > 5000.
  warnings.warn("p-value may not be accurate for N > 5000.")<jupyter_text>p-value less than 0.05 we reject null hypothesis the data is not normally distributed We will proceed with non parametric tests since the dependent variable is not normally distributed<jupyter_code>columns_cat = list(cat_cols.columns)
columns_cat
manwhiteneyy = []
anova = []

for i in columns_cat:
    if (df[i].nunique()>2):
        anova.append(i)
    else:
        manwhiteneyy.append(i)
print("Anova:",anova)
print("TTest:",manwhiteneyy)<jupyter_output>Anova: ['Customer', 'State', 'Coverage', 'Education', 'EmploymentStatus', 'Location Code', 'Marital Status', 'Policy Type', 'Policy', 'Renew Offer Type', 'Sales Channel', 'Vehicle Class', 'Vehicle Size', 'Number of Open Complaints', 'Number of Policies']
TTest: ['Response', 'Gender']<jupyter_text>H0, the distributions of both populations are equalThe alternative hypothesis H1 is that the distributions are not equal.<jupyter_code>yes=df[df['Response']=='No']['CLV']
no=df[df['Response']=='Yes']['CLV']
sns.distplot(yes)
sns.distplot(no)<jupyter_output><empty_output><jupyter_text>If the distributions are identical, which is the null hypothesis of the Mann-Whitney U test, the mean rank will be the same<jupyter_code>f = stats.mannwhitneyu(yes,no)
print(f)<jupyter_output>MannwhitneyuResult(statistic=5060988.0, pvalue=0.2584436122645599)<jupyter_text>p-value is greater than 0.05 we failed to reject null hypothesis i.e the mean rank will be the same for both males and females.<jupyter_code>male=df[df['Gender']=='M']['CLV']
female=df[df['Gender']=='F']['CLV']
b =stats.mannwhitneyu(male,female)
print(b)<jupyter_output>MannwhitneyuResult(statistic=10337134.0, pvalue=0.24374651528224917)<jupyter_text>p-value is greater than 0.05 we failed to reject null hypothesis i.e the mean rank will be the same for both males and females. we test above two features with non-Paramteric unpaired ttest and we concluded that those 2 features is not significant to predict traget - H0:the mean ranks of the groups are the same.- H1:the mean ranks of the groups are not the same.<jupyter_code>ca=df[df['State']=='California']['CLV']
Or=df[df['State']=='Oregon']['CLV']
Ar=df[df['State']=='Arizona']['CLV']
Ne=df[df['State']=='Nevada']['CLV']
Wa=df[df['State']=='Washington']['CLV']
a = stats.kruskal(ca,Or,Ar,Ne,Wa)
print(a)<jupyter_output>KruskalResult(statistic=5.072069583836472, pvalue=0.27998232358818953)<jupyter_text>p-value > 0.05 failed to reject H0<jupyter_code>stats.f_oneway(ca,Or,Ar,Ne,Wa)
print(b)
# df.rename(columns={'CLV':'CLV'},inplace=True)
Ba=df[df['Coverage']=='Basic']['CLV']
Ex=df[df['Coverage']=='Extended']['CLV']
Pr=df[df['Coverage']=='Premium']['CLV']
sns.distplot(Ba,hist=False,label='Basic')
sns.distplot(Ex,hist=False)
sns.distplot(Pr,hist=False)
plt.show()
c =stats.kruskal(Ba,Ex,Pr)
print(c)<jupyter_output>KruskalResult(statistic=502.4996090762843, pvalue=7.648852917746348e-110)<jupyter_text>p-value < 0.05 to reject H0<jupyter_code>stats.f_oneway(Ba,Ex,Pr)
Ba=df[df['Education']=='Bachelor']['CLV']
Co=df[df['Education']=='College']['CLV']
Hi=df[df['Education']=='High School or Below']['CLV']
Ma=df[df['Education']=='Master']['CLV']
Da=df[df['Education']=='Doctor']['CLV']
d= stats.kruskal(Ba,Co,Hi,Ma,Da)
print(d)
plt.figure(figsize=(15,5))
sns.distplot(Ba,hist=False,label='Bachelor')
sns.distplot(Co,hist=False,label='College')
sns.distplot(Hi,hist=False,label='High School or Below')
sns.distplot(Ma,hist=False,label='Master')
sns.distplot(Da,hist=False,label='Doctor')
plt.title("Different Education level comparsion with CLV")
# plt.savefig("Education.png")
plt.show()<jupyter_output><empty_output><jupyter_text>p-value < 0.05 reject H0, Education is significant in predicting CLV<jupyter_code>Ba=df[df['EmploymentStatus']=='Employed']['CLV']
Co=df[df['EmploymentStatus']=='Unemployed']['CLV']
Hi=df[df['EmploymentStatus']=='Medical Leave']['CLV']
Ma=df[df['EmploymentStatus']=='Disabled']['CLV']
Da=df[df['EmploymentStatus']=='Retired']['CLV']
e = stats.kruskal(Ba,Co,Hi,Ma,Da)
print(e)<jupyter_output>KruskalResult(statistic=42.56159302970864, pvalue=1.2758458797542519e-08)<jupyter_text>p-value < 0.05 reject H0, EmploymentStatus is significant in predicting CLV<jupyter_code>Ba=df[df['Location Code']=='Suburban']['CLV']
Co=df[df['Location Code']=='Rural']['CLV']
Hi=df[df['Location Code']=='Urban']['CLV']
g = stats.kruskal(Ba,Co,Hi)
print(g)<jupyter_output>KruskalResult(statistic=2.4638212999332567, pvalue=0.29173464205664124)<jupyter_text>p-value < 0.05 reject H0, Location Code can be significant in predicting CLV<jupyter_code>Ba=df[df['Marital Status']=='Married']['CLV']
Co=df[df['Marital Status']=='Single']['CLV']
Hi=df[df['Marital Status']=='Divorced']['CLV']
h = stats.kruskal(Ba,Co,Hi)
print(h)<jupyter_output>KruskalResult(statistic=20.895798766808223, pvalue=2.9009146434700455e-05)<jupyter_text>p-value < 0.05 reject H0, Marital Status can be significant in predicting CLV<jupyter_code>Ba=df[df['Number of Open Complaints']==0]['CLV']
Co=df[df['Number of Open Complaints']==1]['CLV']
Hi=df[df['Number of Open Complaints']==2]['CLV']
Ma=df[df['Number of Open Complaints']==3]['CLV']
Da=df[df['Number of Open Complaints']==4]['CLV']
As=df[df['Number of Open Complaints']==5]['CLV']
o = stats.kruskal(Ba,Co,Hi,Ma,Da,As)
print(o)<jupyter_output>KruskalResult(statistic=29.160617988301833, pvalue=2.1564409504417123e-05)<jupyter_text>p-value < 0.05 reject H0, Number of Open Complaints can be significant in predicting CLV<jupyter_code>Ba=df[df['Number of Policies']==1]['CLV']
Co=df[df['Number of Policies']==2]['CLV']
Hi=df[df['Number of Policies']==3]['CLV']
p = stats.kruskal(Ba,Co,Hi)
print(p)<jupyter_output>KruskalResult(statistic=5017.7060662040785, pvalue=0.0)<jupyter_text>p-value < 0.05 reject H0, Number of Policies can be significant in predicting CLV<jupyter_code>Ba=df[df['Policy Type']=='Personal Auto']['CLV']
Co=df[df['Policy Type']=='Corporate Auto']['CLV']
Hi=df[df['Policy Type']=='Special Auto']['CLV']
i = stats.kruskal(Ba,Co,Hi)
print(i)<jupyter_output>KruskalResult(statistic=4.6074986354931795, pvalue=0.0998836452694168)<jupyter_text>p-value > 0.05 failed to reject H0, Policy Type may not be significant in predicting CLV<jupyter_code>Ba=df[df['Policy']=='Personal L3']['CLV']
Co=df[df['Policy']=='Personal L2']['CLV']
Hi=df[df['Policy']=='Personal L1']['CLV']
Ma=df[df['Policy']=='Corporate L3']['CLV']
Da=df[df['Policy']=='Corporate L2']['CLV']
Ca=df[df['Policy']=='Corporate L1']['CLV']
s1=df[df['Policy']=='Special L2']['CLV']
s2=df[df['Policy']=='Special L3']['CLV']
s3=df[df['Policy']=='Special L1']['CLV']
j =stats.kruskal(Ba,Co,Hi,Ma,Da,Ca,s1,s2,s3)
print(j)<jupyter_output>KruskalResult(statistic=7.944351884462466, pvalue=0.4389248352544014)<jupyter_text>p-value > 0.05 failed to reject H0, Policy may not be significant in predicting CLV<jupyter_code>Ba=df[df['Renew Offer Type']=='Offer1']['CLV']
Co=df[df['Renew Offer Type']=='Offer2']['CLV']
Hi=df[df['Renew Offer Type']=='Offer3']['CLV']
Ma=df[df['Renew Offer Type']=='Offer4']['CLV']
k =stats.kruskal(Ba,Co,Hi,Ma)
print(k)<jupyter_output>KruskalResult(statistic=168.90219457792458, pvalue=2.196144881912684e-36)<jupyter_text>p-value < 0.05 reject H0, Number of Policies can be significant in predicting CLV<jupyter_code>Ba=df[df['Sales Channel']=='Agent']['CLV']
Co=df[df['Sales Channel']=='Branch']['CLV']
Hi=df[df['Sales Channel']=='Call Center']['CLV']
Ma=df[df['Sales Channel']=='Web']['CLV']
l = stats.kruskal(Ba,Co,Hi,Ma)
print(l)<jupyter_output>KruskalResult(statistic=4.491751501229693, pvalue=0.2130272157209903)<jupyter_text>p-value > 0.05 failed to reject H0, Sales Channel may not be significant in predicting CLV<jupyter_code>Ba=df[df['Vehicle Class']=='Four-Door Car']['CLV']
Co=df[df['Vehicle Class']=='Two-Door Car']['CLV']
Hi=df[df['Vehicle Class']=='SUV']['CLV']
Ma=df[df['Vehicle Class']=='Sports Car']['CLV']
Da=df[df['Vehicle Class']=='Luxury SUV']['CLV']
Ca=df[df['Vehicle Class']=='Luxury Car']['CLV']
m = stats.kruskal(Ba,Co,Hi,Ma,Da,Ca)
print(m)<jupyter_output>KruskalResult(statistic=1310.5494960284416, pvalue=3.309676344717952e-281)<jupyter_text>p-value < 0.05 reject H0, Vehicle Class can be significant in predicting CLV<jupyter_code>Ba=df[df['Vehicle Size']=='Medsize']['CLV']
Co=df[df['Vehicle Size']=='Small']['CLV']
Hi=df[df['Vehicle Size']=='Large']['CLV']
n =stats.kruskal(Ba,Co,Hi)
print(n)<jupyter_output>KruskalResult(statistic=9.564995661218877, pvalue=0.00837505335446059)<jupyter_text>p-value < 0.05 reject H0, Vehicle Class can be significant in predicting CLV<jupyter_code>columns_cat
cat_cols.drop("Customer",axis=1,inplace=True)
cat_cols
p_values = [a[1],b[1],c[1],d[1],e[1],f[1],g[1],h[1],i[1],j[1],k[1],l[1],m[1],n[1],o[1],p[1]]
p_values
len(p_values)
len(cat_cols.columns)
pd.DataFrame(p_values,index=(cat_cols.columns),columns=['P_Value'])<jupyter_output><empty_output><jupyter_text>Coverage,Education,EmploymentStatus,Marital Status, Renew Offer Type, Vehicle class, Vehicle class, Number of open complaints and No of Policies play a significant role in predicting the CLV<jupyter_code>numerical_cols.head()
cat_cols.head()
import statsmodels.api as sm
# X=numerical_cols.drop(['CLV'],axis=1)
# y=numerical_cols['CLV']
# lin_reg=sm.OLS(y,X).fit()
# lin_reg.summary()<jupyter_output><empty_output><jupyter_text>Data Preprocessing<jupyter_code>cat_cols.head()
catg=pd.get_dummies(cat_cols,drop_first=True)
numerical_cols.head()
dfn = pd.concat([numerical_cols,catg],axis=1)
dfn.head()<jupyter_output><empty_output><jupyter_text>OLS Model<jupyter_code>dfn.rename(columns={'CLV':'CLV'},inplace=True)
X=dfn.drop(['CLV'],axis=1)
y=dfn['CLV']
X_constant=sm.add_constant(X)
lin_reg=sm.OLS(y,X_constant).fit()
lin_reg.summary()<jupyter_output>C:\Users\Nantha\Anaconda3\lib\site-packages\numpy\core\fromnumeric.py:2389: FutureWarning: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
  return ptp(axis=axis, out=out, **kwargs)<jupyter_text>Assumptions1. No autocorrelation2. Linearity of residuals3. Normality of error terms4. No Heteroscedacity5. No strong MultiCollinearity **Assumption 1 No Auto correlation.**Test needed : Durbin- Watson Test.It's value ranges from 0-4. If the value of Durbin- Watson is Between 0-2, it's known as Positive Autocorrelation.If the value ranges from 2-4, it is known as Negative autocorrelation.If the value is exactly 2, it means No Autocorrelation.For a good linear model, it should have low or no autocorrelation.we can see here the values of dublin watson test: 1.995 (No AUTO-CORRELATION) From the graph below, we can easily see that there is no autocorrelation.<jupyter_code>import statsmodels.tsa.api as smt #timeseries analysis

acf = smt.graphics.plot_acf(lin_reg.resid, lags=40 , alpha=0.05)#auto_correlation value, lags= previous value
acf.show()<jupyter_output>C:\Users\Nantha\Anaconda3\lib\site-packages\ipykernel_launcher.py:4: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.
  after removing the cwd from sys.path.<jupyter_text>Assumption 2- Normality of Residuals2) The second assumption is the Normality of Residuals / Error terms.For this we prefer the Jarque Bera test. For a good model, the residuals should be normally distributed. The higher the value of Jarque Bera test , the lesser the residuals are normally distributed. We generally prefer a lower value of jarque bera test.The Jarque–Bera test is a goodness-of-fit test of whether sample data have the skewness and kurtosis matching a normal distribution. A large value for the jarque-bera test indicates non normality.The jarque bera test tests whether the sample data has the skewness and kurtosis matching a normal distribution. Note that this test generally works good for large enough number of data samples(>2000) as the test statistics asymptotically has a chi squared distribution with degrees 2 of freedom.<jupyter_code>#H0: the error terms are normally distributed
#dof for residuals is n-2
from scipy import stats
print(stats.jarque_bera(lin_reg.resid))<jupyter_output>(65051.54524138234, 0.0)<jupyter_text>The critical chi square value at the 5% level of significance is 5.99. If the computed value exceeds this value the null hypothesis is rejected.In this case the computed value of the JB statistic 65051.11 is greater than 5.99. Thus we reject the null hypothesis that the error terms are normally distributed.<jupyter_code>import seaborn as sns

sns.distplot(lin_reg.resid)
plt.show()<jupyter_output><empty_output><jupyter_text>From the graph we infer that due to the high value of Jarque Bera test, only few residuals are normally distributed **Asssumption 3 - Linearity of residuals**Here we have 2 options. Either we can plot the observed values Vs predicted values and plot the Residual Vs predicted values and see the linearity of residuals. OR We can go for rainbow test. Let's look both of them one by one.<jupyter_code>import statsmodels.api as sm
sm.stats.diagnostic.linear_rainbow(res=lin_reg, frac=0.5)#the fraction of the data to include in the center model<jupyter_output><empty_output><jupyter_text>H0:The data is linear, pvalue > 0.05 failed to reject H0, here we failed to reject H0 so residuals are normally distributed<jupyter_code>lin_reg.resid.mean()<jupyter_output><empty_output><jupyter_text>Assumption 4 - Homoscedasticity_test(using goldfeld test) OR (Beusch-Wagon Test)Homoscedacity :: If the variance of the residuals are symmetrically distributed across the regression line , then the data is said to homoscedastic.Heteroscedacity :: If the variance is unequal for the residuals across the regression line, then the data is said to be heteroscedastic. In this case the residuals can form an arrow shape or any other non symmetrical shape.This test is based on the hytpothesis testing where null and alternate hypothesis are:𝐻0:𝜎𝑢𝑖 𝑖𝑠 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 𝑎𝑐𝑟𝑜𝑠𝑠 𝑡ℎ𝑒 𝑟𝑎𝑛𝑔𝑒 𝑜𝑓 𝑑𝑎𝑡𝑎𝐻𝑎:𝜎𝑢𝑖 𝑖𝑠 𝑛𝑜𝑡 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 𝑎𝑐𝑟𝑜𝑠𝑠 𝑡ℎ𝑒 𝑟𝑎𝑛𝑔𝑒 𝑜𝑓 𝑑𝑎𝑡𝑎The residuals should be homoscedastic.Goldfeld Test<jupyter_code>import statsmodels.stats.api as sms
name=['F-statistic','p-value']
test=sms.het_goldfeldquandt(lin_reg.resid,lin_reg.model.exog)

test<jupyter_output><empty_output><jupyter_text>p-value below a certain level (like 0.05) indicates we should reject the null in favor of heteroscedasticity. **Assumption 5- NO MULTI COLLINEARITY** The interpretation of VIF is as follows: the square root of a given variable’s VIF shows how much larger the standard error is, compared with what it would be if that predictor were uncorrelated with the other features in the model. If no features are correlated, then all values for VIF will be 1.<jupyter_code>from statsmodels.stats.outliers_influence import variance_inflation_factor
vif = [variance_inflation_factor(X_constant.values, i) for i in range(X_constant.shape[1])]
pd.DataFrame({'vif': vif[1:]}, index=X.columns)<jupyter_output>C:\Users\Nantha\Anaconda3\lib\site-packages\statsmodels\stats\outliers_influence.py:185: RuntimeWarning: divide by zero encountered in double_scalars
  vif = 1. / (1. - r_squared_i)<jupyter_text>we can see that there is high VIF in column Policy & policy typeAlso we have seen in statistical approach these veriables are not significant to predict targethence from above 2 conclusion we can remove them<jupyter_code>from sklearn.model_selection import train_test_split
X=dfn.drop(['CLV','Policy Type_Personal Auto','Policy Type_Special Auto','Policy_Personal L1','Policy_Personal L2','Policy_Personal L3',
'Policy_Special L1','Policy_Special L2','Policy_Special L3'],axis=1)
y=dfn['CLV']          
y =np.log(y)
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.30,random_state=42)
numerical_cols.rename(columns={'CLV':'CLV'},inplace=True)
print(X_train.shape)
print(y_train.shape)
print(X_test.shape)
print(y_test.shape)<jupyter_output>(6393, 42)
(6393,)
(2741, 42)
(2741,)<jupyter_text>Model Building  Linear Regression<jupyter_code>from sklearn.linear_model import LinearRegression
lr = LinearRegression()
model = lr.fit(X_train,y_train)
print(f'R^2 score for train: {lr.score(X_train, y_train)}')
print(f'R^2 score for test: {lr.score(X_test, y_test)}')
y_pred = model.predict(X_test)
from sklearn.metrics import mean_squared_error,r2_score,mean_absolute_error
print('RMSE:',np.sqrt(mean_squared_error(y_test,y_pred)))
print("MAE",mean_absolute_error(y_test,y_pred))
print('R-squared:',r2_score(y_test,y_pred)) 
#Feature Selection by Recursive Backward Elimination
cols = list(X.columns)
pmax = 1
while (len(cols)>0):
    p= []
    X = X[cols]
    Xc = sm.add_constant(X)
    model = sm.OLS(y,Xc).fit()
    p = pd.Series(model.pvalues.values[1:],index = cols)      
    pmax = max(p)
    feature_with_p_max = p.idxmax()
    if(pmax>0.05):
        cols.remove(feature_with_p_max)
    else:
        break
selected_features = cols
print(selected_features)
X_new = X[['Monthly Premium Auto', 'Number of Open Complaints', 'Number of Policies', 'Coverage_Extended', 'Coverage_Premium', 'Education_High School or Below', 'EmploymentStatus_Employed', 'Marital Status_Single', 'Policy_Corporate L2', 'Renew Offer Type_Offer2', 'Renew Offer Type_Offer3', 'Renew Offer Type_Offer4', 'Vehicle Class_SUV', 'Vehicle Class_Sports Car']]
X_constant_new=sm.add_constant(X_new)
lin_reg=sm.OLS(y,X_constant_new).fit()
lin_reg.summary()<jupyter_output><empty_output><jupyter_text>Rigde and Lasso Regression<jupyter_code>from sklearn.linear_model import Ridge, Lasso
ridge=Ridge(alpha=0.1,normalize=True)
ridge.fit(X_train,y_train)
print('Ridge Regression')
y_pred = ridge.predict(X_test)
print('RMSE:', np.sqrt(mean_squared_error(y_test, y_pred)))
print('R-squared:',r2_score(y_test,y_pred)) 
ridge=Ridge(alpha=1.0,normalize=True)
ridge.fit(X_train,y_train)
print('Ridge Regression')
y_pred = ridge.predict(X_test)
print('RMSE:', np.sqrt(mean_squared_error(y_test, y_pred)))
print('R-squared:',r2_score(y_test,y_pred)) 
lasso = Lasso(random_state=1)
lasso.fit(X_train,y_train)
print('Ridge Regression')
y_pred = lasso.predict(X_test)
print('RMSE:', np.sqrt(mean_squared_error(y_test, y_pred)))
print('R-squared:',r2_score(y_test,y_pred))<jupyter_output>Ridge Regression
RMSE: 0.5992937645386806
R-squared: 0.19679167777854611<jupyter_text>Decision Tree<jupyter_code>from sklearn.tree import DecisionTreeRegressor
dt = DecisionTreeRegressor(random_state=1)

dt.fit(X_train,y_train)
y_pred=dt.predict(X_test)
print('RMSE:',np.sqrt(mean_squared_error(y_test,y_pred)))
print("MAE",mean_absolute_error(y_test,y_pred))
print('R-squared:',r2_score(y_test,y_pred)) 
from sklearn.model_selection import train_test_split
X=dfn.drop(['CLV','Policy Type_Personal Auto','Policy Type_Special Auto','Policy_Personal L1','Policy_Personal L2','Policy_Personal L3',
'Policy_Special L1','Policy_Special L2','Policy_Special L3'],axis=1)
y=np.log(dfn['CLV'])         
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.30,random_state=42)<jupyter_output><empty_output><jupyter_text>RandomForest<jupyter_code>from sklearn.ensemble import RandomForestRegressor
rf = RandomForestRegressor(random_state=1)

rf.fit(X_train,y_train)

y_pred=rf.predict(X_test)
print('RMSE:',np.sqrt(mean_squared_error(y_test,y_pred)))
print("MAE",mean_absolute_error(y_test,y_pred))
print('R-squared:',r2_score(y_test,y_pred))<jupyter_output>C:\Users\Nantha\Anaconda3\lib\site-packages\sklearn\ensemble\forest.py:245: FutureWarning: The default value of n_estimators will change from 10 in version 0.20 to 100 in 0.22.
  "10 in version 0.20 to 100 in 0.22.", FutureWarning)<jupyter_text>Hyperparameter tuning of random forest<jupyter_code>from sklearn.model_selection import GridSearchCV
rf = RandomForestRegressor()
params = {
        'max_depth' : [10,20,30],
        
        'n_estimators' : [100,200,50],
        
        "bootstrap" : [True, False],
    
        'max_features': ['auto', 'sqrt', 'log2']
        
        }

grid = GridSearchCV(estimator = rf, param_grid=params, cv = 5, n_jobs = -1, return_train_score = True )
grid.fit(X,y)
grid.best_params_
rf = RandomForestRegressor(**grid.best_params_)

rf.fit(X_train,y_train)

y_pred=rf.predict(X_test)
print('RMSE:',np.sqrt(mean_squared_error(y_test,y_pred)))
print("MAE",mean_absolute_error(y_test,y_pred))
print('R-squared:',r2_score(y_test,y_pred))
from sklearn.ensemble import AdaBoostRegressor
adaboost = AdaBoostRegressor(
    RandomForestRegressor(max_depth=5,criterion='mse'),
    n_estimators=150
)
adaboost.fit(X_train, y_train)
y_pred = adaboost.predict(X_test)
print('RMSE:',np.sqrt(mean_squared_error(y_test,y_pred)))
print("MAE",mean_absolute_error(y_test,y_pred))
print('R-squared:',r2_score(y_test,y_pred))
pd.DataFrame(rf.feature_importances_,index=X.columns,columns=['Score']).sort_values("Score",ascending=False)[:10]
pd.DataFrame(rf.feature_importances_,index=X.columns,columns=['Score']).sort_values("Score",ascending=False).plot(kind='bar',figsize=(20,7))
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14507296278078374, 'fraction_numerical': 0.03508771929824561, 'mean_word_length': 7.380324264907942, 'pattern_counts': {'":': 0, '<': 258, '<?xml version=': 0, '>': 259, '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': '68', 'dir': 'jupyter-structured-clean-dedup', 'id': '229540', 'license': 'MIT', 'n_tokens_mistral': 12059, 'n_tokens_neox': 11025, 'n_words': 2329, 'path': 'Customer Lifetime Value Prediction.ipynb', 'repo_name': 'nanthasnk/Customer-Lifetime-Value-Prediction', 'size': '30495'}