Coursera: Data Analysis and Interpretation Specialization

Assignments for Coursera's Data Analysis and Interpretation Specialization courses

[mlfanatic]

Data Analysis and Interpretation Capstone - Milestone Assignment 3: Preliminary Results

Project title: Predicting the popularity of online content before it is published.

This weeks assignment was to post a draft of our preliminary results section from our final report. The research question can be found in my prior post here.

-----

Preliminary Results

LASSO Regression Analysis

Figure 2 below includes all 58 explanatory variables. This figure shows that as each explanatory variable is added to the model the Mean Square Error (MSE) declines until additional predictors don’t further reduce the MSE.

                   Fig 2. LASSO Mean squared error on each fold

Table 1 below provides definitions for the top explanatory variables that were chosen by LASSO above.

Table 1 - definitions for the top explanatory variables.

Determining the ‘average’ and ‘worst’ keywords was determined by Fernandes, et al. [2] using the following process: “we rank all article keyword average shares (known before publication), in order to get the worst, average and best keywords. For each of these keywords, we extract the minimum, average and maximum number of shares.”. Keywords shared an average number of times are considered ‘average’ keywords. Keywords shared the least number of times are considered the “worst” keywords.

Table 2 below shows the LASSO mean squared errors across 58 explanatory variables. The small difference between the training and test data sets below shows that the model is good since the test result MSE is close to the training result MSE.

Table 2. LASSO MSE and R-squared.

Fig 3 below shows that ‘self reference avg shares’ and ‘avg shares of avg keywords’ were most strongly associated with sharing a Mashable article. When coupled with additional top explanatory variables chosen by LASSO (LDA_topic_0 and avg_shares_of_worst_keywords) implies that articles about a certain topic where particular keywords are present tend to be shared more often than other articles.

Fig 3. LASSO regression coefficients

Table 3 below shows the top chosen explanatory variables. These variables will be the focus for the remainder of this study.

Across all variables below the Pearson correlation coefficient ( r) and coefficient of determination (r^2) are quite low.

Table 3 - The top chosen explanatory variables.

Descriptive Statistics

Table 4 shows descriptive statistics for the top explanatory variables and the online shares response variable. These variables have already been converted via the log transform as mentioned above. Prior to the log transform there were 3,374 average shares per article with a minimum of 1 share for 1 article and a maximum of 690,400 shares for another article. There were 39,644 unique articles in this study.

Below are histograms for the top explanatory variables and the response variable. All variables below show a steep positive or negative skew except for number of hrefs and article subjectivity which appear more normally distributed.

<Note to reviewers: My final report will have more histograms. This is just a sample. >

Fig 4. Distribution of top explanatory variables and the online shares response variable. 

Bivariate Analyses

Overall, the scatter plots below (Figure 5) don’t show a strong correlation among the top explanatory variables and the shares response variable. The articles with ‘avg shares of avg keywords’ had the strongest correlation to online shares (Pearson r=0.18, p< 0.0001). The ‘number of hrefs’ is likely a spurious correlation since the number of hrefs within the primary article shouldn’t sway a reader to share that article more or less.

<Note to reviewers: My final report will have more scatter plots. This is just a sample. >

Fig 5. Association between the top explanatory variables and online shares response variable.

Regression Analysis

It is generally accepted that multiple regression models are used to try to find one that can best predict the response variable. Across all models below the overall p-value was significant. 24 regression models were run each varying the top explanatory variables. The best overall results are included below. The best performing model is marked below in bold.

Table 5. Regression analysis results.

The model above in bold is the best because:

All coefficient p-values are significant (p-values = 0.000 < 0.05 alpha)

The overall model p-value is significant (p-value = 0.000 < 0.05 alpha)

This model has the highest overall adjusted r-squared (r^2 = 0.072) which is a measure of variance.

We can conclude that the chosen regression model can predict online shares fairly well.

The quantile-to-quantile plot in figure 6 below was created for the chosen model above. If the residuals (blue dots) were normally distributed they would follow the red line. However, the residuals generally follow the red line at the middle quantiles but slightly deviate at the lower and higher quantile. This deviation may have several causes:

There may be a curvilinear association between the variables included in this model.

The response variable has several large values well beyond the mean and 3 standard deviations from the mean.

There might be other explanatory variables that could be included in our model, beyond those in the original data set, that could improve estimation.

Fig 6. Quantile to quantile plot for the chosen model: Multiple regression with top explanatory variables.

Conclusions/Limitations

Overview of key findings

To predict the popularity of online articles prior to publication this project used LASSO regression and other regression models. The data set included N = 39,644 samples that describe articles published on Mashable.com from January 7, 2013 to January 7, 2015.

The first step used relaxed LASSO and Pearson correlation to find the best subset of explanatory variables that could predict whether online articles would be ‘shared’ by readers when published. The strongest predictors of online shares ‘self reference avg shares’ and ‘avg shares of avg keywords’ and whether the article was published on the weekend. On one hand this implies that articles about a certain topic have a higher likelihood of being shared.

The best performing model had a low LASSO MSE between the training and test data sets as demonstrated in table 2 which means the model fits the training data well. The chosen model has an overall p-value lower than 0.05 alpha and an adjusted r-squared of .072 which indicates that the model has some good predictive ability. However, the quantile-to-quantile plot in fig 6 shows that this model had a fair number of residuals in the lower and upper quantiles which suggests a curvilinear relationship between the top variables or possibly other issues.

In summary, when comparing the overall results above to prior research:

Article subjectivity was one of the top ten predictors of online shares as seen in table 3. The subjectivity histogram in figure 4 above also shows articles that are slightly less subjective (more objective) is one indicator to determine whether it will be shared. However, according to Bandari, et al. [1] article subjectivity was not a good predictor towards sharing news blogs.

If a Mashable article was published on a weekday was somewhat important to predict subsequent shares. Table 3 shows the variable ‘is weekend’ was ranked #3 by LASSO regression coefficients. Also, a tentative relationship can be seen in the scatter plot in figure 4. If a Mashable article was published on a weekday it had a slightly better chance of being shared. This result seems to match the results from Szabo, et. al [4] in which sharing articles, via Digg, occurred more frequently on weekdays than on weekends.

 Implications

Predicting human behavior is tricky such as predicting whether a user will share an online article. At the outset the ‘right’ variables may be present in the data, however they may not be sufficient to predict human behavior.

In addition, the results may show that a particular audience may prefer articles about a certain topic and thus those articles are more likely to be shared.

 Limitations

The chosen data set was from a single website, Mashable.com, which may draw a particular audience. The results described above may not generalize well to other websites or audiences.

Pearson’s correlation across all 58 explanatory variables and subsequent scatter plots of the top explanatory variables showed that there was no strong correlation between the explanatory variables and the response variable. This of course affects the ability to predict an outcome.

 Future directions

Based on prior research [1], [2], [3], [4], it does seem possible to accurately predict the popularity of online articles before publication. More advanced algorithms may be needed to make accurate predictions, algorithms such as Adaptive Boosting, Support Vector Machines, and Naive Bayes.

It would also be interesting to determine how much influential people acting as intermediaries increase an articles “share-ability”.

Reference

[1] Roja Bandari (University of California, Los Angeles), Sitaram Asur (HP Labs), Bernardo A. Huberman (HP Labs). The Pulse of News in Social Media: Forecasting Popularity 2012, International Conference on Weblogs and Social Media. http://arxiv.org/pdf/1202.0332.pdf

[2] K. Fernandes, P. Vinagre and P. Cortez. A Proactive Intelligent Decision Support System for Predicting the Popularity of Online News. Proceedings of the 17th EPIA 2015 - Portuguese Conference on Artificial Intelligence, September, Coimbra, Portugal. http://repositorium.sdum.uminho.pt/bitstream/1822/39169/1/main.pdf

[3] Sasa Petrovic, Miles Osborne,  Victor Lavrenko. RT to Win! Predicting Message Propagation in Twitter. 2011, National Conference on Artificial Intelligence. http://www.aaai.org/ocs/index.php/ICWSM/ICWSM11/paper/view/2754/3209

[4] Gabor Szabo (Hewlett-Packard), Bernardo A. Huberman (Hewlett-Packard). Predicting the popularity of online content. 2010, Communications of The ACM, volume 53, issue 8, pp 80-88. https://www.researchgate.net/profile/Gabor_Szabo10/publication/23417017_Predicting_the_Popularity_of_Online_Content/links/00463529e2169e339e000000.pdf?disableCoverPage=true

[5] Chen, Edwin. Introduction to Latent Dirichlet Allocation. The link below was last retrieved on August 18, 2016. http://blog.echen.me/2011/08/22/introduction-to-latent-dirichlet-allocation/

[6] “Advantages of doing “double lasso” or performing lasso twice?” (aka Relaxed LASSO). Last accessed the following URL on August 20, 2016. http://stats.stackexchange.com/questions/37989/advantages-of-doing-double-lasso-or-performing-lasso-twice  

[7] Meinshausen, Nicolai. Relaxed Lasso. http://stat.ethz.ch/~nicolai/relaxo.pdf

[mlfanatic]

Data Analysis and Interpretation Capstone - Milestone Assignment 2: Methods

Project title: Predicting the popularity of online content before it is published.

This weeks assignment was to post a draft of our methods section from our final report. The research question can be found in my prior post here.

Note to reviewers: the final report will have a references section that will resolve the references below such as [3].

-----

Method

Sample

The data set was first created for research by K. Fernandes, P. Vinagre and P. Cortez. A Proactive Intelligent Decision Support System for Predicting the Popularity of Online News. Proceedings of the 17th EPIA 2015 - Portuguese Conference on Artificial Intelligence, September, Coimbra, Portugal.

The data set included N = 39,645 samples that describe articles published on Mashable (www.Mashable.com) from January 7, 2013 to January 7, 2015. The actual text of the articles remains the property of Mashable, however the sample contains variables that describe each article.

The UCI Machine Learning repository stores this data set.

Measures

The response variable, number of shares per article, was recorded when each Mashable user shared an article.

There are 58 explanatory variables that represent the attributes that influence Mashable online content “shares” including”:

- attributes about the article such as the number of the words in the title and content body

- keywords describing the article, number of keywords, etc.

- whether the article belongs to one of six data channels (e.g.: lifestyle, entertainment, etc.)

- number of links included within the article

- number of images embedded in the article

- article subjectivity and sentiment were pre-calculated with the Pattern web mining module (http://www.clips.ua.ac.be/pattern)

Analyses

The analysis consisted of several steps to reduce the number of explanatory variables from 58 to the top variables with the most strength to predict online shares. Once the top explanatory variables were determined then more in-depth algorithms were used.

First, several ‘summary’ explanatory variables were removed since they represented another variable already in the data set. For example, the ‘average article sentiment polarity’ was removed since the article sentiment polarity was represented. There were 41 explanatory variables left after the summary explanatory variables were removed.

Next, “relaxed LASSO” [6], [7] and Pearson correlation were used to reduce the 41 explanatory variables to the top variables that best predicted online shares. Data was centered to have a mean = 0 and a standard deviation = 1 then LASSO and Pearson correlation were used. The relaxed LASSO and Pearson correlation used a random sample of 60% training data (N=23,786) and 40% test data (N=15,857). LASSO used K-fold cross validation was used with 10 folds. Prediction accuracy was determined by using the mean square error rate when the training data prediction was applied to observations derived from the test data set.

Once the top explanatory variables were determined, to better understand the shape of the data a bivariate analysis was performed which included the online shares response variable. The mean, standard deviation, minimum and maximum values were calculated. Also, histograms and scatter plots were visually examined.

Because the top explanatory variables and the online shares response variable are quantitative, several regression models were used to determine which one best predicted the online shares response variable.

Logistic regression requires a binary response variable. Since online shares contained continuous values the median was used to separate online shares and create a balanced class distribution as follows; less than 1,400 shares was considered “no share” whereas more than 1,400 shares was considered a “share”. This same approach was also used by Fernandes, et al. [2] on the same data set.

To determine how well each regression model predicted the response variable, the models were evaluated using adjusted R-squared, p-values, confidence intervals and residual plots.

The overall process described above is summarized in figure 1 below.

The Coursera capstone project maximum word limit would have been exceeded if more detailed results were included for selecting the top explanatory variables. The results below will only focus on the top explanatory variables chosen by the relaxed LASSO step as described above. The top variable definitions in table 1 below provide more context for the remainder of this paper.  

[mlfanatic]

Data Analysis and Interpretation Capstone - Week 1 assignment - Identify Your Data and Research Question

Project title: Predicting the popularity of online content before it is published.

The purpose of this study is to identify the best predictors of online content popularity before that content is published. Popularity is determined by the number of online content "shares". The factors that influence online content "shares" include the content, keywords, day of the week the article was published, topic and sentiment.

For decades I have found the ability of a computer to "understand" the written word to be compelling. Further, it is a compelling challenge to predict how the written word will be received by the reader and, if motivated by what they've read, passed onto another person.

The ability to predict the popularity of online content has several implications. It allows content creators the ability to adjust their content before it is published to be more interesting and thus more widely read. It also offers content creators the ability to contribute more effectively to their organizational business goals such as online ad revenue targets.

Advertisers may also be willing to pay more to be included in an online article if it is predicted to receive a certain number of views.

Lastly, news recommendation services and search engines may also use the predicted content popularity as another variable to rank newer content after it has been published.

[mlfanatic]

Machine Learning for Data Analysis - Week 4 assignment - Running a K-means Cluster Analysis

Executive summary

The assignment for this week is to run a K-means cluster analysis. I used the same data set that I described in my prior post here.

In summary, K-means cluster analysis is an unsupervised, exploratory technique that attempts to find similar sub-groups of data within a larger dataset. 'Unsupervised' simply means that a pre-determined outcome has not been provided in advance for the algorithm to 'learn' what the outcome should look like. Practical applications of clustering include customer segmentation, demographic studies, grouping similar documents, (e.g.: customer review feedback), anomaly detection (e.g.: potentially fraudulent credit-card transactions), grouping similar images, etc.

For this particular assignment, a K-means cluster analysis was conducted to provide more insight into demographic and payment history data that may have an impact on whether customers will default on their next credit card payment. The clustering results below are largely inconclusive. The standardized data is not easily separable into clusters. Also, some of the cluster centers (aka: centroids) overlap. In addition, an Analysis of Variance (ANOVA) was conducted to determine if there were significant differences between the clusters and credit card default rate. Significant differences were found, however the R-squared was quite low.

Also, per the assignment instructions the cluster analysis was run once on a training data set.

Details

The following explanatory variables were included in the cluster analysis: credit limit, gender, education, marital status, age and a variable called 'X6_PAY_0'. X6_PAY_0 tracks repayment status for the prior 9+ months in the data set. This variable was included versus other payment variables because my prior LASSO regression results and my random forest results showed that this variable had the greatest ability to predict credit card payment default.

The data below provides some background behind the variables we are using. This summary was created before the data was standardized or split into training and test sets as described below.

Fig 1.

These variables were standardized to have a mean of zero and a standard deviation of 1 prior to clusters being formed.

The data were randomly split into a training set that included 70% of the observations (N=21,000) and a test set that included 30% of the observations (N=9,000). Only the training set was used.

Since K-means clustering requires the number of clusters to be specified in advance, a series of cluster analysis were performed on the training set to determine the optimal number of clusters.

To improve the results, the K-means algorithm was initialized to perform 10 iterations where the starting centroid location was randomly chosen at the start of each iteration. Within each iteration the centroid was recalculated 100 times at most to best represent the center of a cluster. At the end of each iteration the best centroid was chosen programmatically.

Next, to visualize a few clusters a maximum number of 3 clusters were chosen. 3 clusters were chosen since more clusters would be difficult to visually represent.

Lastly, an analysis of variance (ANOVA) was run to to test for significance between clusters and to predict credit card default payment next month.

Results

Fig 2.

The results above show that 6 or 7 seems to be the optimal number of clusters. Typically a clear "elbow" appears in the line that marks a point of diminishing returns. However, in this case the slope of the line becomes flatter when it reaches 6 or 7 clusters.

The remaining results below use 3 clusters.

Fig 3.

The means of the clustering variables shows that customers in cluster 0:

had the lowest credit limit balance (-0.44)

tend not to be married (0.90)

were the youngest (-0.66)

The means of the clustering variables shows that customers in cluster 1:

had the highest credit limit balance (1.26)

the highest education (-0.55)

the lowest mean for the X6_PAY_0 variable (-0.66). This variable tracks the number of months that the customer is behind payment. Each prior month where there was a default is tracked as a sequentially increasing positive number. Customers in this cluster tended to default the least.

The means of the clustering variables shows that customers in cluster 2:

had the lowest education (0.62)

tend to be married (-0.86)

tend to be the oldest (0.61)

tend to default the most (0.30)

Fig 4.

The scatter plot above is a little disappointing. All variables were standardized thus their values ranged from -2.97 to 7.13 so they tended to overlap. The clustering results are inconclusive since 2 of the 3 cluster centroids overlap - the middle star actually marks 2 centroids. 

 ANOVA results

Fig 5.

Fig 6.

The 'Prob (F-statistic)' (p-value<0.001) is much lower than alpha (0.05) which leads us to reject the null hypothesis.

A Tukey post hoc comparison between clusters was performed. The results also indicate significant differences between clusters towards predicting credit card default (p<0.001).

However, the R-squared is quite low (0.018), which means that the 6 explanatory variables mentioned above only explain 1.8% of the variation when predicting the credit card default response variable.

Fig 7.

When analyzing the mean and standard deviation of payment default by cluster, cluster 0 had the lowest credit card defaults (mean=0.13, SD=0.34). Cluster 1 had the highest credit card defaults (mean=0.27, SD=0.45) with customers in cluster 2 having a very close defaults (mean=0.24, SD=0.43). These results are not consistent with the means of the clustering variables analysis in Fig 3 above.

How can 2 different analysis of credit card default rates show 2 different results across the same clusters? Well, the results in Fig 3 above uses prior historical data represented in the explanatory variable X6_PAY_0. On the other hand, the results in Fig 7 use the 6 explanatory variables to predict the forward looking response variable Y_default_payment_next_month.

Appendix

Python code

=====

from pandas import DataFrame import pandas as pd import numpy as np import matplotlib.pylab as plt from sklearn.cross_validation import train_test_split from sklearn.preprocessing import StandardScaler from sklearn.cluster import KMeans from sklearn.decomposition import PCA from scipy.spatial.distance import cdist import statsmodels.formula.api as smf import statsmodels.stats.multicomp as multi

# Set numerical display formats. pd.options.display.float_format = '{:0,.2f}'.format

# Remove the preprocessing.scale warning since I know clusterDF is a copy. pd.options.mode.chained_assignment = None

# --------------- # Load the dataset, Drop NA's and show overall dataframe summary. defaultRateloadDF = pd.read_csv("default_rate_credit_card_clients.csv")

defaultRateDF = defaultRateloadDF.dropna()

# --------------- # Get a subset of variables I need for clustering. clusterDF = defaultRateDF[['X1_LIMIT_BAL', 'X2_SEX', 'X3_EDUCATION', 'X4_MARRIAGE', 'X5_AGE', 'X6_PAY_0']] print('Variable overview prior to standardizing:\n', clusterDF.describe())

# Standardize clustering variables to have mean=0 and sd=1. # Rename the column names to be more readable. # All columns are already float 64.

stdsc = StandardScaler() clusterDF = pd.DataFrame(stdsc.fit_transform(clusterDF)) clusterDF = clusterDF.rename(columns={0: 'X1_LIMIT_BAL',                                      1: 'X2_SEX',                                      2: 'X3_EDUCATION',                                      3: 'X4_MARRIAGE',                                      4: 'X5_AGE',                                      5: 'X6_PAY_0'}) print("clusterDF types:\n", clusterDF.dtypes)

# --------------- # Split data into train and test sets and run K-means clustering. clus_train, clus_test = train_test_split(clusterDF, test_size=.3, random_state=123)

# k-means cluster analysis                                                           MIN_CLUSTERS = 1 MAX_CLUSTERS = 11 clusters = range(MIN_CLUSTERS, MAX_CLUSTERS) meandist = []

print("\nIteration number and inertia result") print("-----------------------------")

for k in clusters:

   # n_init will run KMeans the specified number times and randomly choose    # the starting centroid location each time.    # max_iter specifies the maximum number of times KMeans will recalc the    # centroid in each iteration.    model = KMeans(n_clusters=k, n_init=10, max_iter=100)    #model = KMeans(n_clusters=k)    model.fit(clus_train)    print(k, model.inertia_)

   # model.predict obtains the results of the cluster analyses. It has      # the closest cluster that each observation belongs to.    clusassign = model.predict(clus_train)

   meandist.append(sum(np.min(cdist(clus_train, model.cluster_centers_, 'euclidean'), axis=1))    / clus_train.shape[0])

""" Plot average distance from observations from the cluster centroid to use the Elbow Method to identify number of clusters to choose """

plt.plot(clusters, meandist) plt.xlabel('Number of clusters') plt.ylabel('Average distance') plt.title('Selecting k with the Elbow Method') plt.show()

# Interpret 3 cluster solution. # This is necessary since a 3-D space can be visualized but anything larger # would be difficult/impossible to visualize. #model3 = KMeans(n_clusters=3) model3 = KMeans(n_clusters=3, n_init=10, max_iter=100) model3.fit(clus_train) clusassign = model3.predict(clus_train)

# plot clusters pca_2 = PCA(2) plot_columns = pca_2.fit_transform(clus_train)

plt.scatter(x=plot_columns[:,0], y=plot_columns[:,1], c=model3.labels_) plt.scatter(model3.cluster_centers_[:,0], model3.cluster_centers_[:,1], c='red', marker='*', s=200, label='Centroids') plt.xlabel('Canonical variable 1') plt.ylabel('Canonical variable 2') plt.title('Scatterplot of Canonical Variables for 3 Clusters') plt.legend() plt.grid() plt.show()

""" BEGIN multiple steps to merge cluster assignment with clustering variables to examine cluster variable means by cluster

Note, most of the code below is unchanged since that’s just the way you do it.

""" # create a unique identifier variable from the index for the # cluster training data to merge with the cluster assignment variable clus_train.reset_index(level=0, inplace=True)

# create a list that has the new index variable cluslist=list(clus_train['index'])

# create a list of cluster assignments labels=list(model3.labels_)

# combine index variable list with cluster assignment list into a dictionary newlist=dict(zip(cluslist, labels)) #print("Newlist dictionary results:\n", newlist)

# convert newlist dictionary to a dataframe newclus=DataFrame.from_dict(newlist, orient='index') #print("\nNewclus results:\n", newclus.head(10))

# rename the cluster assignment column newclus.columns = ['Cluster']

# Next, do the same for the cluster assignment variable # create a unique identifier variable from the index for the # cluster assignment dataframe # to merge with cluster training data newclus.reset_index(level = 0, inplace = True)

# merge the cluster assignment dataframe with the cluster training variable # dataframe by the index variable. merged_train=pd.merge(clus_train, newclus, on = 'index') merged_train.head(n=100)

# cluster frequencies merged_train.cluster.value_counts()

""" END multiple steps to merge cluster assignment with clustering variables to examine cluster variable means by cluster """

# FINALLY calculate clustering variable means by cluster clustergrp = merged_train.groupby('cluster').mean() print ("Clustering variable means by cluster") print(clustergrp)

# validate clusters in training data by examining cluster differences in credit card

# payment default using ANOVA. First merge credit card payment default with clustering # variables and cluster assignment data. CCdefaultDF = defaultRateDF['Y_default_payment_next_month']

# split credit card payment default data into train and test sets CCdefault_train, CCdefault_test = train_test_split(CCdefaultDF, test_size=.3, random_state=123) CCdefault_train1 = pd.DataFrame(CCdefault_train) CCdefault_train1.reset_index(level=0, inplace=True) merged_train_all=pd.merge(CCdefault_train1, merged_train, on='index') sub1 = merged_train_all[['Y_default_payment_next_month', 'cluster']].dropna()

CCdefaultmod = smf.ols(formula='Y_default_payment_next_month ~ C(cluster)', data=sub1).fit() print ("\n", CCdefaultmod.summary())

print ('\nMean credit card payment default by cluster') m1= sub1.groupby('cluster').mean() print (m1)

print ('\nStandard deviations for credit card payment default by cluster') m2= sub1.groupby('cluster').std() print (m2)

mc1 = multi.MultiComparison(sub1['Y_default_payment_next_month'], sub1['cluster']) res1 = mc1.tukeyhsd() print("\n", res1.summary())

# Debug: write the results to output file. print("\nWrite results to the output file.") sub1.to_csv('sub1.csv', sep=',', encoding='utf-8')

 =====

 Additional resources (not necessary for this assignment)

This blog presents how to use K-means inertia to find the optimal number of clusters. Below you’ll find my results using this method with the same data above.

In addition to the method provided above, I also used the K-means 'inertia_' variable for each of the 10 iterations described above to find the optimal number of clusters to use. Both methods agree on 6-7 clusters. According to the SKLean docs inertia_  is the "sum of distances of samples to their closest cluster center". The arrow shows the point of diminishing returns.  

If you don't have access to this Coursera course, this PDF provides an excellent overview of K-means clustering.

[mlfanatic]

Machine Learning for Data Analysis - Week 3 assignment - Running a LASSO Regression Analysis

Executive summary

This week’s assignment is to create a LASSO regression analysis.

LASSO means Least Absolute Selection Shrinkage Operator. As the name implies LASSO selects the best predictors that minimize prediction error. Also, it places a constraint on predictors such that less meaningful predictors towards zero, thus decreasing their importance in the overall linear regression equation.

Details

I used the same data that I described in my prior post. In summary, the data has 23 explanatory\predictor variables that attempt to predict credit card default. The explanatory variables include categorical and quantitative variables.

The original data set was split into two different sets; test set that contained 30% of the observations (N=9,000) and the remaining 70% into a training set (N=21,000).

All explanatory variables were standardized to have a mean = 0 and standard deviation = 1.

The LASSO least angle regression with k=10 fold cross validation was used with the training set first then validated against the test set.

Results

Out of the original 23 explanatory variables 18 were retained (in other words their regression coefficient was non-zero). These 18 variables accounted for 26.5% of the variance in the credit card payment default response variable.

The explanatory variables that were removed were: X13_BILL_AMT2, X14_BILL_AMT3, X15_BILL_AMT4, X16_BILL_AMT5, X17_BILL_AMT6 . These variables represent the credit card debt that a customer was billed in 5 of the last 6 months. This implies that with the exception of the most recently billed amount (see X12_BILL_AMT1 below) all other credit card billed amount variables do not predict whether a customer will default on their next payment.

Of the 18 explanatory variables that were retained the variable X6_PAY_0 had the largest positive predictive ability on credit card payment default (0.107). This variable tracks repayment status for the prior 9+ months. Essentially, did the customer pay on-time last month? If their payments have been delayed how many months were their payments delayed? This result makes sense since this variable tracks all prior payment states in our sample.

The explanatory variable X12_BILL_AMT1 had the largest negative predictive ability (-0.034). This variable tracks the dollar amount billed for the most recent month in our data.

In the chart above the topmost line is X6_PAY_0 . The purple line on the bottom is X12_BILL_AMT1. Both variables are described above.

The chart above shows that there is variability in each individual cross-validation fold in the training data set.

As explanatory variables are added the change in the mean square error follows a similar pattern for each fold. It eventually levels off to a point where adding more explanatory variables doesn't further reduce the mean square error.

The Mean Squared Error (MSE) is necessary to evaluate the performance of a predictor. The training and test set have a MSE that are very close (.152 and .149).

The R-square (aka, coefficient of determination) tells us how well changes in explanatory variable(s) can predict the response variable. Generally, it says that the explanatory variable(s) shares some portion of the variability with the response variable.

The training and test set R-Square are very close (.125 and .120). Essentially, the 18 explanatory variables retained in the training and test sets explain only about 12% of the variability in the response variable.

Alternatively, R-square can be used to show how much variability is unaccounted for simply by subtracting 1 from the R-square and multiplying by 100 to receive a percent. For example, using the results above: 1 - .12 = 0.88 * 100 = 88% of the variability between the explanatory variables and the response variable is unaccounted for.

In summary, our R-square results show us that changes in the explanatory variables cannot be used to accurately predict changes in the response variable with a high degree of confidence.

Appendix - Python code

===== import pandas as pd import numpy as np import matplotlib.pylab as plt from sklearn.cross_validation import train_test_split from sklearn.linear_model import LassoLarsCV from sklearn import preprocessing from sklearn.metrics import mean_squared_error

# Set numerical display formats. pd.options.display.float_format = '{:0,.3f}'.format

# Remove the preprocessing.scale warning. pd.options.mode.chained_assignment = None  

# --------------- # Load the dataset, Drop NA's and show overall dataframe summary.

defaultRateDF = pd.read_csv("default_rate_credit_card_clients.csv")

defaultRateDF = defaultRateDF.dropna()

# The data is already INT64 so no need to convert to_numeric #print('\nDefaultRateDF data types\n', defaultRateDF.dtypes) #print ("\nOriginal Dataframe describe() results\n", defaultRateDF.describe())

# ------------------------- # Prepare the data # Predictors are explanatory vars. # Target is the response variable.

predictors = defaultRateDF[['X1_LIMIT_BAL', 'X2_SEX', 'X3_EDUCATION', 'X4_MARRIAGE', 'X5_AGE', 'X6_PAY_0', 'X7_PAY_2', 'X8_PAY_3', 'X9_PAY_4',                     'X10_PAY_5', 'X11_PAY_6', 'X12_BILL_AMT1', 'X13_BILL_AMT2',                            'X14_BILL_AMT3', 'X15_BILL_AMT4', 'X16_BILL_AMT5',                            'X17_BILL_AMT6', 'X18_PAY_AMT1', 'X19_PAY_AMT2',                            'X20_PAY_AMT3', 'X21_PAY_AMT4', 'X22_PAY_AMT5',                            'X23_PAY_AMT6']]

target = defaultRateDF.Y_default_payment_next_month

# Standardize predictors to have mean = 0 and Std Dev = 1.

predictors['X1_LIMIT_BAL'] = preprocessing.scale(predictors['X1_LIMIT_BAL'].astype('float64')) predictors['X2_SEX'] = preprocessing.scale(predictors['X2_SEX'].astype('float64')) predictors['X3_EDUCATION'] = preprocessing.scale(predictors['X3_EDUCATION'].astype('float64')) predictors['X4_MARRIAGE'] = preprocessing.scale(predictors['X4_MARRIAGE'].astype('float64')) predictors['X5_AGE'] = preprocessing.scale(predictors['X5_AGE'].astype('float64')) predictors['X6_PAY_0'] = preprocessing.scale(predictors['X6_PAY_0'].astype('float64')) predictors['X7_PAY_2'] = preprocessing.scale(predictors['X7_PAY_2'].astype('float64')) predictors['X8_PAY_3'] = preprocessing.scale(predictors['X8_PAY_3'].astype('float64')) predictors['X9_PAY_4'] = preprocessing.scale(predictors['X9_PAY_4'].astype('float64')) predictors['X10_PAY_5'] = preprocessing.scale(predictors['X10_PAY_5'].astype('float64')) predictors['X11_PAY_6'] = preprocessing.scale(predictors['X11_PAY_6'].astype('float64')) predictors['X12_BILL_AMT1'] = preprocessing.scale(predictors['X12_BILL_AMT1'].astype('float64')) predictors['X13_BILL_AMT2'] = preprocessing.scale(predictors['X13_BILL_AMT2'].astype('float64')) predictors['X14_BILL_AMT3'] = preprocessing.scale(predictors['X14_BILL_AMT3'].astype('float64')) predictors['X15_BILL_AMT4'] = preprocessing.scale(predictors['X15_BILL_AMT4'].astype('float64')) predictors['X16_BILL_AMT5'] = preprocessing.scale(predictors['X16_BILL_AMT5'].astype('float64')) predictors['X17_BILL_AMT6'] = preprocessing.scale(predictors['X17_BILL_AMT6'].astype('float64')) predictors['X18_PAY_AMT1'] = preprocessing.scale(predictors['X18_PAY_AMT1'].astype('float64')) predictors['X19_PAY_AMT2'] = preprocessing.scale(predictors['X19_PAY_AMT2'].astype('float64')) predictors['X20_PAY_AMT3'] = preprocessing.scale(predictors['X20_PAY_AMT3'].astype('float64')) predictors['X21_PAY_AMT4'] = preprocessing.scale(predictors['X21_PAY_AMT4'].astype('float64')) predictors['X22_PAY_AMT5'] = preprocessing.scale(predictors['X22_PAY_AMT5'].astype('float64')) predictors['X23_PAY_AMT6'] = preprocessing.scale(predictors['X23_PAY_AMT6'].astype('float64'))

print ("\nPredictors. Verify mean = 0 and Std Dev = 1.\n", predictors.describe())

# --------------- # Split into training and testing sets.

pred_train, pred_test, tar_train, tar_test = train_test_split(predictors, target,                                                              test_size=.3, random_state=123)

# Specify and run the lasso regression model. # cv = 10 is cross validation with 10 random folds. # precompute = False if the model is not very large. model = LassoLarsCV(cv=10, precompute=False).fit(pred_train,tar_train)

# Print variable names and regression coefficients coef_results = list(zip(predictors.columns, model.coef_))

def getKey(item):    return item[1]

print('\nSorted coefficient results (importance low to high)') for item in sorted(coef_results, key=getKey):    print('%s: %0.3f' % (item[0], item[1])) print("\n")

# --------------- # Plot coefficient progression # Used negative log to make the results easier to read.

m_log_alphas = -np.log10(model.alphas_) ax = plt.gca() plt.plot(m_log_alphas, model.coef_path_.T) plt.axvline(-np.log10(model.alpha_), linestyle='--', color='k',            label='alpha CV')         plt.legend(defaultRateDF.columns[1:-1], loc='upper center', bbox_to_anchor=(0.5, -0.1)) plt.ylabel('Regression Coefficients') plt.xlabel('-log(alpha)') plt.title('Regression Coefficients Progression for Lasso Paths') plt.show()

# Plot mean square error for each fold m_log_alphascv = -np.log10(model.cv_alphas_) plt.figure() plt.plot(m_log_alphascv, model.cv_mse_path_, ':') plt.plot(m_log_alphascv, model.cv_mse_path_.mean(axis=-1), 'k',         label='Average across the folds', linewidth=2)

# Draw a dashed vertical line at the chosen model. plt.axvline(-np.log10(model.alpha_), linestyle='--', color='k',            label='alpha CV') plt.legend() plt.xlabel('-log(alpha)') plt.ylabel('Mean squared error') plt.title('Mean squared error on each fold') plt.show()

# MSE from training and test data train_error = mean_squared_error(tar_train, model.predict(pred_train)) test_error = mean_squared_error(tar_test, model.predict(pred_test))

print('') print ('Training data MSE: %0.3f' % train_error) print ('Test data MSE: %0.3f' % test_error)

# R-square from training and test data rsquared_train=model.score(pred_train,tar_train) rsquared_test=model.score(pred_test,tar_test)

print ('Training data R-square: %0.3f' % rsquared_train) print ('Test data R-square: %0.3f' % rsquared_test)

=====

[mlfanatic]

"Data Munge-O-Matic" - Machine Learning for Data Analysis - Week 2 assignment - Running a Random Forest

Executive summary

I've titled this blog post "Data Munge-O-Matic" since our assignment this week was to run a random forest. Random forests can use categorical and descriptive data as inputs/explanatory variables. Based on my previously chosen data set I used all 23 explanatory variables and shoved them through a random forest to find the best predictor of credit card default (a binary response variable). The results below show that some explanatory variables have higher predictive scores (aka importances) than others.

In summary, random forests are like an ensemble of individual decision trees working together to find the most important explanatory variables to predict an outcome. Random forests can’t be easily interpreted like decision trees, however they offer several advantages. For example, the only hyperparameter that *really* needs to be considered is the number of trees for the random forest. [Source, Python Machine Learning by Sebastian Raschka, pg 90.]

The results below were created by running a random forest 10 separate times with a maximum of 25 trees. Each run randomly pulled different data from the overall sample. The results across 10 runs are:

The variable 'X6_PAY_0' remained the most significant predictor with a score around 0.09. As our ‘best’ predictor this is quite low. This variable is the repayment status for the most recent month in our sample. Overall we can say that this variable is our best predictor to determine whether the credit card holder will default next month.

Model prediction accuracy remained around 0.80 or 80%, however details are provided below that describe why this result isn’t particularly impressive.

Only 7-8 trees were needed to provide the prediction accuracy above whereas 25 trees were used in each run.

A validation approach was used to split the original data set into 2 parts; 60% training and 40% test. The  random forest model was developed against the training set then the test set was used to predict the binary values in the response variable. There are several drawbacks with the validation approach that won't be discussed in this post.

 Results

Accuracy

In my previous post I mentioned that accuracy is the ratio of instances the classifier correctly classified. The random forest model was run 10 times and consistently achieved ~80% accuracy. However, accuracy best represents how well a classifier performs when using balanced classes/features/data columns. Overall, the data used by our random forest was unbalanced.

~80% accuracy sounds impressive, however when looking at the overall sample with 30,000 individual customer records 78% of the customers did not default, thus only 22% actually did default. Our model performed slightly better than randomly selecting a result from the sample data set.

To balance our response variable, credit card defaults, I could have deleted 16,728 non-default records so that the data comprised of 6,636 customers who did default and 6,636 who did not default, however this opens up a whole new set of questions outside the scope of this assignment (e.g.: Which records to delete? Will this create an imbalance in other features/columns like customer age, education level, gender, etc.)

Also, in the future we will likely look beyond accuracy and possibly use precision (the ratio of correct positive classifications), recall (the ratio of correctly predicted positive classifications), F1 score and other measures to determine how well a model performs.

 Confusion matrix

A confusion matrix shows the number of true positives, false negatives, false positives and true negatives. As you can see in the matrix below, from 1 run, the model produced a fair number of true positives (e.g.: the model predicted that the customer would not default and the true result was that the customer did not default). These results were mostly consistent across all 10 runs.

Feature Importances

Below is a table from 1 run of the model that shows the feature names and their corresponding feature importances. Feature importance ranks how well the explanatory variables can predict the response variable, in this case the possibility of credit card default. Feature importances always sum to 1, thus individual feature importances can be seen as a percent of their importance towards predicting the response variable. More information on how the Python feature importances are calculated can be found here.

Max trees

In the chart below we can see that only 7-8 trees were needed to provide ~80% prediction accuracy mentioned above.

To produce the results below a random forest was created 25 times with the same training and prediction data. This process started with 1 tree. The number of trees increased by 1 until 25 trees were reached.

Appendix - Python code

===== import pandas as pd import numpy as np import os import matplotlib.pylab as plt from sklearn.cross_validation import train_test_split import sklearn.metrics

# Feature Importance from sklearn.ensemble import ExtraTreesClassifier from sklearn.ensemble import RandomForestClassifier

# --------------- # Set constants

# This is where the classification tree image will be placed. os.chdir("C:\Python\Coursera-DataMan")

# Set numerical display formats. pd.options.display.float_format = '{:0,.3f}'.format

# Set the number of trees to create MAXTREES = 25

# --------------- # Load the dataset, Drop NA's and show overall dataframe summary. # defaultRateDF = pd.read_csv("default_rate_credit_card_clients.csv")

defaultRateDF = defaultRateDF.dropna()

# The data is already INT64 so no need to convert to_numeric #print('\nDefaultRateDF data types\n', defaultRateDF.dtypes) print ("\nDataframe describe results\n", defaultRateDF.describe())

# ------------------------- # Modeling and Prediction # Predictors are explanatory vars. # Target is the response variable. # Split into training and testing sets. predictors = defaultRateDF[['X1_LIMIT_BAL', 'X2_SEX', 'X3_EDUCATION', 'X4_MARRIAGE', 'X5_AGE', 'X6_PAY_0', 'X7_PAY_2', 'X8_PAY_3', 'X9_PAY_4', 'X10_PAY_5', 'X11_PAY_6', 'X12_BILL_AMT1', 'X13_BILL_AMT2',  'X14_BILL_AMT3', 'X15_BILL_AMT4', 'X16_BILL_AMT5',   'X17_BILL_AMT6', 'X18_PAY_AMT1', 'X19_PAY_AMT2',   'X20_PAY_AMT3', 'X21_PAY_AMT4', 'X22_PAY_AMT5',   'X23_PAY_AMT6']]

targets = defaultRateDF.Y_default_payment_next_month

pred_train, pred_test, tar_train, tar_test  = train_test_split(predictors, targets, test_size=.4)

# For kicks, show the shape of the training and test sets. print("\npred_train.shape:", pred_train.shape) print("\npred_test.shape:", pred_test.shape) print("\ntar_train.shape:", tar_train.shape) print("\ntar_test.shape:", tar_test.shape)

# --------------- # Build model on training data and show results. # N_estimators is the number of trees to build with the random forest algorithm. classifier=RandomForestClassifier(n_estimators = MAXTREES) classifier=classifier.fit(pred_train, tar_train)

# Use the test data on the previously created classifier to predict the outcome. predictions=classifier.predict(pred_test)

# Get and print the following metrics. accuScore = sklearn.metrics.accuracy_score(tar_test, predictions) print("\nAccuracy score %.3f\n" % accuScore)

confMatrix = sklearn.metrics.confusion_matrix(tar_test, predictions) print("\nConfusion Matrix") print("----------------") print("True Positive   %i | False Negative %i" % (confMatrix[0,0], confMatrix[0,1])) print("False Positive  %i | True Negative  %i" % (confMatrix[1,0], confMatrix[1,1])) #print(confMatrix)

# Visually show the confusion matrix. plt.matshow(confMatrix) plt.title('Confusion matrix') plt.ylabel("True Label") plt.xlabel("Predicted Label") plt.colorbar() plt.show()

# --------------- # Fit an Extra Trees model to the data. # Since we don't interpret the individual trees, the most helpful info is the # measured importance for each explanatory variable based on how many votes or # splits each tree has produced in our 25 tree ensemble. # The following generates these importance scores. model = ExtraTreesClassifier() model.fit(pred_train, tar_train)

# Display the relative importance of each attribute #print("Raw Feature Importances") #print(model.feature_importances_) #print(model.feature_importances_[0])

# Create a DF with Feature names and associated feature importances so we can # sort the output and clearly label the features. FeatureNames = ['X1_LIMIT_BAL', 'X2_SEX', 'X3_EDUCATION', 'X4_MARRIAGE',                            'X5_AGE', 'X6_PAY_0', 'X7_PAY_2', 'X8_PAY_3', 'X9_PAY_4',                            'X10_PAY_5', 'X11_PAY_6', 'X12_BILL_AMT1', 'X13_BILL_AMT2',                            'X14_BILL_AMT3', 'X15_BILL_AMT4', 'X16_BILL_AMT5',                            'X17_BILL_AMT6', 'X18_PAY_AMT1', 'X19_PAY_AMT2',                            'X20_PAY_AMT3', 'X21_PAY_AMT4', 'X22_PAY_AMT5',                            'X23_PAY_AMT6']

FeatureImportancesDF = pd.DataFrame(data=FeatureNames, columns=['FeatureName']) FeatureImportancesDF['FeatureImportances'] = Series(model.feature_importances_, index=FeatureImportancesDF.index) FeatureImportancesDF = FeatureImportancesDF.sort_values(by='FeatureImportances', ascending=False) print("\n      Features and Importances (high to low)\n", FeatureImportancesDF)

# ---------------

# Run a different number of trees and see the effect of that on the accuracy # of the prediction. Essentially did we need MAXTREES to get the % # or could we have used fewer trees? trees = range(MAXTREES) accuracy = np.zeros(MAXTREES)

for idx in range(len(trees)):   classifier=RandomForestClassifier(n_estimators = idx + 1)   classifier=classifier.fit(pred_train, tar_train)   predictions=classifier.predict(pred_test)   accuracy[idx]=sklearn.metrics.accuracy_score(tar_test, predictions)

plt.cla() plt.plot(trees, accuracy)

=====

[mlfanatic]

Machine Learning for Data Analysis - Week 1 assignment - Create a classification tree

Executive summary

This week we'll focus on classification trees/decision trees. Classification trees allow exploration of non-linear explanatory variables.

Also, this week we will forego hypothesis testing and instead separate the data into training and test data sets in order to test the ability of a classification tree to predict an outcome. I used the same data that I described in my prior post here.

This week we'll focus on the 2 explanatory variables from last weeks experiment. Last weeks experiment showed that marriage and education levels had high significance and the highest odd ratios towards predicting credit card default. Specifically:

Binary explanatory variable: SingleMarried (0 = single, 1 = married)

Binary explanatory variable: Education (0 = high school, 1 = university or graduate school; )

Binary response variable: credit card payment default (0 = no default, 1 = default)

Results

Figure 1. Classification tree (N = 29,163) for credit card default (response variable)

I manually added the “Node” numbers above to make it easier to explain the results.

Also, the term "sample" refers to the training data set that was used to create the classification tree above. 40% of the data set was used to create the training data set.

Decision tree analysis was performed to test nonlinear relationships among the explanatory variables and a binary response variable previously mentioned.

Marital status and education level explanatory variables were included in a classification tree model to evaluate predicting the credit card default response variable.

Node 1 - Education level was the first variable to separate the sample into two subgroups; 2,875 customers in our sample have achieved a high school education (Node 2) and 14,622 customers have received a university or graduate level education (Node 3).

Node 2 - Of the customers in our sample who achieved a high school education 74% were single (2136 / 2875 = 0.743) and 26% were married  (739 / 2875 = 0.257).

Next, another split is made on marital status; customers who had high school education and were single (Node 4), and customers who had a high school education and were married (Node 5).

Node 4 - Of these single customers 861 have not defaulted on their credit card debt, whereas 267 have defaulted.

Node 5 - Of these married customers 1,275 have not defaulted on their credit card debt whereas 472 have defaulted.

Node 3 - On the right side of the tree are customers who have a university or graduate level education. Of these customers 78% were single (11,419 / 14,622 = 0.7809) and 22% were married (3,203 / 14,622 = 0.2191).

Next, another split is made on marital status; customers who have a university or graduate level education and were single (Node 6) and customers with this education level who were married (Node 7).

Node 6 - Of these single customers 6,526 have not defaulted on their credit card debt whereas 1,751 have defaulted.  

Node 7 - Of these married customers 4,893 have not defaulted on their credit card debt, however 1,452 have defaulted.

To summarize 1 pathway we could follow through nodes 1, 3, and 6 and say that 6,526 customers in our sample who have a university or graduate level education and were single did not default on their credit card debt.

Accuracy is the ratio of instances the classifier correctly classified. The model was run 10 times and consistently achieved ~77% accuracy. However, accuracy best represents how well a classifier performed when using balanced classes/features/columns like education or marital status above.

Overall, our data is imbalanced. 16% of the customers in our original data set achieved a high school degree while 83% received a university or graduate level education. Likewise, 46% were single and 54% were married.

In the future we will likely look into using precision (the ratio of correct positive classifications), recall (the ratio of correctly predicted positive classifications), F1 score and other measures to determine how well a model performs.

Given the 2 explanatory variables above, marital status was a slightly better predictor of credit card default when compared to education level. The explanatory variable importance predictors were:

Education          0.40

SingleMarried   0.60

Appendix - Python code

===== import pandas as pd #import numpy as np import os from sklearn.cross_validation import train_test_split from sklearn.tree import DecisionTreeClassifier import sklearn.metrics from sklearn import tree from io import StringIO import pydotplus

# This is where the classification tree image will be placed. os.chdir("C:\Apps\Docs\Python\Coursera-DataMan")

# Set numerical display formats. pd.options.display.float_format = '{:0,.2f}'.format

# ----------------- # Load the dataset defaultRateDF = pd.read_csv("default_rate_credit_card_clients.csv")

defaultRateDF = defaultRateDF.dropna()

# Marriage and education are explanatory variables (X) # Payment default is the binary response variable (Y) defaultRateDF['X3_EDUCATION'] = pd.to_numeric(defaultRateDF['X3_EDUCATION'], errors='coerce') defaultRateDF['X4_MARRIAGE'] = pd.to_numeric(defaultRateDF['X4_MARRIAGE'], errors='coerce') defaultRateDF['Y_default_payment_next_month'] = pd.to_numeric(defaultRateDF['Y_default_payment_next_month'], errors='coerce')

# Create a new dataframe and add the columns that I intend to use. marriageCreditDF = defaultRateDF[['X3_EDUCATION', 'X4_MARRIAGE', 'Y_default_payment_next_month']]

# Rename the column names to be more readable. marriageCreditDF = marriageCreditDF.rename(columns={'X3_EDUCATION': 'Education', 'X4_MARRIAGE': 'SingleMarried', 'Y_default_payment_next_month': 'Y_default'})

# ----------------- # FILTER DATA SET # Filter the data for marriage = 1 or 2 ONLY (married or single) marriageCreditDF = marriageCreditDF[(marriageCreditDF['SingleMarried'] == 1) | (marriageCreditDF['SingleMarried'] == 2)]

# Filter the data for education = 1 through 3. # Other educational values aren't defined so can't really include them. marriageCreditDF = marriageCreditDF[(marriageCreditDF['Education'] == 1) | (marriageCreditDF['Education'] == 2) | (marriageCreditDF['Education'] == 3)]

# --------------- # RECODE the data set # Recode Education. # Original values: Education (1 = graduate school; 2 = university; 3 = high school # Recode the education variable to: # 1 = university and graduate school; 0 = high school; recode1 = {1: 1, 2: 1, 3: 0} marriageCreditDF['Education'] = marriageCreditDF['Education'].map(recode1)

# Recode the marital status variable to: 0 = single and 1 = married. recode1 = {2: 0, 1:1} marriageCreditDF['SingleMarried'] = marriageCreditDF['SingleMarried'].map(recode1)

# Optional summary of the data we're going to use. print ('\nDefault rate data information from our data set:\n', marriageCreditDF.describe())

# ------------------------- # Modeling and Prediction

# Y_default is from the CSV file above. This is our response variable. targets = marriageCreditDF.Y_default

# This is the subset data columns. predictors = marriageCreditDF[['Education', 'SingleMarried']]

#Split into training and testing sets. predictor_train, predictor_test, target_train, target_test = train_test_split(predictors, targets, test_size=.4)

predictor_train.shape predictor_test.shape target_train.shape target_test.shape

# Build model from the training data classifier = DecisionTreeClassifier() classifier = classifier.fit(predictor_train, target_train)

predictions = classifier.predict(predictor_test)

# ----------------------- # Use common metrics to determine how well the classification model performed.

important_predictors = DataFrame(classifier.feature_importances_, columns = ['Imp'], index = predictors.columns).Imp print ('\nExplanatory variable importance predictors:\n', important_predictors)

accuracyScore = sklearn.metrics.accuracy_score(target_test, predictions) print('\nAccuracy Score: %.3f \n' % accuracyScore )

# ------------------------ # Print the classification tree print('\nGenerating decision tree png...\n') out = StringIO() tree.export_graphviz(classifier, out_file=out, feature_names=predictors.columns) graph = pydotplus.graph_from_dot_data(out.getvalue()).write_png('ClassificationTree1.png')

=====

[mlfanatic]

Course: Regression Modeling in Practice - Week 4 assignment - Test a logistic regression model

Executive summary

The assignment for this week is to test a logistic regression model. I used the same data that I described in my prior post here.

In our multiple regression experiment last week we found that after adjusting for potential confounding factors (gender and education level) education level (Beta=55,761, p=.0001) was significantly and positively associated with a larger credit limit. Secondarily, marital status was also significantly associated with a larger credit limit (Beta= 16,854, p=.0001). However, even though the results above allowed us to reject last weeks null hypothesis the inaccurate model did not support rejecting the null hypothesis.

This week we will create a logistic regression model. We will add a new (and exciting!) response variable to determine if customers in the sample have defaulted on a credit card payment. We will continue to use the explanatory variables education level, marital status, gender and customer age to determine the odds of default. Let's begin!

Experiment

Overall, I hypothesis that married customers in our sample will default significantly less on their credit card payments than credit card customers who are single. I hypothesize that marriage, and likely family, bring personal responsibility which may also be present in fulfilling credit card payments.

Null Hypothesis: There is no significant difference in credit card default rates between married customers and single customers.

Alternate Hypothesis: Married customers will default significantly less on credit card payments when compared to single customers.

The alpha chosen for this experiment is 0.05. Alpha is the probability of committing a Type 1 error or rejecting the null hypothesis when it should not have been rejected.

While running this experiment we will also determine if education, gender and age confound (influence) the marital status variable.

  Details behind the data

After filtering we have data from 29,163 individual customers to use in our experiment. The reason for filtering can be found in my prior blog post here in the details section.

Overall, our sample data has 6,514 customers who default, or 22.34%.

In the logistic regression model that we will create below, we need to provide explanatory variables and a response variable that are either 0 or 1.

Education had the original values; 1 = graduate school; 2 = university; 3 = high school. However, these have been recoded to 2 binary groups; 1 = higher education (university and graduate school) or 0 = high school.

Marital status has already been coded as 1 = married, 0 = single

Gender has already been coded as 1 = Female, 0 = Male.

Since customer age represented a range of values from 21 to 79 with 34 being the median, I recoded age as 1 = less than 35, 0 = 35 and older. A new variable named 'CustomerAge_34median' was created to hold this binary value. This new variable holds a nice split in the data:

15,386 customers were less than 35 years of age.

13,777 customer were aged 35 and older.

The response variable 'Y_default' represents whether or not the customer defaulted on their credit card payment; 1 = Yes, 0 = No.

Below are initial results that show the probability of default among single customers is a little lower than married customers (21% vs. 24%).

Final experiment results

As you’ll see in the code in the appendix the experiment was run 4 times. Each time we added an explanatory variable; SingleMarried, Education, Gender and CustomerAge_34median. The final results for the last experiment are in the image below.

After adjusting for potential confounding variables (education, gender and age), the odds of married customers defaulting on credit card debt was only 1.17 times more likely than single customers (OR = 1.17, CI=1.10 to 1.25, P=.0001). Since our p-value of 0.0001 is less than 0.05 alpha we can reject the null hypothesis.

Results for additional explanatory variables include:

Customers in the sample who achieved a university or graduate school education were 0.84 times less likely to default than high-school educated customers (OR=0.84, CI=0.78 to 0.90, P=0.0001). In other words, the probability of default is slightly lower among those with university or graduate school education than those with a high-school education. Given the low p-value this result is significant.

Female customers in the sample were 0.82 times less likely to default than men (OR=0.82, CI=0.78 to 0.87, P=0.0001). This is equivalent to saying that the probability of default is slightly lower among women than among men. Given the low p-value this result is significant.

Customers who were less than 35 years of age were only 1.03 times more likely to default as customers who were 35 and older (OR=1.03, CI=0.97 to 1.10, P=0.293). Given the high p-value this result is not significant.

Also, between the 4 variables above I didn't find any strong confounding variables. For example, when marital status was tested alone it received an odds ratio of 1.16. Next, when the education explanatory variable was added the marital status odds ratio decreased to 1.14, a decrease of only 0.02. Similarly, marital status changed by +/- 0.02 when gender and the median customer age were added sequentially. 

In addition, as I added explanatory variables the prior p-values for marital status remained at 0.000 which is another sign that confounding didn’t exist in this experiment.

Beyond the odds ratio, the plot below shows the probability of default for each of the 4 explanatory variables.

For example:

There are 15,738 total single customers, 3,325 of them default. 3,325 / 15,738 = 0.21 or 21% chance of default.

There are 13,425 total married customers, 3,189 of them default. 3,189 / 13,425 = 0.24 or 24% chance of default.

The other variables in the chart above should be easily understood given this example.

Overall, the odds results are fairly small, nothing special... The results above imply that to more accurately predict credit card payment default:

other variables in the sample data, such as payment history, may represent a stronger relationship.

there may be additional data outside of the current sample that can better predict default.

  Appendix - Python code

import pandas as pandas import statsmodels.api as sm import statsmodels.formula.api as smf import matplotlib.pyplot as plt import seaborn as sns import numpy

# ------------------------------------ # Function to group customers into 2 age groups. # 1 = customers age 34 and under. 0 = all other ages 35 and over.

def CustomerAge_34median (row):   if (row['CustomerAge'] < 35):      return 1   else:      return 0

# ------------------------------------

#Read my local data into a data frame (DF). defaultRateDF = pandas.read_csv('default_rate_credit_card_clients.csv', low_memory=False)

# Set numerical display formats. pandas.options.display.float_format = '{:0,.2f}'.format

# convert variables I'm using to numeric format. # Marriage and customer age are explanatory variables (X) # Credit Limit is the quantitative response variable (Y) defaultRateDF['X4_MARRIAGE'] = pandas.to_numeric(defaultRateDF['X4_MARRIAGE'], errors='coerce') defaultRateDF['X5_AGE'] = pandas.to_numeric(defaultRateDF['X5_AGE'], errors='coerce') defaultRateDF['X1_LIMIT_BAL'] = pandas.to_numeric(defaultRateDF['X1_LIMIT_BAL'], errors='coerce') defaultRateDF['X2_SEX'] = pandas.to_numeric(defaultRateDF['X2_SEX'], errors='coerce') defaultRateDF['X3_EDUCATION'] = pandas.to_numeric(defaultRateDF['X3_EDUCATION'], errors='coerce') defaultRateDF['Y_default_payment_next_month'] = pandas.to_numeric(defaultRateDF['Y_default_payment_next_month'], errors='coerce')

# Create a new dataframe and add the columns that I intend to use. marriageCreditDF = defaultRateDF[['X2_SEX', 'X3_EDUCATION', 'X4_MARRIAGE', 'X5_AGE', 'X1_LIMIT_BAL', 'Y_default_payment_next_month']]

# Rename the column names to be more readable. marriageCreditDF = marriageCreditDF.rename(columns={'X2_SEX': 'Gender', 'X3_EDUCATION': 'Education', 'X4_MARRIAGE': 'SingleMarried', 'X5_AGE': 'CustomerAge', 'X1_LIMIT_BAL': 'CreditLimit', 'Y_default_payment_next_month': 'Y_default'})

# --- FILTER --- # Filter the data for marriage = 1 or 2 ONLY (married or single) marriageCreditDF = marriageCreditDF[(marriageCreditDF['SingleMarried'] == 1) | (marriageCreditDF['SingleMarried'] == 2)]

# Filter the data for education = 1 through 3. # Other educational values aren't defined so can't really include them. marriageCreditDF = marriageCreditDF[(marriageCreditDF['Education'] == 1) | (marriageCreditDF['Education'] == 2) | (marriageCreditDF['Education'] == 3)]

# --- RECODE --- # Recode Education. # Original values: Education (1 = graduate school; 2 = university; 3 = high school # Recode the education variable to: # 1 = university and graduate school; 0 = high school; recode1 = {1: 1, 2: 1, 3: 0} marriageCreditDF['Education'] = marriageCreditDF['Education'].map(recode1)

# Recode the marital status variable to: 0 = single and 1 = married. recode1 = {2: 0, 1:1} marriageCreditDF['SingleMarried'] = marriageCreditDF['SingleMarried'].map(recode1)

# Recode the gender variable to: 0 = Male and 1 = Female. recode1 = {1: 0, 2: 1} marriageCreditDF['Gender'] = marriageCreditDF['Gender'].map(recode1)

# Recode Customer Age print ('\nCustomer Age details:\n', marriageCreditDF["CustomerAge"].describe()) print ('\nCustomer median age: %.2f\n' % marriageCreditDF["CustomerAge"].median())

# The median customer age is 34. The mean customer age is 35. For Logistic regression # need to create 2 categories; 1 = customers age 34 and under. Everyone else is 0.

marriageCreditDF['CustomerAge_34median'] = marriageCreditDF.apply (lambda row: CustomerAge_34median (row), axis=1)

chk2 = marriageCreditDF['CustomerAge_34median'].value_counts(sort=False, dropna=False) print("\nCustomerAge_34median count. 0 = 35 and older, 1 = less than 35: \n", chk2)

print('\n###############################')

print(' LOGISTIC REGRESSION EXPERIMENT') print('################################')

# ------------------------------- # logistic regression with SingleMarried lreg1 = smf.logit(formula = 'Y_default ~ SingleMarried', data = marriageCreditDF).fit() print (lreg1.summary())

# odds ratios print("\nOdds Ratios") print(numpy.exp(lreg1.params))

#print(numpy.exp(lreg1.params[0])) #print(numpy.exp(lreg1.params[1]))

# odd ratios with 95% confidence intervals params = lreg1.params conf = lreg1.conf_int() conf['Odds Ratio'] = params conf.columns = ['Lower CI', 'Upper CI', 'Odds Ratio'] print("\n", numpy.exp(conf)) print('')

# ------------------------------- # logistic regression with marital status and Education lreg2 = smf.logit(formula = 'Y_default ~  SingleMarried + Education', data = marriageCreditDF).fit() print (lreg2.summary())

# odds ratios print("\nOdds Ratios") print(numpy.exp(lreg2.params))

# odd ratios with 95% confidence intervals params = lreg2.params conf = lreg2.conf_int() conf['Odds Ratio'] = params conf.columns = ['Lower CI', 'Upper CI', 'Odds Ratio'] print("\n", numpy.exp(conf)) print('')

# ------------------------------- # logistic regression with marital status, Education and Gender. lreg3 = smf.logit(formula = 'Y_default ~ SingleMarried + Education + Gender', data = marriageCreditDF).fit() print (lreg3.summary())

# odds ratios print("\nOdds Ratios") print(numpy.exp(lreg3.params))

# odd ratios with 95% confidence intervals params = lreg3.params conf = lreg3.conf_int() conf['Odds Ratio'] = params conf.columns = ['Lower CI', 'Upper CI', 'Odds Ratio'] print("\n", numpy.exp(conf)) print('')

# ------------------------------- # logistic regression with marital status, Education, Gender, customer age. lreg4 = smf.logit(formula = 'Y_default ~ SingleMarried + Education + Gender + CustomerAge_34median', data = marriageCreditDF).fit() print (lreg4.summary())

# odds ratios print("\nOdds Ratios") print(numpy.exp(lreg4.params))

# odd ratios with 95% confidence intervals params = lreg4.params conf = lreg4.conf_int() conf['Odds Ratio'] = params conf.columns = ['Lower CI', 'Upper CI', 'Odds Ratio'] print("\n", numpy.exp(conf)) print('')

# ------------------------------ # Visualize the data overall.

# 2 bars scat1 = sns.barplot(x="SingleMarried", y="Y_default", data=marriageCreditDF) plt.show(scat1)

# Pair plot the 4 variables. This shows the probability of default for each # of the 4 explanatory variables in our experiment.  

sns.set(style="whitegrid") PairGrid1 = sns.PairGrid(marriageCreditDF, y_vars="Y_default",                 x_vars=['Education', 'SingleMarried', 'Gender', 'CustomerAge_34median'],                 size=5, aspect=.5)

PairGrid1.map(sns.pointplot, color=sns.xkcd_rgb["blue"]) PairGrid1.set(ylim=(0, 1)) sns.despine(fig=PairGrid1.fig, left=True)

[mlfanatic]

Course: Regression Modeling in Practice - Week 3 assignment - Test a multiple regression model

Executive summary

The assignment for this week is to test a multiple regression model. I used the same data described in my prior post here.

In our prior experiment we were able to reject the null hypothesis and show that there is a significant difference in the credit limit extended to Taiwanese customers who are married vs. those who are single. Remember this does not mean we "accept" the alternate hypothesis or the alternative hypothesis is "true". It just means that we have support for our alternate hypothesis.

This week we will include customer age to our experiment so that a multiple regression model can be tested. We will also add 2 confounding variables.

Experiment

Overall, I hypothesize that older married customers in our sample will have more credit allocated than older single customers. My rationale is that more credit is needed by older, married couples to buy shared items for the household.

Null Hypothesis: there is no significant difference in the credit limit extended to Taiwanese customers who are married or single regardless of age.

Alternate Hypothesis: more credit is allocated to older Taiwanese customers who are married vs. those who are single.

The alternate hypothesis states that just one of the independent explanatory variables, marital status or age, will be useful to predict the response variable, credit limit.

The alpha chosen for this experiment is 0.05. Alpha is the probability of committing a Type 1 error or rejecting the null hypothesis when it should not have been rejected.

Details behind the data:

Overall, the experiment used data for 29,163 customers.

Please see my prior post for information behind the credit limit response variable and marital status explanatory variable.

Age is a new explanatory variable added to test the hypothesis above.

In addition, we'll look at possible confounding explanatory variables gender and education level (high school, university and graduate level).

Lastly, all "dollar" results mentioned below are in Taiwanese dollars.

Hypothesis Results:

Based on the information above we can see that there is a slight positive relationship between customer age and the amount of credit extended to the customers in our sample. Below we will further investigate this potential relationship.

The areas highlighted on the image above are discussed below.

We'll focus on adjusted R-squared because "Adjusted R-squared explicitly adjusts for the sample size n and the number of explanatory variables. It is common to use adjusted R-squared for model selection because it imposes a penalty for any additional explanatory variable that is included in the analysis." [Source: Business Statistics page 446].

The adjusted R-squared results from our experiment are: 0.024 or 2.4%. This means that only 2.4% of the response variable, credit limit, is explained by the explanatory variables marital status and customer age. Conversely, this means that 97.6% of the affect on the response variable is not described by the model presented here. Not a very good result so far.

The Prob (F-statistic) is 2.12e-158, this is our p-value. Since our p-value is well below our alpha of 0.05 we could reject the null hypothesis above, however further results are provided below before we make a final decision.

Hypothesis regression coefficient results

The intercept is 173,344.19 dollars and shows the credit limit (response variable) when both marriage and customer age equal zero.

The MarriedSingle coefficient  is -9,352.63. The model estimates that if marital status increases by 1 unit (from Married = 0 to Single = 1) then it is predicted that the customer will have 9,352 less credit, holding all other explanatory variables constant. Another way to state this coefficient is that married customers could have on average 9,352 more credit than single customers.

The CustomerAge_centered coefficient is 1,919.04. The model estimates that if customer age increases by 1 unit then credit allocated to the customer is predicted to increase by 1,919 dollars, holding all other explanatory variables constant.

Hypothesis coefficient P-values

MarriedSingle and CustomerAge_centered p-values are zero which means that each significantly contributes to the ability to predict the credit limit, however further results are provided below before we make a final decision.

Hypothesis 95% confidence interval

Overall, confidence intervals show a range of scores whereby the population mean will fall in that range in 95% of the samples.

The MarriedSingle confidence interval ranges from -12,700 to -6,041. With 95% confidence we can say that the population mean amount of credit provided to single customers is between -12,700.00 to -6,041 less than people who are married.

The CustomerAge_centered confidence interval ranges from 1,739.29 to 2,098.78. With 95% confidence we can say that the population mean amount of credit provided to customers will increase between 1,739 to 2,098 dollars for each year they advance in age.

Hypothesis plots to evaluate model fit

The quantile to quantile (Q-Q) plot above shows that the data are not normally distributed. If the data were normally distributed then the blue dots (our data) would fall on the red line. This is an indication that our model is not sufficient to represent the data.

Note that to produce the plot above I used CustomerAge that was not centered. This was necessary to produce a readable plot. CustomerAge_Centered had too many values too close to each other to produce a readable plot.

The standardized residuals plot shows the difference between the predicted and outcome values represented as a standard deviation. The red horizontal line at zero is the sample mean. It's best to try and visualize this chart like a normal distribution tiled on its right-side.

Our plot above shows that we seem to have a high number of residuals beyond 3 standard deviations from the mean (3+ on the Y-axis above). This means that the level of error in our model not acceptable and the model is a poor fit of the sample data.

Nice plot! :-)

A leverage\influence plot above is used to identify observations that are outliers and have a large influence on model estimation. The leverage plot isn't useful when you get more than a few hundred data points. My sample had just over 28,000 data points.

For this assignment we can say that I know how to produce an influence/leverage plot however there was too much data for it to be useful.

Results for confounding variables

While performing the experiment above I also investigated confounding variables education level and gender.There was evidence of confounding as described below.

Education level had the largest effect on credit of all explanatory variables tested, well beyond age as originally hypothesized above. Each unit increase in education level (high-school to university to graduate school) resulted in 55,760 dollars of additional credit.

Gender was the third most influential explanatory variable. Women in the study had 11,060 more credit allocated than men.

A summary of coefficients can be found below.

Explanatory variable: Coefficients showing credit allocated per unit increase

education:          55,761

married single*: 16,854

gender:              11,060

customer age:    2,580

* Note: I made this number positive such that the variable shows marital status of single customers followed by married customers to make the coefficient sorting easier to compare.

Conclusion

After adjusting for potential confounding factors (gender and education level) education level (Beta=55,761, p=.0001) was significantly and positively associated with a larger credit limit. Secondarily, marital status was also significantly associated with a larger credit limit (Beta= 16,854, p=.0001).

In addition, even though the p-value and coefficient p-values above allows us to reject the null hypothesis the adjusted R-squared, Q-Q plot and standardized residuals plot show that our multiple regression model is not sufficient to represent the data. The inaccurate model above does not support rejecting the null hypothesis.

Appendix - Python code

 =====

import pandas as pandas import statsmodels.api as sm import statsmodels.formula.api as smf import matplotlib.pyplot as plt import seaborn

#Read my local data into a data frame (DF). defaultRateDF = pandas.read_csv('default_rate_credit_card_clients.csv', low_memory=False)

# Set numerical display formats. pandas.options.display.float_format = '{:0,.2f}'.format

# convert variables I'm using to numeric format. # Marriage and customer age are explanatory variables (X) # Credit Limit is the quantitative response variable (Y) defaultRateDF['X4_MARRIAGE'] = pandas.to_numeric(defaultRateDF['X4_MARRIAGE'], errors='coerce') defaultRateDF['X5_AGE'] = pandas.to_numeric(defaultRateDF['X5_AGE'], errors='coerce') defaultRateDF['X1_LIMIT_BAL'] = pandas.to_numeric(defaultRateDF['X1_LIMIT_BAL'], errors='coerce') defaultRateDF['X2_SEX'] = pandas.to_numeric(defaultRateDF['X2_SEX'], errors='coerce') defaultRateDF['X3_EDUCATION'] = pandas.to_numeric(defaultRateDF['X3_EDUCATION'], errors='coerce')

# Create a new dataframe and add the columns that I intend to use. marriageCreditDF = defaultRateDF[['X2_SEX', 'X3_EDUCATION', 'X4_MARRIAGE', 'X5_AGE', 'X1_LIMIT_BAL']]

# Rename the column names to be more readable. marriageCreditDF = marriageCreditDF.rename(columns={'X2_SEX': 'Gender', 'X3_EDUCATION': 'Education', 'X4_MARRIAGE': 'MarriedSingle', 'X5_AGE': 'CustomerAge', 'X1_LIMIT_BAL': 'CreditLimit'})

# --- FILTER --- # Filter the data for marriage = 1 or 2 ONLY (married or single) marriageCreditDF = marriageCreditDF[(marriageCreditDF['MarriedSingle'] == 1) | (marriageCreditDF['MarriedSingle'] == 2)]

# Filter the data for education = 1 through 3. # Other educational values aren't defined so can't really include them. marriageCreditDF = marriageCreditDF[(marriageCreditDF['Education'] == 1) | (marriageCreditDF['Education'] == 2) | (marriageCreditDF['Education'] == 3)]

# --- RECODE --- # Recode the marital status variable to: 0 = married and 1 = single. recode1 = {1: 0, 2: 1} marriageCreditDF['MarriedSingle'] = marriageCreditDF['MarriedSingle'].map(recode1)

# Recode the gender variable to: 0 = Male and 1 = Female. recode1 = {1: 0, 2: 1} marriageCreditDF['Gender'] = marriageCreditDF['Gender'].map(recode1)

# Recode the education variable to: # 2 = graduate school; 1 = university; 0 = high school; recode1 = {1: 10, 2: 20, 3: 30} marriageCreditDF['Education'] = marriageCreditDF['Education'].map(recode1) recode1 = {10: 2, 20: 1, 30: 0} marriageCreditDF['Education'] = marriageCreditDF['Education'].map(recode1)

# Center quantitative explanatory variable for regression analysis

# and verify it has been centered. marriageCreditDF['CustomerAge_centered'] = (marriageCreditDF['CustomerAge'] - marriageCreditDF['CustomerAge'].mean()) print ('Verify the Customer Age has been centered, the mean should be zero:\n', marriageCreditDF["CustomerAge_centered"].describe())

print("\nCredit limit details\n", marriageCreditDF['CreditLimit'].describe())

# -------------------------------- # Simple linear regression for comparison

reg0 = smf.ols('CreditLimit ~ MarriedSingle', data=marriageCreditDF).fit() print ("\nSimple OLS regression model for the association between marital status and credit limit") print (reg0.summary())

# --------------------------------- # Multiple Regression #1 - X vars = MarriedSingle + CustomerAge_centered

reg1 = smf.ols('CreditLimit ~ MarriedSingle + CustomerAge_centered', data=marriageCreditDF).fit() #print (reg1.summary())

print ("\nOLS regression model for the association between marital status, customer age and credit limit") # Show the important sections of the OLS Regression Results table. print(reg1.summary().tables[0]) print(reg1.summary().tables[1]) print ("\nEasier to read regression coefficient results:\n", reg1.params)

# --------------------------------- # Multiple Regression #2 - X vars = MarriedSingle + CustomerAge_centered + Gender

reg2 = smf.ols('CreditLimit ~ MarriedSingle + CustomerAge_centered + Gender', data=marriageCreditDF).fit() #print (reg1.summary())

print ("\nOLS regression model for the association between marital status, customer age, gender and credit limit") # Show the important sections of the OLS Regression Results table. print(reg2.summary().tables[0]) print(reg2.summary().tables[1]) print ("\nEasier to read regression coefficient results:\n", reg2.params)

# --------------------------------- # Multiple Regression #3 - X vars = MarriedSingle + CustomerAge_centered + Gender + Education level

reg3 = smf.ols('CreditLimit ~ MarriedSingle + CustomerAge_centered + Gender + Education', data=marriageCreditDF).fit() #print (reg1.summary())

print ("\nOLS regression model for the association between marital status, customer age, gender, education level and credit limit") # Show the important sections of the OLS Regression Results table. print(reg3.summary().tables[0]) print(reg3.summary().tables[1]) print ("\nEasier to read regression coefficient results:\n", reg3.params)

# Charts for the hypothesis only. scat1 = seaborn.regplot(x="CustomerAge", y="CreditLimit", scatter=True, data=marriageCreditDF) plt.xlabel('CustomerAge') plt.ylabel('CreditLimit') plt.title ('Scatterplot for the Association Between Customer Age and Credit Limit') plt.show(scat1)

# Q-Q plot for normality # NOTE, this won't appear unless you comment out all remaining code below. QQPlot=sm.qqplot(reg1.resid, line='r')

# Can't use the centered age since the resulting plot is unuseable. # Instead used the CustomerAge that is NOT centered just for the residual plot. reg2 = smf.ols('CreditLimit ~ MarriedSingle + CustomerAge', data=marriageCreditDF).fit()

# Standardized residuals for all observations. stdres=pandas.DataFrame(reg2.resid_pearson) plt.plot(stdres, 'o', ls='None') l = plt.axhline(y=0, color='r') plt.title('Standardized residuals for all observations') plt.ylabel('Standardized Residual') plt.xlabel('Observation Number')

# Additional regression diagnostic plots # Not required for homework assignment. RegDiagPlots = plt.figure(figsize=(12,8)) RegDiagPlots = sm.graphics.plot_regress_exog(reg1,  "CustomerAge_centered", fig=RegDiagPlots)

# leverage plot / influence plot LeveragePlot=sm.graphics.influence_plot(reg1, size=8) print(LeveragePlot)

 =====

Notable resources:

While doing research for this assignment I came across the following resources that were valuable.

This page describes the COMMON MISTAKES IN USING STATISTICS: Spotting and Avoiding Them.

This video describes how to interpret multiple regression results, confidence intervals, etc.

[mlfanatic]

Course: Regression Modeling in Practice - Week 2 assignment - Test a basic linear regression model

Executive summary

The assignment for this week is to test a linear regression model. I used the same data that I described in my prior post here. Overall, married customers in our sample have more credit allocated than customers who are single.

Experiment

Taiwanese customers in our sample who are married will likely have more credit allocated than customers who are single. My rationale is that more credit is needed by married couples to buy shared items for the household.

Null Hypothesis: there is no significant difference in the credit limit extended to Taiwanese customers who are married vs. those who are not married. 

Alternate Hypothesis: there is a significant difference in the credit limit extended to Taiwanese customers who are married vs. those who are not married.  

The alpha chosen for this experiment is 0.05.

Details

Marital status is a categorical explanatory variable as required by this assignment. This variable contains original status values: 0 = unknown; 1 = married; 2 = single; 3 = others. For this assignment I will subset the data to select observations from 2 categories; married and single. Removing the other marriage categories only removes 377 rows or 1.3% of the data.

Also, per the assignment instructions the original marital status variable has been recoded to; 0 = married and 1 = single. A frequency table for this variable has been generated below:

The credit limit is a quantitative response variable as required by this assignment. The values range from 10,000 to 1 million Taiwanese dollars with 168,304 as the mean.

The data that I used in my experiment didn't have any extreme outliers. There was 1 person in the study who had a 1,000,000 Taiwanese dollar credit limit which appears as a slight outlier, however a steady set of credit limits appear in the data starting at 800,000. For this reason I decided not to exclude any outliers. 

The linear regression assumptions of normality of the residuals, linearity between the explanatory variable and response variable, homoscedasticity and independence were out of scope for this assignment.

Experiment Results

The results of the linear regression model indicated single customers in our study (Beta = -25,787.23, p = 2.46e-65) had a negative association with the amount of credit allocated vs. married couples. In other words, married customers in our study had more credit allocated than single customers. Since the p-value is less than our alpha of 0.05 we reject the null hypothesis.

If we take the general regression equation: Y = intercept + slope(X) and adjust it for this experiment: 

Credit limit = intercept + slope(Marital status) and enter the results above for a customer with a 'single' marital status (X = 1):

Credit limit = 182,200 + -25,787.23 (1) = 156,412.77 which is also the average credit limit for a single customer as shown above (156,413).

Further, if we take the mean credit limit of a single customer (156,413) and subtract the mean credit limit of a married customer (182,200) we get -25,787 with is the slope shown above.

Appendix - Python code

import pandas as pandas import statsmodels.formula.api as smf import matplotlib.pyplot as plt import seaborn

#Read my local data into a data frame (DF). defaultRateDF = pandas.read_csv('default_rate_credit_card_clients.csv', low_memory=False)

# Set numerical display formats. pandas.options.display.float_format = '{:0,.2f}'.format

# convert variables I'm using to numeric format. defaultRateDF['X4_MARRIAGE'] = pandas.to_numeric(defaultRateDF['X4_MARRIAGE'], errors='coerce') defaultRateDF['X1_LIMIT_BAL'] = pandas.to_numeric(defaultRateDF['X1_LIMIT_BAL'], errors='coerce')

# Create a new dataframe and add the columns that I intend to use. marriageCreditDF = defaultRateDF[['X4_MARRIAGE', 'X1_LIMIT_BAL']]

# Rename the column names to be more readable. marriageCreditDF = marriageCreditDF.rename(columns={'X4_MARRIAGE': 'MarriedSingle', 'X1_LIMIT_BAL': 'CreditLimit'})

# Filter the data for marriage = 1 or 2 ONLY (married or single) marriageCreditDF = marriageCreditDF[(marriageCreditDF['MarriedSingle'] == 1) | (marriageCreditDF['MarriedSingle'] == 2)]

# Recode the marital status variable to: 0 = married and 1 = single. recode1 = {1: 0, 2: 1} marriageCreditDF['MarriedSingle'] = marriageCreditDF['MarriedSingle'].map(recode1)

#print ("DEBUG: marriageCreditDF\n", marriageCreditDF.head(20))

# Optional: Use .describe() on credit limit to show the range and prove it's quantitative. print("\nQuantitative credit limit response variable:\n", marriageCreditDF['CreditLimit'].describe())

# Post a frequency table for (recoded) categorical explanatory # variable; marital status. marriageCreditCT = pandas.crosstab(marriageCreditDF['MarriedSingle'], marriageCreditDF['CreditLimit'], margins=True) print ("\nFrequency table for recoded categorical explanatory variable (Married = 0, Single = 1):\n\n", marriageCreditCT)

# --------------------------------- # Linear Regression

print ("\nOLS regression model for the association between marital status and credit limit") reg1 = smf.ols('CreditLimit ~ C(MarriedSingle)', data=marriageCreditDF).fit() print (reg1.summary())

print ("\nEasier to read regression results:\n", reg1.params)

# group means and standard deviation. print ("\nMean credit limit for married(0) and single(1) customers") ds1 = marriageCreditDF.groupby('MarriedSingle').mean() print (ds1)

print ("\nCredit limit standard deviation for married(0) and single(1) customers") ds2 = marriageCreditDF.groupby('MarriedSingle').std() print (ds2)

# bivariate bar graph barChart = seaborn.factorplot(x="MarriedSingle", y="CreditLimit", data=marriageCreditDF, kind="bar", ci=None) plt.title("Mean credit for married vs. single customers\n") plt.xlabel('Married                      Single') plt.ylabel('Credit Limit') plt.show(barChart)

# Box and whisker plot. boxplot1 = seaborn.boxplot(x="MarriedSingle", y="CreditLimit", data=marriageCreditDF) plt.title("Dispersion of credit allocated to married vs. single customers\n") plt.xlabel('Married                                                              Single') plt.ylabel('Credit Limit') plt.show(boxplot1)

[mlfanatic]

Course: Regression Modeling in Practice - Week 1 assignment - Writing about your data

This weeks assignment was to write about a data set that I've chosen, describe how the data was collected and the data management steps.

Sample

The sample is the "default of credit card clients data set". It is available at the University of California, Irvine Machine Learning Repository.

The credit card data contains individual customers in Taiwan (N=30,000). The data contains men (39.6%) and women (60.3%). Various education levels such as high school (16%), university (46.7%), graduate school (35%) with the remainder in other educational categories (1.5%). Marital status includes married (45.5%), single (53%) and other (1%). The customers are aged 21 to 79 with 34 being the median age.

Procedure

Data were collected randomly from a bank database, however there was an attempt to balance the data set by including more payments that were in default.

The original purpose of collecting this data was to compare the "predictive accuracy of probability of default among six data mining methods". The original study mentioned that the estimated probability of default was more valuable since the real probability of default is unknown. Also, the probability of default was more valuable than binary classification that might determine whether a client should be granted credit or not.

The data contains monthly credit card payment records from April 2005 through September 2005.

Measures

The data was collected by the Department of Information Management, Chung Hua University, Taiwan and the Department of Civil Engineering, Tamkang University, Taiwan.  

All monetary information is presented in Taiwanese dollars (NT$). 

The measures are described below.

Unique ID identifies the credit card payment history for the credit card holder but yet doesn't specifically identify the credit card holder. 

The explanatory variables are:

Amount of credit provided to the customer 

Gender (male, female)

Education level (high school, university, graduate, other) 

Marital status (married, single, other) 

Credit card customers age in years 

History of past credit card payments (unique states show whether the payment was made in full or if there were payment delays for each month there was a delay).

Credit card bill amount in Taiwanese dollars for each month during the time period mentioned above. 

Credit card payments in Taiwanese dollars for each month during the time period mentioned above.

The response variable is:

a binary variable to denote whether a credit card payment default occurred at any point during the time period mentioned above.

The sample data will be managed in a local file store. Data will be binned as appropriate for future experiments as defined in this course.

Overall, the data is a large, thorough sample with enough variables that can be used in future my experiments. Also, the data doesn't have any missing values that will require me to reduce the data set or use other methods to fill-in the missing data.

[mlfanatic]

Data Analysis Tools course - Week 4 assignment - Testing a Potential Moderator

This post provides the details to test a potential moderator. A moderator is a third variable that affects the direction and or strength of the relation between an explanatory variable and a response variable. Another term for moderator is tertium quid. The Free Dictionary defines tertium quid as "an unknown or indefinite thing related in some way to two known or definite things, but distinct from both."

I chose to use ANOVA in the experiments below. In summary, ANOVA takes a categorical explanatory variable and a quantitative response variable. The moderator can be categorical.

Below are 2 experiments and results that are titled "overall ANOVA" and "ANOVA with moderation". Also, the Python code is in the appendix below.

Overall ANOVA experiment

Null hypothesis: there is no relationship between the age at onset of cannabis dependence and the number of subsequent cannabis dependence episodes.

Alternate hypothesis: there is a relationship between the age at onset of cannabis dependence and the number of subsequent cannabis dependence episodes.

This experiment uses the NESARC data set. The chosen variables are below:

- Categorical explanatory variable:  S3CD5Q15B - AGE AT ONSET OF CANNABIS DEPENDENCE

For this variable the median age is 18. I chose the median age to create 2 age categories since the S3CD5Q15C variable below contains scores of mostly 1 or 2 episodes. Age seemed like a fair way to split the data in half. The plot in the appendix below provides more details.

- Quantitative Response variable: S3CD5Q15C - NUMBER OF EPISODES OF CANNABIS DEPENDENCE

- Moderators: MAJORDEPLIFE - MAJOR DEPRESSION LIFETIME (NON-HIERARCHICAL). The values can be either No or Yes.

Overall ANOVA results

F-statistic: 0.2117

Prob (F-statistic): 0.646

These results are not significant, thus we will not reject the null hypothesis. If we did reject the null hypothesis we would be wrong ~64% of the time.

In addition, the information and image below show that the mean number of dependence episodes for the 2 age categories aren't significantly different. Also, given that the number of dependence episodes range from 1 to 50 the standard deviation around the mean isn't very wide which shows the mean is a good representative model for this variable.

ANOVA with moderation experiment

Since we didn’t reject the null hypothesis above we typically wouldn't run an additional experiment to get more details. However, since this is an assignment I needed to run another experiment with a moderator. The moderator is the NESARC variable MAJORDEPLIFE as described above. I hypothesize that major lifetime depression may encourage the use of cannabis earlier in life and/or increase the number of cannabis dependence episodes.This particular study summarizes a number of other studies that describe this potential relationship.

Null hypothesis: the presence of major lifetime depression does not affect the age of cannabis dependence or number of cannabis dependence episodes.

Alternate Hypothesis: the presence of major lifetime depression does affect the age of cannabis dependence or number of cannabis dependence episodes.

ANOVA with moderation results

Moderation #1

The association between the cannabis dependency age category and number of episodes of dependence for those with NO major lifetime depression is not significant. The F-statistic is 0.1318 and Prob (F-statistic) is 0.717.

In addition, the chart below clearly shows that there isn't a linear relationship between those with NO major lifetime depression, the age at onset of cannabis dependence and number of dependence episodes.

Note, dots on the charts below represent the frequency of occurrence; frequent occurrences are darker (the row of 1′s on the Y-axis) vs. lighter dots that have less frequent occurrences (the single dot on 50 on the Y-axis).

Moderation #2

The association between the cannabis dependency age category and number of episodes of dependence for those WITH major lifetime depression is not significant. The F-statistic 0.06004 and Prob (F-statistic) 0.807.

In addition, the chart below shows that there isn't much of a linear relationship between those WITH major lifetime depression, the age at onset of cannabis dependence and number of dependence episodes.

Appendix:

Python code

=====

import numpy import pandas import statsmodels.formula.api as smf import seaborn import matplotlib.pyplot as plt

# Removes the runtime warning... numpy.seterr(divide='ignore', invalid='ignore')

nesarcDF = pandas.read_csv('nesarc_pds_subset.csv', low_memory=False)

# This bug fix was provided for display formats to avoid run time errors # Put after code for loading data above. pandas.set_option('display.float_format', lambda x:'%.2f'%x)

# Categorical explanatory: S3CD5Q15B AGE AT ONSET OF CANNABIS DEPENDENCE # Quantitative Response variable = S3CD5Q15C NUMBER OF EPISODES OF CANNABIS DEPENDENCE # Moderators: MAJORDEPLIFE MAJOR DEPRESSION - LIFETIME (NON-HIERARCHICAL)

# Create a new dataframe and add the 3 columns that I intend to use for ANOVA. anovaDF = nesarcDF[['S3CD5Q15B', 'S3CD5Q15C', 'MAJORDEPLIFE']]

# Rename 2 of the column names to be more readable. anovaDF = anovaDF.rename(columns={'S3CD5Q15B': 'AgeOfDependence', 'S3CD5Q15C': 'NumEpisodesDependence'})

# convert to numeric. anovaDF['AgeOfDependence'] = pandas.to_numeric(anovaDF['AgeOfDependence'], errors='coerce') anovaDF['NumEpisodesDependence'] = pandas.to_numeric(anovaDF['NumEpisodesDependence'], errors='coerce') anovaDF['MAJORDEPLIFE'] = pandas.to_numeric(anovaDF['MAJORDEPLIFE'], errors='coerce')

# Handle missing data and spaces. 99 are unknown in the variables below. # NumEpisodesDependence = 98 is an outlier with only 4 responses so I removed it. # Lastly drop NAN's. anovaDF['AgeOfDependence'] = anovaDF['AgeOfDependence'].replace(99, numpy.nan) anovaDF['NumEpisodesDependence'] = anovaDF['NumEpisodesDependence'].replace(99, numpy.nan) anovaDF['NumEpisodesDependence'] = anovaDF['NumEpisodesDependence'].replace(98, numpy.nan) anovaDF['MAJORDEPLIFE'] = anovaDF['MAJORDEPLIFE'].replace(99, numpy.nan) anovaDF = anovaDF.replace(' ',numpy.nan) anovaDF = anovaDF.dropna()

# Show why I chose to categorize the data by age print("\nThe following shows why I chose a median age of 18 for 2 categorical variables.") bellc= seaborn.distplot(anovaDF['AgeOfDependence']) plt.xlabel('Age at onset of cannabis dependence') plt.show(bellc)

# Create 2 categories from the ages. def ageOfDependency (row):                                                                                  if row['AgeOfDependence'] <= 18:     return "18 or younger"  if row['AgeOfDependence'] > 18:     return "19 or older"

anovaDF['AgeOfDependencyCategory'] = anovaDF.apply (lambda row: ageOfDependency (row),axis=1)  

# Debug #print(anovaDF.head(10))

print ("\n===== Run #1 =====") print ("===== Overall ANOVA results =====\n") model1 = smf.ols(formula='NumEpisodesDependence ~ C(AgeOfDependencyCategory)', data=anovaDF).fit() print (model1.summary())

subset1 = anovaDF[['NumEpisodesDependence', 'AgeOfDependencyCategory']].dropna()

#print ("\nMeans for number of episodes of cannabis dependence by categorized ages (A vs. B)") print ("\nCannabis dependence episode means by age of dependency categories (A vs. B)") mean1= subset1.groupby('AgeOfDependencyCategory').mean() print (mean1)

#print ("\nStandard deviation for number of episodes of cannabis dependence by categorized ages (A vs. B)") print ("\nCannabis dependence episode standard deviations by age of dependency categories (A vs. B)") stdDev1= subset1.groupby('AgeOfDependencyCategory').std() print (stdDev1)

# bivariate bar graph seaborn.factorplot(x="AgeOfDependencyCategory", y="NumEpisodesDependence", data=anovaDF, kind="bar", ci=None, size=5) plt.xlabel('Cannabis dependency age category') plt.ylabel('Mean episodes of dependence') plt.show()

print ("\n===== Run #2 =====") print ("\n===== ANOVA with moderation results (a third variable) =====\n") #======================================================================= # Testing moderation in the context of ANOVA # Run two separate ANOVAS, one for each level of a third/moderator variable - major lifetime depression. subset2=anovaDF[(anovaDF['MAJORDEPLIFE']==0)] subset3=anovaDF[(anovaDF['MAJORDEPLIFE']==1)]

print("Moderation #1\n") print ('\nAssociation between cannabis dependency age category and number of episodes of dependence for those with NO major lifetime depression') model2 = smf.ols(formula='NumEpisodesDependence ~ C(AgeOfDependencyCategory)', data=subset2).fit() print (model2.summary())

print("Moderation #2\n") print ('\nAssociation between cannabis dependency age category and number of episodes of dependence for those WITH major lifetime depression') model3 = smf.ols(formula='NumEpisodesDependence ~ C(AgeOfDependencyCategory)', data=subset3).fit() print (model3.summary())

#print ("\nMeans for number of episodes of dependence by cannabis dependency age category (A vs. B) for those with NO major lifetime depression") #print ("\nAge of dependence mean for each number of cannabis dependence episodes (A vs. B) for those with NO major lifetime depression") print ("\nFor those with NO major lifetime depression - Each cannabis dependence episode has a mean age of dependence (A vs. B)") mean2= subset2.groupby('NumEpisodesDependence').mean() print (mean2)

print ("\nFor those WITH major lifetime depression - Each cannabis dependence episode has a mean age of dependence (A vs. B)") mean3 = subset3.groupby('NumEpisodesDependence').mean() print (mean3)

# Note: I left these other charts in since they are cool, but not what I wanted. #seaborn.factorplot(x="AgeOfDependence", y="NumEpisodesDependence", data=subset2, kind="point", ci=None, size=8) #seaborn.swarmplot(x="AgeOfDependence", y="NumEpisodesDependence", data=subset2, size=9) seaborn.lmplot(x="AgeOfDependence", y="NumEpisodesDependence", data=subset2, size=7) plt.xlabel('Age at onset of cannabis dependence') plt.ylabel('Number of cannabis depedence episodes') plt.title('For those with NO major lifetime depression show the association\n between the age at onset of cannabis dependence and number of dependence episodes\n') plt.show()

# Note: I left these other charts in since they are cool, but not what I wanted. #seaborn.factorplot(x="AgeOfDependence", y="NumEpisodesDependence", data=subset3, kind="point", ci=None, size=8) #seaborn.swarmplot(x="AgeOfDependence", y="NumEpisodesDependence", data=subset3, size=8) seaborn.lmplot(x="AgeOfDependence", y="NumEpisodesDependence", data=subset3, ci=None, size=7) plt.xlabel('Age at onset of cannabis dependence') plt.ylabel('Number of cannabis depedence episodes') plt.title('For those WITH major lifetime depression show the association\n between the age at onset of cannabis dependence and number of dependence episodes\n') plt.show()

=====

The chart below shows why I chose to use the median age of 18 to create 2 age categories for the variable S3CD5Q15B “Age at onset of cannabis dependence”. This seemed like a fair way to split the data in half for 2 age categories.

[mlfanatic]

Data Analysis Tools course - Week 3 assignment - Generating a Correlation Coefficient

The Pearson correlation coefficient measures the linear relationship between 2 quantitative variables. It seeks to show the strength of a relationship, but doesn't imply causation, between the 2 variables.

In summary, Pearson's correlation is a parametric test and thus has 4 assumptions that must be met for the results to be valid. Those assumptions are; normally distributed data, homogeneity of variance, interval data and independent data. [Source: Discovering Statistics using R, page 167-168]. It is beyond the scope of this assignment to determine if the data used below met these 4 assumptions.

More about Pearson's correlation coefficient can be found here.

 The experiment

This week I am continuing to use the NESARC data set, although the variables had to change to be quantitative.

Null Hypothesis: there is no relationship between the age when cannabis was first used and the number of episodes of cannabis dependence.

Alternate Hypothesis: there is a relationship between the age when cannabis was first used and the number of episodes of cannabis dependence.

I chose 0.05 as the alpha. Alpha is the probability of a Type 1 error or rejecting the null hypothesis when it should not have been rejected.

Specifically, my explanatory variable from the NESARC data is S3BD5Q2A AGE FIRST USED CANNABIS. My response variable is S3CD5Q15C NUMBER OF EPISODES OF CANNABIS DEPENDENCE (1-98 separate periods). 

Python Code

=====

import pandas

import numpy

import seaborn

import matplotlib.pyplot as plt

from scipy.stats import pearsonr

# ---------------------------------------------------

def RunPearsonCorrelation(dataframe1Col, dataframe2Col, chosenAlpha):

    r_value, p_value = pearsonr(dataframe1Col, dataframe2Col)

    coef_determination = r_value**2

    coef_determination_calc = 1 - coef_determination

    coef_determination_percent = coef_determination_calc * 100

    print("\nPearson correlation coefficient: %.2f \nCoefficient of determination (r squared): %.3f \np-value: %.2f \n" % (r_value, coef_determination, p_value))

    if p_value <= chosenAlpha:

        print ("Reject the null hypothesis since the p-value %.2f is less than alpha %s." % (p_value, chosenAlpha))

    else:

        print ("Do not reject the null hypothesis since the p-value %.2f is greater than alpha %s." % (p_value, chosenAlpha))

    print ("\n%.2f percent of the variability between the variables analyzed is unaccounted for." % (coef_determination_percent))

# ---------------------------------------------------

# Removes the runtime warning...

numpy.seterr(divide='ignore', invalid='ignore')

#Read my local NESARC data into a data frame (DF).

nesarcDF = pandas.read_csv('nesarc_pds_subset.csv', low_memory=False)

# This bug fix was provided for display formats to avoid run time errors

# Put after code for loading data above.

pandas.set_option('display.float_format', lambda x:'%f'%x)

# Set variables to numeric.

# S3BD5Q2A AGE FIRST USED CANNABIS

# S3CD5Q15C NUMBER OF EPISODES OF CANNABIS DEPENDENCE

nesarcDF['S3BD5Q2A'] = pandas.to_numeric(nesarcDF['S3BD5Q2A'], errors='coerce')

nesarcDF['S3CD5Q15C'] = pandas.to_numeric(nesarcDF['S3CD5Q15C'], errors='coerce')

# Handle missing data. 99 are unknown in the variables below.

nesarcDF['S3BD5Q2A'] = nesarcDF['S3BD5Q2A'].replace(99, numpy.nan)

nesarcDF['S3CD5Q15C'] = nesarcDF['S3CD5Q15C'].replace(99, numpy.nan)

# Create a new dataframe and add the 2 columns that I intend to use for correlation.

# Further prepare these columns.

# Must drop NAN's otherwise pearsonr() will give an error.

ageDF = nesarcDF[['S3BD5Q2A', 'S3CD5Q15C']]

ageDF = ageDF.replace(' ',numpy.nan)

ageDF = ageDF.dropna(subset = ['S3BD5Q2A', 'S3CD5Q15C'])

print ("Age DF\n",ageDF )

# Pearson correlation and coefficient of determination

print("\n===== Pearson correlation and Coefficient of determination =====")

# Create a constant for the alpha value.

ALPHA = 0.05

print ("\nComparison #1")

print ("--------------")

print("\nVariables analyzed: Age first used cannabis compared to number of episodes of cannabis dependence.")

RunPearsonCorrelation(ageDF['S3BD5Q2A'], ageDF['S3CD5Q15C'], ALPHA)

# Create and show a scatter plot.

scat1 = seaborn.regplot(x="S3BD5Q2A", y="S3CD5Q15C", fit_reg=True, data=ageDF)

plt.xlabel('Age first used cannabis')

plt.ylabel('number of episodes of cannabis dependence')

plt.title('Age first used cannabis compared to number of episodes of cannabis dependence')

plt.show()

=====

Pearson correlation and coefficient of determination results

 Interpretation of results

P-value:

In the results above the 0.08 p-value is not below the 0.05 alpha so we do not reject the null hypothesis.

To be clear, this does not mean that the null hypothesis is "true" or that we "accept" the null hypothesis because the variables used may not be inconsistent with the null hypothesis. Based on limited information we are unable to detect a relationship between the 2 variables that is significant enough to reject the null hypothesis [Source: Business Statistics page 272].

 Pearson's correlation coefficient:

The result I received in the comparison above is close to zero (-0.07) which shows that there is no significant linear relationship between the variables used in this analysis.

 Coefficient of determination:

The coefficient of determination is calculated by squaring Pearson's correlation coefficient. The coefficient of determination tells us how well changes in the explanatory variable can predict the response variable. Generally, it says that the explanatory variable shares some portion of the variability with the response variable.

I received very low coefficient of determination results (0.005) which means that changes in the explanatory variables can't be used to accurately predict changes in the response variable.

Alternatively, the coefficient of determination can be used to show how much variability is unaccounted for between variables simply by subtracting 1 from the coefficient of determination and multiplying by 100 to receive a percent. (For example, using the results above: 1 - 0.005 = 0.995 * 100 = 99.5% of the variability between 2 variables is unaccounted for).

[mlfanatic]

Data Analysis Tools course - Week 2 assignment - Running a Chi-Square Test of Independence

A Chi-Square test is a test of independence between 2 categorical variables. This test compares "the frequencies you observe in certain categories to the frequencies you might expect to get in those categories by chance." [Source: Discovering Statistics using R, page 814]

This week I am continuing to use the NESARC data set.

My hypothesis for the Chi-Square test are:

Null Hypothesis: there is no relationship between the two categorical variables cannabis use and dependency over a 12 month period.

Alternate Hypothesis: there is a relationship between the two categorical variables cannabis use and dependency over a 12 month period.

I chose 0.05 as the alpha. Alpha is the probability of a Type 1 error or rejecting the null hypothesis when it is actually true.

Specifically, my explanatory variable from the NESARC data is "S3BD5Q2C" which is how often cannabis was used over a 12 month period. My response variable is "MAR12ABDEP" which is cannabis dependence over a 12 month period.

The variable "MAR12ABDEP" has the following response values:

Frequency of responses, item value and description

---------------------------------------------------

42,533 0. No cannabis diagnosis

427      1. Cannabis abuse only

23        2. Cannabis dependence only

110      3. Cannabis abuse and dependence

 Since we are looking for cannabis dependence I chose to use item values "0 - No cannabis diagnosis" and "2 - Cannabis dependence only".

 Python Code

========================

import pandas import numpy import scipy.stats import seaborn import matplotlib.pyplot as plt

# ----------------------------------------------------------------------- # The following procedure will run the Chi-Square test using the params below. # group 1 and group 2 are the combinations for my UseFreqPerMo variable.

def runTest(group1, group2, subsetDF):

   recodeDict = {}

   if group1 != group2:        print("---------------------------------------------------\n")        print("Tables and Chi-Square Results for explanatory variable groups %s and %s\n" % (group1, group2))

       # Dictionary should look like this: {1:1, 2.5:2.5}        recodeDict = {group2:group2, group1:group1}

       subsetDF['Comparison1v2'] = subsetDF['UseFreqPerMo'].map(recodeDict)

       # Contingency table with observed counts        contTable2 = pandas.crosstab(subsetDF['MAR12ABDEP'], subsetDF['Comparison1v2'])        print("Contingency Table\n", contTable2)

       #prepare column percentages        columnSum = contTable2.sum(axis=0)        columnPct = contTable2/columnSum        print("\nColumn Percentages\n", columnPct)

       # Run the Chi-square and print results.        chi2Result, pvalue, degreesOfFreedom, expectedFreq = scipy.stats.chi2_contingency(contTable2, correction=False)        print("\n Chi-Square test statistic: %s\n Degrees of freedom: %s\n p-value: %s\n" % (chi2Result, degreesOfFreedom, pvalue))

# -----------------------------------------------------------------------

#Read my local NESARC data into a data frame (DF). nesarcDF = pandas.read_csv('nesarc_pds_subset.csv', low_memory=False)

# This bug fix was provided for display formats to avoid run time errors # Put after code for loading data above. pandas.set_option('display.float_format', lambda x:'%f'%x)

# Set variables to numeric. nesarcDF['S3BD5Q2C'] = pandas.to_numeric(nesarcDF['S3BD5Q2C'], errors='coerce') nesarcDF['S3BQ3'] = pandas.to_numeric(nesarcDF['S3BQ3'], errors='coerce') nesarcDF['MAR12ABDEP'] = pandas.to_numeric(nesarcDF['MAR12ABDEP'], errors='coerce')

# Get a subset of data for young adults age 18 to 25. # These variables include the last month only, hence 'S3BD5Q2C' < 7. # Also, since we are looking for cannabis abuse and dependence I chose to  #use MAR12ABDEP variable with responses 0 "No cannabis diagnosis" and #response 2 Cannabis dependence only".

youngAdultDF = nesarcDF[(nesarcDF['AGE']>=18) & (nesarcDF['AGE']<=25) & (nesarcDF['S3BD5Q2C'] < 7) & ((nesarcDF['MAR12ABDEP']==0 ) | (nesarcDF['MAR12ABDEP'] == 2))]

# Handle missing data. 9 and 99 are unknown in the variables below. youngAdultDF['S3BD5Q2C'] = youngAdultDF['S3BD5Q2C'].replace(99, numpy.nan) youngAdultDF['S3BQ3'] = youngAdultDF['S3BQ3'].replace(99, numpy.nan) youngAdultDF['MAR12ABDEP'] = youngAdultDF['MAR12ABDEP'].replace(99, numpy.nan)

# Recode number of days smoked in the past month into a new variable Use Frequency Per Month. # 1: 30 days smoked, 2: 22 days per mo or 5-6 days per week smoked, etc. recode1 = {1: 30, 2: 22, 3: 14, 4: 6, 5: 2.5, 6: 1} youngAdultDF['UseFreqPerMo']= youngAdultDF['S3BD5Q2C'].map(recode1)

# Convert the new variable UseFreqPerMo to numeric. youngAdultDF['UseFreqPerMo'] = pandas.to_numeric(youngAdultDF['UseFreqPerMo'], errors='coerce')

# Create a secondary variable to hold the Number of joints smoked per month. # Multiply the days smoked per month by the number of joints smoked per day. youngAdultDF['NumJointSmokedPerMo_est'] = youngAdultDF['UseFreqPerMo'] * youngAdultDF['S3BQ3'] youngAdultDF['NumJointSmokedPerMo_est'] = pandas.to_numeric(youngAdultDF['NumJointSmokedPerMo_est'], errors='coerce')

# Run initial Chi-Square test

print("\n===== Initial Chi-Square Test results =====\n")

subset4 = youngAdultDF[['UseFreqPerMo', 'MAR12ABDEP']].dropna()

contTable = pandas.crosstab(subset4['MAR12ABDEP'], subset4['UseFreqPerMo'])

print("\nContingency Table\n", contTable)

# Create column percentages columnSum = contTable.sum(axis=0) columnPct = contTable/columnSum print("\nColumn Percentages\n", columnPct)

# run the Chi-square print("\nInitial Chi-Square Results") print("---------------------")

chi2Result, pvalue, degreesOfFreedom, expectedFreq = scipy.stats.chi2_contingency(contTable, correction=False) print(" Chi-Square test statistic: %s\n Degrees of freedom: %s\n p-value: %s\n" % (chi2Result, degreesOfFreedom, pvalue))

subset4['UseFreqPerMo'] = subset4['UseFreqPerMo'].astype('category') subset4['MAR12ABDEP'] = pandas.to_numeric(youngAdultDF['MAR12ABDEP'], errors='coerce')

seaborn.factorplot(x="UseFreqPerMo", y="MAR12ABDEP", data=subset4, kind="bar", ci=None) plt.xlabel("Days smoked Cannabis per month") plt.ylabel("Proportion of Cannabis dependence")

# Run post-hoc tests for each explanatory variable group. print("\n===== Post-hoc Chi-Square Tests =====\n")

# I know, nested for loops would work, however this works too...

# Invoke the runTest procedure defined above for each group.

runTest(1, 2.5, subset4) runTest(1, 6, subset4) runTest(1, 14, subset4) runTest(1, 22, subset4) runTest(1, 30, subset4)

runTest(2.5, 6, subset4) runTest(2.5, 14, subset4) runTest(2.5, 22, subset4) runTest(2.5, 30, subset4)

runTest(6, 14, subset4) runTest(6, 22, subset4) runTest(6, 30, subset4)

runTest(14, 22, subset4) runTest(14, 30, subset4)

runTest(22, 30, subset4)

========================

Note: Some interesting information on scipy.stats.chi2_contingency() can be found here.

For example, in the code above I chose to call this function with the 4 return types so that I could better format the output:

chi2Result, pvalue, degreesOfFreedom, expectedFreq = scipy.stats.chi2_contingency(contTable2, correction=False)

print(" Chi-Square test statistic: %s\n p-value: %s\n Degrees of freedom: %s\n" % (chi2Result, pvalue, degreesOfFreedom))

  Initial Chi-Square Test Results

Chi-Square test statistic: 3.41

Degrees of freedom: 5

p-value: 0.64

The Chi-Square test statistic above and the degrees of freedom are used together to determine the p-value, or probability, that the two categorical variables are independent. Degrees of freedom are the number of observations that are free to vary.

The degrees of freedom are calculated using the contingency table by taking the number of rows minus one times the number of columns minus one or (R-1)*(C-1).

In the initial contingency table below we have 2 rows and 6 columns, thus the degrees of freedom are (2-1) * (6-1) = 5.

The p-value is the probability of observing the Chi-Square test statistic. Essentially, there is a 64% probability that any deviation from expected results is due to chance. Since the 0.64 p-value is well above our 0.05 alpha I won't reject the null hypothesis. Based on the NESARC data there is no statistically significant association between cannabis use and cannabis dependency over a 12 month period.

 Note: On the chart above the value for 14.0 "Days cannabis smoked per month" looks odd since the bar doesn't appear, however this is accurate. In the NESARC data the variable "MAR12ABDEP" under the previously chosen option 2 "Cannabis dependence only" has just 23 responses. These responses are eliminated when we filter S3BD5Q2C "HOW OFTEN USED CANNABIS IN THE LAST 12 MONTHS" by "3" or per the "...map(recode1)” in the Python code above which sets this to 14 days.

Overall, there are only 23 cannabis dependence responses, well below 30 which is typically used for the central limit theorem and could affect the results.

Chi-Squared post-hoc paired comparison tests

My explanatory variable, cannabis use, has more than two groups so I conducted Chi-Squared post-hoc tests to determine which groups are different from the others. Do to these tests I used the Bonferroni Adjustment to protect against Type I error. Next, I ran the Chi Square Test of Independence for each paired group.

The Bonferroni Adjustment was calculated as follows: My chosen alpha was 0.05. There were 15 comparisons so my Bonferroni Adjustment is: 0.05/15 = 0.003.The Bonferroni Adjustment evenly divides the alpha among each pair to be tested.

I ran the post-hoc paired tests. The results can be found in the table below.

Next, I looked at the p-values from each paired post-hoc test to determine if it was < 0.003. If it was then I can reject the null hypothesis for that particular pair. However, per the table above there are no p-values that meet this criteria thus we do not reject the null hypothesis for any of the pairs above.

[mlfanatic]

Coursera - Data Analysis tools - Week 1 assignment - Running an analysis of variance (ANOVA)

In the previous course "Data Management and Visualization" I was using the GapMinder data. Since this is the 2nd class in Coursera's Data Analysis and Interpretation Specialization I decided to use the NESARC data set.

Details

I created the following hypothesis for the part 1 test below:

Null hypothesis: There is no relationship between young adults (age 18-25) with major depression and the average number of cannabis joints smoked per month in the last 12 months. 

Alternate Hypothesis: There is a positive relationship between the average number of cannabis joints smoked per month among young adults (age 18-25) and major depression over the last 12 months.

I used ANOVA since I have a categorical variable -> quantitative variable relationship. Major depression is my categorical variable since the NESARC variable MAJORDEP12 is either a 0 or 1. The number of cannabis joints smoked per month is my quantitative variable since it shows a frequency.

All of the NESARC variables I used are below:

S3BD5Q2C -- HOW OFTEN USED CANNABIS IN THE LAST 12 MONTHS

S3BQ3 -- NUMBER OF CANNABIS JOINTS USUALLY SMOKED IN A DAY IN THE LAST 12 MONTHS

MAJORDEP12 -- MAJOR DEPRESSION IN LAST 12 MONTHS (NON-HIERARCHICAL)

ETHRACE2A -- IMPUTED RACE/ETHNICITY [the values are 1 through 5. Also, I chose to re-use this variable from the sample code since I really couldn’t easily find another categorical variable in the NESARC data set.]

AGE -- Age of respondent  

My code

=======================

import numpy import pandas import statsmodels.formula.api as smf import statsmodels.stats.multicomp as multi

#Read my local subset of the NESARC data into a data frame (DF). nesarcDF = pandas.read_csv('nesarc_pds.csv', low_memory=False)

# This bug fix was provided for display formats to avoid run time errors # Put after code for loading data above. pandas.set_option('display.float_format', lambda x:'%f'%x)

# Set variables to numeric. nesarcDF['S3BD5Q2C'] = pandas.to_numeric(nesarcDF['S3BD5Q2C'], errors='coerce') nesarcDF['S3BQ3'] = pandas.to_numeric(nesarcDF['S3BQ3'], errors='coerce')

# Get a subset of data for young adults age 18 to 25. # These variables include the last month only, hence 'S3BD5Q2C' < 7. youngAdultDF=nesarcDF[(nesarcDF['AGE']>=18) & (nesarcDF['AGE']<=25) & (nesarcDF['S3BD5Q2C'] < 7)]

# Handle missing data. 9 and 99 are unknown in the variables below. youngAdultDF['S3BD5Q2C']=youngAdultDF['S3BD5Q2C'].replace(99, numpy.nan) youngAdultDF['S3BQ3']=youngAdultDF['S3BQ3'].replace(99, numpy.nan)

# Recode number of days smoked in the past month into a new variable Use Frequency Per Month. # 1: 30 days smoked, 2: 22 5-6 days smoked, etc. recode1 = {1: 30, 2: 22, 3: 14, 4: 5, 5: 2.5, 6: 1} youngAdultDF['UseFreqPerMo']= youngAdultDF['S3BD5Q2C'].map(recode1)

# Convert the new variable UseFreqPerMo to numeric. youngAdultDF['UseFreqPerMo']= pandas.to_numeric(youngAdultDF['UseFreqPerMo'], errors='coerce')

# Create a secondary variable to hold the Number of Joints Smoked Per Month. # Multiply the days smoked per month by the number of joints smoked per day. youngAdultDF['NumJointSmokedPerMo_est'] = youngAdultDF['UseFreqPerMo'] * youngAdultDF['S3BQ3'] youngAdultDF['NumJointSmokedPerMo_est'] = pandas.to_numeric(youngAdultDF['NumJointSmokedPerMo_est'], errors='coerce')

ct1 = youngAdultDF.groupby('NumJointSmokedPerMo_est').size() print ("\nNumber of joints smoked per month (est) | The number of users who smoked that amount:\n", ct1)

# Use the Ordinary Least Squares function and calculate the F-statistic and associated p value. # Use capital C to tell Python that the following param is a categorical variable. model1 = smf.ols(formula='NumJointSmokedPerMo_est ~ C(MAJORDEP12)', data=youngAdultDF) results1 = model1.fit() print (results1.summary())

subset2 = youngAdultDF[['NumJointSmokedPerMo_est', 'MAJORDEP12']].dropna()

print ('\nMeans for NumJointSmokedPerMo_est by major depression status\n') mean1 = subset2.groupby('MAJORDEP12').mean() print (mean1)

print ('\nStandard deviations for NumJointSmokedPerMo_est by major depression status\n') stdDev1 = subset2.groupby('MAJORDEP12').std() print (stdDev1)

print("===================== Cannabis use and Ethnicity ==========================")

# subset3 represents smoking cannabis per month when compared to ethnicity. subset3 = youngAdultDF[['NumJointSmokedPerMo_est', 'ETHRACE2A']].dropna()

# Use the Ordinary Least Squares function and calculate the F-statistic and associated p value. model2 = smf.ols(formula='NumJointSmokedPerMo_est ~ C(ETHRACE2A)', data=subset3).fit() print (model2.summary())

print ('\nMeans for NumJointSmokedPerMo_est by ethnicity\n') mean2= subset3.groupby('ETHRACE2A').mean() print (mean2)

print ('\nStandard deviations for NumJointSmokedPerMo_est by ethnicity\n') stdDev2 = subset3.groupby('ETHRACE2A').std() print (stdDev2)

# Post hoc ANOVA results # NumJointSmokedPerMo_est is my quantitative response variable # ETHRACE2A is my categorical explanatory variable. multiComp1 = multi.MultiComparison(subset3['NumJointSmokedPerMo_est'], subset3['ETHRACE2A']) result1 = multiComp1.tukeyhsd() print(result1.summary())

=======================

Note: If you run the code you'll notice that NumJointSmokedPerMo_est has 1 value for 1,920 cannabis joints smoked per month. I triple-checked the data and my calculations. There was one guy in the NESARC data set who said that he usually smoked 64 joints a day in the last 12 months!! Multiply this by 30 days in a month and you get an abnormally high (get it?!) 1,920 joints smoked per month (which is an estimate of course).

 Part 1 - Cannabis use and major depression in the last 12 months results

I used the Ordinary Least Squares function to calculate the F-statistic and associated p-value for NumJointSmokedPerMo_est and MAJORDEP12. My results are:

F-statistic:            0.002492

Prob (F-statistic): 0.960

We get an amazingly high p-value of .960. Since this p-value is not less than our .05 alpha we do not reject the null hypothesis. Based on the NESARC data over the last 12 months there is no relationship between young adults with major depression and the average number of cannabis joints smoked per month. Perhaps the other measures of depression in the NESARC data would have provided slightly different results.

This result is confirmed since the mean number of joints smoked per month (estimate) vs. people with and without major depression in the last 12 months are very close.

                        NumJointSmokedPerMo_est

MAJORDEP12                        

0                         65.660606

1                         66.474490

  Part 2 - Cannabis use and Ethnicity Results

Null Hypothesis: There is no relationship between cannabis use and ethnic group.

Alternate Hypothesis: There is a relationship between cannabis use an ethnic group.

I used the Ordinary Least Squares function to calculate the F-statistic and associated p-value for NumJointSmokedPerMo_est and ETHRACE2A. My results are:

F-statistic:                     2.177

Prob (F-statistic):         0.0708

Since the p-value .0708 is not less than our 0.05 alpha we do not reject the null hypothesis.

Looking at the means for number of joints smoked per month (estimate) by ethnicity, groups 2 and 3 seems to have the highest usage, however the post hoc ANOVA results below show that we shouldn’t reject the null hypothesis.

                       Mean NumJointSmokedPerMo_est

ETHRACE2A                        

1                           53.321918 2                           101.711765

3                           106.807692

4                           66.090909 5                           57.440000

The standard deviations for estimated number of joints smoked per month by ethnicity is very large, thus showing that the means above aren’t an accurate representation of the data.

                         Std Dev NumJointSmokedPerMo_est ETHRACE2A                         1                       103.206606 2                       224.858585 3                       202.608977 4                       113.879282 5                       109.772659

Part 3 - post hoc ANOVA results

I ran a post hoc test to determine which ethnic groups are different from the others. I used Tukey HSD for this analysis. From the results below we can see that we should not reject the null hypothesis. Post hoc comparisons of mean number of cannabis joints smoked per day by ethnicity shows no relationship.

[mlfanatic]

Coursera - Data Management and Visualization - Week 4 assignment - Creating graphs

New code snippets can be found below. The code to import the data and convert it can be found in my earlier posts.

Internet User Rate

print ("\nInternet user rate details:") print (gapminderDF['internetuserate'].describe())

print ("\nInternet use rate:") print (gapminderDF['internetuserate'].hist(bins=5))

The chart shows the Internet user rate per 100 people (x-axis) across various countries (y-axis). The data is skewed to the right.

It's interesting to note that after the tallest bar on the left (showing the lowest internet user rates) the rate seems to show less variability. This suggests that countries in this study can be separated into low Internet use (79 countries) and high Internet use (113 countries). The causes for this are beyond the scope of this study.

Income per person

print ("\nIncome per person details:") print (gapminderDF['incomeperperson'].describe())

print ("\nIncome per person:") print (gapminderDF['incomeperperson'].hist(bins=5))

The X-axis is income per person. Y-axis is number of countries. For example, the longest bar shows that there are ~161 countries with an income per person of ~$20,000 or less.

The data is highly skewed to the right with a vast majority of the income per person falling into the first bar: ~$100 to ~$20,000 range.

Polity Score

print ("\nPolity score details:") print (gapminderDF['polityscore'].describe())

print ("\nPolity scores:") print (gapminderDF['polityscore'].hist(bins=5))

The X-axis is the polity score. Y-axis is number of countries. For example, the first bar on the left shows that there are 20 countries with a polity score of -7 or less (very autocratic and less free).

When comparing the tallest bar on the right (showing 80 countries with the highest democratic and free nature) to the remaining bars on the left (81 countries with lower democratic and free) suggests that countries can be clearly separated into 2 groups. The causes for this are beyond the scope of this study.

Income per person and Internet use rate

scat1 = seaborn.regplot(x="incomeperperson", y="internetuserate", data=gapminderDF) plt.xlabel('income per person') plt.ylabel('internet use rate') plt.title('Scatterplot for Income per Person and Internet Use Rate')

This scatter plot shows the Income per person compared to Internet use rate per 100 people. The chart above mostly disproves my original hypothesis.

My original hypothesis stated that there should a positive relationship between Internet use and income. I assumed that wealthier countries would have greater Internet use due to better infrastructure, more disposable income to buy Internet connected devices and pay for services available over the Internet (e.g.: streaming TV shows, etc.).

However, the chart above clearly shows that countries with the highest incomes have roughly the same ratio of citizens who use the internet as those in the ~$10,000 to $40,000 income range. Although, countries with the lowest income (~$10,000 and below) have the lowest Internet usage.

Internet use rate and polity score

scat2 = seaborn.regplot(x="internetuserate", y="polityscore", data=gapminderDF) plt.xlabel('internet use rate') plt.ylabel('polity score') plt.title('Scatterplot for internet use rate and polity score')

The scatter plot above mostly disproves my secondary hypothesis.

My secondary hypothesis stated that there should be a positive relationship between Internet use and a countries polity score. The polity score reflects the democratic and free nature (autocracy/democracy) of each country in the GapMinder data. I assumed there would be a positive relationship between these variables because citizens in more free and democratic countries are allowed to self-organize and form “news” outlets that speak with viewpoints different from government approved news outlets.

The chart above doesn’t show a strong relationship between these variables. We can see that in the column representing internet use rate of 0 to 20 the polity score varies from -9 (autocratic and less free) to 10 (democratic and more free). There are similar observations for Internet use rates from 20 to 60.

[mlfanatic]

Coursera - Data Management and Visualization - Week 3 assignment - Managing Data

This week we focused on various approaches to clean, summarize and group data at a higher level.   My primary research question is whether Internet use is related to income per person. My secondary research question is whether Internet use is related to the countries democratic and free nature (autocracy/democracy) score.   Technical note for Python users, if the input has spaces where you expect numeric values when you convert the input using “to_numeric” with the “errors=coerce” param it will convert the spaces to NaN’s for you. For example the following will take care of spaces in my polityscore column:         pandas.to_numeric(gapminderDF[‘polityscore’], errors='coerce’)   I don’t see any need to add the newest code in the week 3 videos to convert spaces to NaN’s. Also, you can pass the equivalent of a space with a Python regular expression.         gapminderDF['incomeperperson’].replace('  ’, numpy.nan)   Each variable below was grouped by quartiles using the pandas.qcut function and the pandas.cut function. I’ll let Wikipedia explain percentiles, quartiles, etc.   The pandas.cut function provides a different perspective of the data because it separates the variable observations into 4 equal-width bins.    Note, I took out the code from last week and only show the new code added this week. The frequency distributions can be found in my prior post here. Only new information is presented below.

My code

# ————————————————————————— # Polity score frequency distribution # quartile split (use qcut function & ask for 4 groups - gives you quartile split) print ('Polity score distributions’) gapminderDF['polityscorePercentileGroup’]=pandas.qcut(gapminderDF.polityscore, 4, labels=[“0th percentile”,“25th percentile”,“50th percentile”,“75th percentile”]) gapminderQuartiles = gapminderDF['polityscorePercentileGroup’].value_counts(sort=False, dropna=True) print(“\nPolity Score percentiles:\n”,gapminderQuartiles)

# categorize quantitative variable based on customized splits using cut function # splits into 4 groups (-10 to -5, -4 to 0, 1 to 5 and 6 to 10) gapminderDF['polityscoreGroups’] = pandas.cut(gapminderDF.polityscore, [-11, -5, 0, 5, 10]) gapminderGroups = gapminderDF['polityscoreGroups’].value_counts(sort=False, dropna=False) print(“\nPolity Score groups\n”,gapminderGroups)

# ————————————————————————— # Internet use rate distribution # quartile split (use qcut function & ask for 4 groups - gives you quartile split) print (’\nInternet use rate distributions’) gapminderDF['internetuseratePercentileGroup’]=pandas.qcut(gapminderDF.internetuserate, 4, labels=[“0th percentile”,“25th percentile”,“50th percentile”,“75th percentile”]) gapminderQuartiles = gapminderDF['internetuseratePercentileGroup’].value_counts(sort=False, dropna=True) print(“\nInternet use rate percentiles:\n”,gapminderQuartiles)

# categorize quantitative variable based on customized splits using cut function # splits into 4 groups/bins. gapminderDF['internetuserateGroups’] = pandas.cut(gapminderDF.internetuserate, 4) gapminderGroups = gapminderDF['internetuserateGroups’].value_counts(sort=False, dropna=False) print(“\nInternet use rate groups\n”,gapminderGroups)

# ————————————————————————— # Income per person distribution # quartile split (use qcut function & ask for 4 groups - gives you quartile split) print (’\nIncome per person distributions’) gapminderDF['incomeperpersonPercentileGroup’]=pandas.qcut(gapminderDF.incomeperperson, 4, labels=[“0th percentile”,“25th percentile”,“50th percentile”,“75th percentile”]) gapminderQuartiles = gapminderDF['incomeperpersonPercentileGroup’].value_counts(sort=False, dropna=True) print(“\nIncome per person percentiles:\n”,gapminderQuartiles)

# categorize quantitative variable based on customized splits using cut function # splits into 4 group/bins gapminderDF['incomeperpersonGroups’] = pandas.cut(gapminderDF.incomeperperson, 4) gapminderGroups = gapminderDF['incomeperpersonGroups’].value_counts(sort=False, dropna=False) print(“\nIncome per person groups\n”,gapminderGroups)

# ————————————————————————— print (“\nExta information”) print (“—————-”) gapminderstats = gapminderDF.groupby('polityscore’).agg({'internetuserate’: [numpy.size, numpy.mean, numpy.std]}) print(“\nInternet use rate grouped by polity score:\n\n”, gapminderstats)

The polity score was pre-calculated by GapMinder. It is calculated by subtracting an autocracy score from a democracy score and is the summary measure of a country’s democratic and free nature. -10 is the lowest value, 10 the highest.

The polity score percentiles are mostly uniform. However, the polity score groups are unevenly distributed. Less than half of the countries, 71, are ranked 4 and below with the remaining 90 countries ranked 5 and above. This means that a large portion of countries in the GapMinder data are considered democratic and free.

The Internet use rate is Internet users (per 100 people) in 2010 across various countries.

You’ll note that the internet use rate percentiles are exactly 48. I double-checked the data and this is correct.

In the Internet use rate groups the most common range of variables range from between zero and 47. Only about 15% of countries reach an Internet use rate above 70%. More research would need to be done to determine the main causes that impact this variable. I would guess tele-comm infrastructure, wide spread reliable sources of electricity and disposable income might be some main determinants.  

Income per person is 2010 Gross Domestic Product per capita in constant 2000 US dollars.  

The percentiles seem too uniform but they are accurate.

The income per person groups are very uneven with ~89% of the observations in the first group.  

Beyond the assignment…

Beyond the assignment I added a distribution that groups internet use rate by polity score.

In the example circled above there are 2 countries with a -10 polity score. In those 2 countries the mean internet use rate is 61.30 (41 and 81). That mean has a standard deviation of 28.70.  

In this particular example there are too few observations to draw any conclusions. Also, the standard deviation is fairly broad so the mean isn’t a very good representation of the internet use rate.