<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Use the cleaner chaining method for transforming the data https
Step1:  The sale price is in hte hundreds of thousands, so let's divide the price by 1000 to get more manageable numbers.
Step2:  The distribution is skewed (as demonstrated by the large z-score (and small pvalue) of teh skewtest). It is right skewed (the skew is positive). Skewed distribution are not ideal for linear models, which often assume a normal distribution. One way to correct for right-skewness is to take the log [1,2,3]
[1] http
Step4:  Merge the training and test datasets for data preparation
We're going to explore the training dataset and apply some transformations to it (fixing missing values, transforming columns etc). We'll need to apply the same transformations to the test dataset. To make that easy, let's put the training and test datasets into one dataframe.
Step5:  Initialize pipeline with raw data. We can always get the data and apply all the transformations in the pipeline by calling pp().
Step6:  Select Numerical features
The dataset is wide with 78 features. Create dataframe containing the numerical features only.
Step7:  We've got 3 data types
Step8:  Split the data between categorical and numerical features
Step9:  We've got 36 numerical features. We can use the describe method to get some statistics
Step10:  But that's a lot of numbers to digest. Better get started plotting!
Step11:  Deal with NaN values
Step12:  Based on the description, the null values for the MasVnrArea should be 0 (no massonry veneer type)
Step13:  Order columns in alphabetical order
Step14:  Plot violinplots for each feature
The violin plots give us some idea of the distribution of data for each feature. We can look for things like skewness, non-normality, and the presence of outliers.
Step15:  Many of the features are higly skewed and some have very long tails. Some have discrete values (YrSold, Fireplaces).
The features with very long and thin tails, such as ScreenPorch, are almost constant (blobs with long tail) as can be seen below
Step16:  Drop nearly constant features
Step17:  Log transform the other features if they have a high skewness
Using a log transformation for some of the skewed features should help, as illustrated below. We use the raw data (not the standardized one) because we need positive values for the log function (we'll standardize the transformed variables later).
Step20:  Use dataframe & series whenever possible for maximum flexibility (see below)
Step22:  Let's apply a log1p transform to all these and plot the distributions again
Step23:  Now our originally skewed features look more symmetric. 
Check that the distributions are less skewed
Step24:  Save transformed numerical data
Use the storage magic to communicate between notebooks.
Step25:  Correlations
We're now in a good position to identify the key numerical features. Those should be hightly correlated with the sale price.
Step26:  Sort columns in dfnum2 by correlation.
Step27:  Scatter plots
The correlation can be high even if there is no linear relationship between a feature and the target, so we should check the scatter plots.
Step28:  Some features have some sort of bi-modal distribution with a lots of 0 values.
Step29:  Dealing with the zeros
Let's recompute the correlations after dropping the zeros.
Step30:  With the zeros
Step31:  Without the zeros
Step32:  Data imputation
Impute the zero values from the rest of the dataset for the "not ok" columns.
<ASSISTANT_TASK:>
Python Code:
target = pd.read_csv('../data/train_target.csv')
target.describe()
Explanation: Use the cleaner chaining method for transforming the data https://tomaugspurger.github.io/method-chaining.html
Sale price distribution
First step is to look at the target sale price for the training data set, i.e. the column we're trying to predict.
End of explanation
target = target / 1000
sns.distplot(target);
plt.title('SalePrice')
import scipy as sp
sp.stats.skew(target)
sp.stats.skewtest(target)
Explanation: The sale price is in hte hundreds of thousands, so let's divide the price by 1000 to get more manageable numbers.
End of explanation
logtarget = np.log1p(target)
print('skewness of logtarget = ', sp.stats.skew(logtarget)[0])
print('skewness test of logtarget = ', sp.stats.skewtest(logtarget))
sns.distplot(logtarget)
plt.title(r'log(1 + SalePrice)')
Explanation: The distribution is skewed (as demonstrated by the large z-score (and small pvalue) of teh skewtest). It is right skewed (the skew is positive). Skewed distribution are not ideal for linear models, which often assume a normal distribution. One way to correct for right-skewness is to take the log [1,2,3]
[1] http://fmwww.bc.edu/repec/bocode/t/transint.html 
[2] https://www.r-statistics.com/2013/05/log-transformations-for-skewed-and-wide-distributions-from-practical-data-science-with-r/
[3] Alexandru Papiu's notebook https://www.kaggle.com/apapiu/house-prices-advanced-regression-techniques/regularized-linear-models/commentsnotebook 
We apply the function $x \rightarrow \log(1 + x)$ because it is always positive for $x \geq 0$
End of explanation
def read():
    Read training and test data and return a dataframe with ['Dataset','Id'] multi-index
    
    raw_train = pd.read_csv('../data/train_prepared_light.csv')
    raw_test = pd.read_csv('../data/test_prepared_light.csv')
    df = pd.concat([raw_train, raw_test], keys=['train', 'test'])
    df.index.names = 'Dataset', 'Id'
    return df
    
df = read()
df.shape
df.head()
df.tail()
Explanation: Merge the training and test datasets for data preparation
We're going to explore the training dataset and apply some transformations to it (fixing missing values, transforming columns etc). We'll need to apply the same transformations to the test dataset. To make that easy, let's put the training and test datasets into one dataframe.
End of explanation
pp = samlib.Pipeline(df.copy())  
assert pp == df  # the pipeline output equals df
Explanation: Initialize pipeline with raw data. We can always get the data and apply all the transformations in the pipeline by calling pp().
End of explanation
df.columns, len(df.columns)
Explanation: Select Numerical features
The dataset is wide with 78 features. Create dataframe containing the numerical features only.
End of explanation
df.dtypes.value_counts()
Explanation: We've got 3 data types: int, float and object
End of explanation
is_categorical = (df.dtypes == object)
is_numerical = ~is_categorical
dfnum = df.loc[:, is_numerical].copy()
dfnum.columns, len(dfnum.columns)
Explanation: Split the data between categorical and numerical features
End of explanation
dfnum.describe()
Explanation: We've got 36 numerical features. We can use the describe method to get some statistics:
End of explanation
def select_numerical_features(df):
    return df.loc[:, df.dtypes != object]
pp.append(select_numerical_features)
# Check the pipline
pp == dfnum
Explanation: But that's a lot of numbers to digest. Better get started plotting!
End of explanation
cols_with_nulls = dfnum.columns[dfnum.isnull().sum() > 0]
cols_with_nulls
dfnum[cols_with_nulls].isnull().sum().sort_values(ascending=False)
Explanation: Deal with NaN values
End of explanation
# We may want to refine this in the future. Perhaps build a model to predict the missing GarageCars from the other features?
median_list = 'LotFrontage', 'BsmtFullBath','BsmtHalfBath', 'GarageCars', 'GarageArea', 'MasVnrArea', 'BsmtFinSF1', 'BsmtFinSF2', 'TotalBsmtSF', 'BsmtUnfSF'
zero_list = []
def fillnans(dfnum):
    return dfnum.pipe(samlib.fillna, 'median', median_list)\
     .pipe(samlib.fillna, lambda df: 0, zero_list)\
     .assign(GarageYrBlt=dfnum.GarageYrBlt.fillna(
                dfnum.YearBuilt[dfnum.GarageYrBlt.isnull()]))  # fill with year garage was built
dfnum = fillnans(dfnum)
# Check that we got rid of the nulls
assert not samlib.has_nulls(dfnum)
pp.append(fillnans)
# Check the pipline
pp == dfnum
Explanation: Based on the description, the null values for the MasVnrArea should be 0 (no massonry veneer type)
End of explanation
def order_columns(df):
    return df.reindex_axis(df.columns.sort_values(), 1)
pp.append(order_columns)
pp().head()
dfnum = pp()
dfnum.head()
Explanation: Order columns in alphabetical order
End of explanation
dfnum.shape
samlib.featureplot(dfnum, ncols=6, nrows=6, figsize=(12, 4))
Explanation: Plot violinplots for each feature
The violin plots give us some idea of the distribution of data for each feature. We can look for things like skewness, non-normality, and the presence of outliers.
End of explanation
fig, ax = plt.subplots(1,1, figsize=(4, 4))
sns.distplot(dfnum.ScreenPorch, ax=ax)
ax.set_title('Distribution of ScreenPorch')
Explanation: Many of the features are higly skewed and some have very long tails. Some have discrete values (YrSold, Fireplaces).
The features with very long and thin tails, such as ScreenPorch, are almost constant (blobs with long tail) as can be seen below
End of explanation
def test_nearly_constant(series):
    counts = series.value_counts()
    max_val_count = max(counts)
    other_val_count = counts.drop(counts.argmax()).sum()
    return other_val_count / max_val_count < 0.25
is_nearly_constant = dfnum.apply(test_nearly_constant)
is_nearly_constant.value_counts()
dropme = dfnum.columns[is_nearly_constant]
dropme
def drop_constant_features(df):
    return df.drop(df.columns[df.apply(test_nearly_constant)], axis=1)
pp.append(drop_constant_features)
pp == dfnum.drop(dropme, axis=1) 
dfnum = dfnum.drop(dropme, axis=1)
Explanation: Drop nearly constant features
End of explanation
fig, axes = plt.subplots(1,2, figsize=(8, 4))
sns.distplot(dfnum['LotArea'], ax=axes[0])
sns.distplot(np.log1p(dfnum['LotArea']), ax=axes[1])
Explanation: Log transform the other features if they have a high skewness
Using a log transformation for some of the skewed features should help, as illustrated below. We use the raw data (not the standardized one) because we need positive values for the log function (we'll standardize the transformed variables later).
End of explanation
def skewtest(train, sort=True, ascending=True):
    Return dataframe of zfactor and pvalue for skew test
    test = sp.stats.skewtest(train)
    zfactor = test[0]
    pvalue = test[1]
    df = pd.DataFrame(dict(zfactor=zfactor, pvalue=pvalue), index=train.columns)
    if sort:
        return df.sort_values(by='zfactor', ascending=ascending)
    else:
        return df
skewtest(dfnum).head()
def is_skewed(train, min_zfactor=10, plot=False):
    Return series of booleans indicating whether a column is skewed or not.
    
    sk = skewtest(train)
    if plot:
        plt.figure(1)
        plt.title('Z-factor distribution from skewtest')
        plt.xlabel('Z-factor')
        sns.distplot(sk.zfactor)
        plt.figure(2)
        sk.zfactor.plot(kind='barh')
        plt.title('Z-factor for skewtest')
    return sk.zfactor > min_zfactor
is_skewed(dfnum, min_zfactor=10, plot=True)
Explanation: Use dataframe & series whenever possible for maximum flexibility (see below)
End of explanation
def transform_skewed_colums(dfnum):
    
    dfnum: dataframe to transform
    dropme: columns to drop
    is_skewed: iterable of length dfnum.columns indicating if a column is skewed
    
    dfnum2 = dfnum.copy()
    skewed_colz = is_skewed(dfnum)
    dfnum2.loc[:, skewed_colz] = dfnum2.loc[:, skewed_colz].apply(np.log1p)
    return dfnum2
pp.append(transform_skewed_colums)
# the transformed dataset has fewer columns and we only want those
dfnum2 = pp()
dfnum2.columns
is_skewed(dfnum2)
sorted(sp.stats.skewtest(dfnum2)[0])
zfactors2 = sp.stats.skewtest(dfnum2)[0]
pd.Series(data=zfactors2, index=dfnum2.columns)[is_skewed(dfnum)].sort_values().plot(kind='barh')
Explanation: Let's apply a log1p transform to all these and plot the distributions again
End of explanation
skewed = is_skewed(dfnum)
skewed.value_counts()
dfnum.shape
samlib.featureplot(dfnum2.loc[:, skewed], nrows=3, ncols=6, figsize=(10,3))
samlib.featureplot(dfnum2.loc[:, ~skewed], nrows=2, ncols=5, figsize=(10, 3))
Explanation: Now our originally skewed features look more symmetric. 
Check that the distributions are less skewed
End of explanation
dfnum2.to_csv('transformed_dataset_dfnum2.csv', index=True)
Explanation: Save transformed numerical data
Use the storage magic to communicate between notebooks.
End of explanation
def correlation(train, target_t):
    corr = pd.DataFrame(data=train.apply(lambda feature: sp.stats.pearsonr(feature, target_t['SalePrice'])), 
                        columns=['pearsonr'])
    corr = corr.assign(correlation=corr.applymap(lambda x: x[0]),
                       pvalue=corr.applymap(lambda x: x[1]))
    corr = corr.drop('pearsonr', axis=1)
    return corr.sort_values('pvalue', ascending=False)['correlation']
correlation(dfnum2.loc['train', :], logtarget).plot(kind='barh')
Explanation: Correlations
We're now in a good position to identify the key numerical features. Those should be hightly correlated with the sale price.
End of explanation
def sort_columns_by_correlation(dfnum2, target_t=logtarget):
    corr = correlation(dfnum2.loc['train',:], target_t)
    return dfnum2.reindex_axis(corr[::-1].index, axis=1)
#sort_columns_by_correlation(dfnum2)
pp.append(sort_columns_by_correlation)
pp().head()
pp().to_csv('transformed_numerical_dataset.csv', index=True)
Explanation: Sort columns in dfnum2 by correlation.
End of explanation
train = pp().loc['train'].assign(target=logtarget)
samlib.featureplot2(train, ncols=4, size=3, aspect=1.0, plotfunc=sns.regplot, y="target", data=train)
#(train.iloc[:, :-1], ncols=4, nrows=7, plotfunc=scatter, figsize=(12,3))
Explanation: Scatter plots
The correlation can be high even if there is no linear relationship between a feature and the target, so we should check the scatter plots.
End of explanation
cols_with_zeros = ['OpenPorchSF', 'MasVnrArea', 'TotalBsmtSF', 'WoodDeckSF', 'BsmtUnfSF', 'BsmtFinSF1', '2ndFlrSF']
not_oktrain = train.loc[:, cols_with_zeros + ["target"]]
samlib.featureplot2(not_oktrain, ncols=4, size=3, aspect=1.0, 
                    plotfunc=sns.regplot, y="target", data=not_oktrain)
Explanation: Some features have some sort of bi-modal distribution with a lots of 0 values.
End of explanation
notok = not_oktrain[not_oktrain != 0].drop('target', axis=1)
#correlation(notok, logtarget)
notok.head()
def correlation2(train, target=logtarget['SalePrice'], ignorena=True):
    
    def corr(series):
        if ignorena:
            mask = ~series.isnull()
            return sp.stats.pearsonr(series[mask], target[mask])
    
    df = pd.DataFrame(data=train.apply(corr), columns=['pearsonr'])
    return df.assign(pearson=df.applymap(lambda x: x[0]), pvalue=df.applymap(lambda x: x[1])).drop('pearsonr', axis=1)
notok_corrs = correlation2(notok).sort_values('pearson')
corrs_with_zeros = correlation(not_oktrain, logtarget).reindex(notok_corrs.index)
corrs_without_zeros = notok_corrs['pearson']
pd.DataFrame(dict(with_zeros=corrs_with_zeros, no_zeros=corrs_without_zeros)).plot(kind='barh')
plt.title('Effect of removing zeros on correlation')
not_oktrain2 = not_oktrain.reindex_axis(notok_corrs.index[::-1], axis=1);
def regplot_dropzeros(data=not_oktrain, drop_zeros=False, **kwargs):
    col = data.columns[0]
    if drop_zeros:
        mask = data[col] != 0
        xt = data[mask].squeeze()
        yt = logtarget[mask].squeeze()
    else:
        xt = data.squeeze()
        yt = logtarget.squeeze()
    sns.regplot(xt, yt, **kwargs)
Explanation: Dealing with the zeros
Let's recompute the correlations after dropping the zeros.
End of explanation
samlib.featureplot(not_oktrain2, ncols=4, nrows=2, figsize=(12, 3), plotfunc=regplot_dropzeros, drop_zeros=False)
Explanation: With the zeros
End of explanation
samlib.featureplot(not_oktrain2, ncols=4, nrows=2, figsize=(12, 3), plotfunc=regplot_dropzeros, drop_zeros=True)
Explanation: Without the zeros
End of explanation
import statsmodels.api as sm
from statsmodels.imputation import mice
df = pp()
# Replace zeros by NaNs
df[df.loc[:, cols_with_zeros] == 0] = np.nan
df = df.rename_axis({'1stFlrSF':'FrstFlrSF', '2ndFlrSF':'SndFlrSF'}, axis=1)
df.loc['train','SalePrice'] = logtarget.values
df.head()
samlib.has_nulls(df)
imp = mice.MICEData(df)
imp.update_all()
imp.data.head()
cols = ['OpenPorchSF', 'MasVnrArea', 'TotalBsmtSF', 'WoodDeckSF', 'BsmtUnfSF', 'BsmtFinSF1', 'SndFlrSF','SalePrice']
imputed_notok = imp.data.loc[:, cols]
imputed_notok.columns
samlib.featureplot2(imputed_notok, ncols=4, size=3, aspect=1.0, 
                    plotfunc=sns.regplot, y="SalePrice", data=imputed_notok)
imp.data.shape
df.shape
imp.data.index=df.index
imp.data.head()
imp.data.reindex_like(df).to_csv('transformed_numerical_dataset_imputed.csv', index=True)
Explanation: Data imputation
Impute the zero values from the rest of the dataset for the "not ok" columns.
End of explanation
<END_TASK>