Attribution Models in Marketing

April 18, 2019/0 Comments/in Big Data, Business Analytics, Business Intelligence, Data Mining, Data Science, Insights, Machine Learning, Main Category, Predictive Analytics, Statistics, Use Case, Use Cases /by Barbara Staron

Attribution Models

A Business and Statistical Case

INTRODUCTION

A desire to understand the causal effect of campaigns on KPIs

Advertising and marketing costs represent a huge and ever more growing part of the budget of companies. Studies have found out this share is as high as 10% and increases with the size of companies (CMO study by American Marketing Association and Duke University, 2017). Measuring precisely the impact of a specific marketing campaign on the sales of a company is a critical step towards an efficient allocation of this budget. Would the return be higher for an euro spent on a Facebook ad, or should we better spend it on a TV spot? How much should I spend on Twitter ads given the volume of sales this channel is responsible for?

Attribution Models have lately received great attention in Marketing departments to answer these issues. The transition from offline to online marketing methods has indeed permitted the collection of multiple individual data throughout the whole customer journey, and allowed for the development of user-centric attribution models. In short, Attribution Models use the information provided by Tracking technologies such as Google Analytics or Webtrekk to understand customer journeys from the first click on a Facebook ad to the final purchase and adequately ponderate the different marketing campaigns encountered depending on their responsibility in the final conversion.

Issues on Causal Effects

A key question then becomes: how to declare a channel is responsible for a purchase? In other words, how can we isolate the causal effect or incremental value of a campaign ?

1. A/B-Tests

One method to estimate the pure impact of a campaign is the design of randomized experiments, wherein a control and treated groups are compared. A/B tests belong to this broad category of randomized methods. Provided the groups are a priori similar in every aspect except for the treatment received, all subsequent differences may be attributed solely to the treatment. This method is typically used in medical studies to assess the effect of a drug to cure a disease.

Main practical issues regarding Randomized Methods are:

Assuring that control and treated groups are really similar before treatment. Uually a random assignment (i.e assuring that on a relevant set of observable variables groups are similar) is realized;
Potential spillover-effects, i.e the possibility that the treatment has an impact on the non-treated group as well (Stable unit treatment Value Assumption, or SUTVA in Rubin’s framework);
The costs of conducting such an experiment, and especially the costs linked to the deliberate assignment of individuals to a group with potentially lower results;
The number of such experiments to design if multiple treatments have to be measured;
Difficulties taking into account the interaction effects between campaigns or the effect of spending levels. Indeed, usually A/B tests are led by cutting off temporarily one campaign entirely and measuring the subsequent impact on KPI’s compared to the situation where this campaign is maintained;
The dynamical reproduction of experiments if we assume that treatment effects may change over time.

In the marketing context, multiple campaigns must be tested in a dynamical way, and treatment effect is likely to be heterogeneous among customers, leading to practical issues in the lauching of A/B tests to approximate the incremental value of all campaigns. However, sites with a lot of traffic and conversions can highly benefit from A/B testing as it provides a scientific and straightforward way to approximate a causal impact. Leading companies such as Uber, Netflix or Airbnb rely on internal tools for A/B testing automation, which allow them to basically test any decision they are about to make.

References:

Books:

A/B testing: the most powerful way to turn clicks into customers. Dan Siroker, Pete Koomen; Wiley, 2013.

Blogs:

https://eng.uber.com/xp

https://medium.com/airbnb-engineering/growing-our-host-community-with-online-marketing-9b2302299324

Study:

https://cmosurvey.org/wp-content/uploads/sites/15/2018/08/The_CMO_Survey-Results_by_Firm_and_Industry_Characteristics-Aug-2018.pdf

2. Attribution models

Attribution Models do not demand to create an experimental setting. They take into account existing data and derive insights from the variability of customer journeys. One key difficulty is then to differentiate correlation and causality in the links observed between the exposition to campaigns and purchases. Indeed, selection effects may bias results as exposure to campaigns is usually dependant on user-characteristics and thus may not be necessarily independant from the customer’s baseline conversion probabilities. For example, customers purchasing from a discount price comparison website may be intrinsically different from customers buying from FB ad and this a priori difference may alone explain post-exposure differences in purchasing bahaviours. This intrinsic weakness must be remembered when interpreting Attribution Models results.

2.1 General Issues

The main issues regarding the implementation of Attribution Models are linked to

Causality and fallacious reasonning, as most models do not take into account the aforementionned selection biases.
Their difficult evaluation. Indeed, in almost all attribution models (except for those based on classification, where the accuracy of the model can be computed), the additionnal value brought by the use of a given attribution models cannot be evaluated using existing historical data. This additionnal value can only be approximated by analysing how the implementation of the conclusions of the attribution model have impacted a given KPI.
Tracking issues, leading to an uncorrect reconstruction of customer journeys
- Cross-device journeys: cross-device issue arises from the use of different devices throughout the customer journeys, making it difficult to link datapoints. For example, if a customer searches for a product on his computer but later orders it on his mobile, the AM would then mistakenly consider it an order without prior campaign exposure. Though difficult to measure perfectly, the proportion of cross-device orders can approximate 20-30%.
- Cookies destruction makes it difficult to track the customer his the whole journey. Both regulations and consumers’ rising concerns about data privacy issues mitigate the reliability and use of cookies.1 – From 2002 on, the EU has enacted directives concerning privacy regulation and the extended use of cookies for commercial targeting purposes, which have highly impacted marketing strategies, such as the ‘Privacy and Electronic Communications Directive’ (2002/58/EC). A research was conducted and found out that the adoption of this ‘Privacy Directive’ had led to 64% decrease in advertising methods compared to the rest of the world (Goldfarb et Tucker (2011)). The effect was stronger for generalized sites (Yahoo) than for specialized sites.2 – Users have grown more and more conscious of data privacy issues and have adopted protective measures concerning data privacy, such as automatic destruction of cookies after a session is ended, or simply giving away less personnal information (Goldfarb et Tucker (2012) ) .Valuable user information may be lost, though tracking technologies evolution have permitted to maintain tracking by other means. This issue may be particularly important in countries highly concerned with data privacy issues such as Germany.
- Offline/Online bridge: an Attribution Model should take into account all campaigns to draw valuable insights. However, the exposure to offline campaigns (TV, newspapers) are difficult to track at the user level. One idea to tackle this issue would be to estimate the proportion of conversions led by offline campaigns through AB testing and deduce this proportion from the credit assigned to the online campaigns accounted for in the Attribution Model.
- Touch point information available: clicks are easy to follow but irrelevant to take into account the influence of purely visual campaigns such as display ads or video.

2.2 Today’s main practices

Two main families of Attribution Models exist:

Rule-Based Attribution Models, which have been used for in the last decade but from which companies are gradualy switching.

Attribution depends on the individual journeys that have led to a purchase and is solely based on the rank of the campaign in the journey. Some models focus on a single touch points (First Click, Last Click) while others account for multi-touch journeys (Bathtube, Linear). It can be calculated at the customer level and thus doesn’t require large amounts of data points. We can distinguish two sub-groups of rule-based Attribution Models:

One Touch Attribution Models attribute all credit to a single touch point. The First-Click model attributes all credit for a converion to the first touch point of the customer journey; last touch attributes all credit to the last campaign.
Multi-touch Rule-Based Attribution Models incorporate information on the whole customer journey are thus an improvement compared to one touch models. To this family belong Linear model where credit is split equally between all channels, Bathtube model where 40% of credit is given to first and last clicks and the remaining 20% is distributed equally between the middle channels, or time-decay models where credit assigned to a click diminishes as the time between the click and the order increases..

The main advantages of rule-based models is their simplicity and cost effectiveness. The main problems are:

– They are a priori known and can thus lead to optimization strategies from competitors
– They do not take into account aggregate intelligence on customer journeys and actual incremental values.
– They tend to bias (depending on the model chosen) channels that are over-represented at the beggining or end of the funnel, according to theoretical assumptions that have no observationnal back-ups.

Data-Driven Attribution Models

These models take into account the weaknesses of rule-based models and make a relevant use of available data. Being data-driven, following attribution models cannot be computed using single user level data. On the contrary values are calculated through data aggregation and thus require a certain volume of customer journey information.

References:

https://dspace.mit.edu/handle/1721.1/64920

3. Data-Driven Attribution Models in practice

3.1 Issues

Several issues arise in the computation of campaigns individual impact on a given KPI within a data-driven model.

Selection biases: Exposure to certain types of advertisement is usually highly correlated to non-observable variables which are in turn correlated to consumption practices. Differences in the behaviour of users exposed to different campaigns may thus only be driven by core differences in conversion probabilities between groups whether than by the campaign effect.
Complementarity: it may be that campaigns A and B only have an effect when combined, so that measuring their individual impact would lead to misleading conclusions. The model could then try to assess the effect of combinations of campaigns on top of the effect of individual campaigns. As the number of possible non-ordered combinations of k campaigns is 2k, it becomes clear that inclusing all possible combinations would however be time-consuming.
Order-sensitivity: The effect of a campaign A may depend on the place where it appears in the customer journey, meaning the rank of a campaign and not merely its presence could be accounted for in the model.
Relative Order-sensitivity: it may be that campaigns A and B only have an effect when one is exposed to campaign A before campaign B. If so, it could be useful to assess the effect of given combinations of campaigns as well. And this for all campaigns, leading to tremendous numbers of possible combinations.
All previous phenomenon may be present, increasing even more the potential complexity of a comprehensive Attribution Model. The number of all possible ordered combination of k campaigns is indeed :

3.2 Main models

A) Logistic Regression and Classification models

If non converting journeys are available, Attribition Model can be shaped as a simple classification issue. Campaign types or campaigns combination and volume of campaign types can be included in the model along with customer or time variables. As we are interested in inference (on campaigns effect) whether than prediction, a parametric model should be used, such as Logistic Regression. Non paramatric models such as Random Forests or Neural Networks can also be used though the interpretation of campaigns value would be more difficult to derive from the model results.

A common pitfall is the usual issue of spurious correlations on one hand and the correct interpretation of coefficients in business terms.

An advantage if the possibility to evaluate the relevance of the model using common model validation methods to evaluate its predictive power (validation set \ AUC \pseudo R squared).

B) Shapley Value

Theory

The Shapley Value is based on a Game Theory framework and is named after its creator, the Nobel Price Laureate Lloyd Shapley. Initially meant to calculate the marginal contribution of players in cooperative games, the model has received much attention in research and industry and has lately been applied to marketing issues. This model is typically used by Google Adords and other ad bidding vendors. Campaigns or marketing channels are in this model seen as compementary players looking forward to increasing a given KPI.
Contrarily to Logistic Regressions, it is a non-parametric model. Contrarily to Markov Chains, all results are built using existing journeys, and not simulated ones.

Channels are considered to enter the game sequentially under a certain joining order. Shapley value try to The Shapley value of channel i is the weighted sum of the marginal values that channel i adds to all possible coalitions that don’t contain channel i.
In other words, the main logic is to analyse the difference of gains when a channel i is added after a coalition Ck of k channels, k<=n. We then sum all the marginal contributions over all possible ordered combination Ck of all campaigns excluding i, with k<=n-1.

Subsets framework

A first an most usual way to compute the Shapley Vaue is to consider that when a channel enters coalition, its additionnal value is the same irrelevant of the order in which previous channels have appeared. In other words, journeys (A>B>C) and (B>A>C) trigger the same gains.
Shapley value is computed as the gains associated to adding a channel i to a subset of channels, weighted by the number of (ordered) sequences that the (unordered) subset represents, summed up on all possible subsets of the total set of campaigns where the channel i is not present.
The Shapley value of the channel ???????? is then:

where |S| is the number of campaigns of a coalition S and the sum extends over all subsets S that do not not contain channel j. ????(????) is the value of the coalition S and ????(???? ∪ {????????}) the value of the coalition formed by adding ???????? to coalition S. ????(???? ∪ {????????}) − ????(????) is thus the marginal contribution of channel ???????? to the coalition S.

The formula can be rewritten and understood as:

This method is convenient when data on the gains of on all possible permutations of all unordered k subsets of the n campaigns are available. It is also more convenient if the order of campaigns prior to the introduction of a campaign is thought to have no impact.

Ordered sequences

Let us define ????((A>B)) as the value of the sequence A then B. What is we let ????((A>B)) be different from ????((B>A)) ?
This time we would need to sum over all possible permutation of the S campaigns present before ???????? and the N-(S+1) campaigns after ????????. Doing so we will sum over all possible orderings (i.e all permutations of the n campaigns of the grand coalition containing all campaigns) and we can remove the permutation coefficient s!(p-s+1)!.

This method is convenient when the order of channels prior to and after the introduction of another channel is assumed to have an impact. It is also necessary to possess data for all possible permutations of all k subsets of the n campaigns, and not only on all (unordered) k-subsets of the n campaigns, k<=n. In other words, one must know the gains of A, B, C, A>B, B>A, etc. to compute the Shapley Value.

Differences between the two approaches

We simulate an ordered case where the value for each ordered sequence k for k<=3 is known. We compare it to the usual Shapley value calculated based on known gains of unordered subsets of campaigns. So as to compare relevant values, we have built the gains matrix so that the gains of a subset A, B i.e ????({B,A}) is the average of the gains of ordered sequences made up with A and B (assuming the number of journeys where A>B equals the number of journeys where B>A, we have ????({B,A})=0.5( ????((A>B)) + ????((B>A)) ). We let the value of the grand coalition be different depending on the order of campaigns-keeping the constraints that it averages to the value used for the unordered case.

Note: mvA refers to the marginal value of A in a given sequence.
With traditionnal unordered coalitions:

With ordered sequences used to compute the marginal values:

We can see that the two approaches yield very different results. In the unordered case, the Shapley Value campaign C is the highest, culminating at 20, while A and B have the same Shapley Value mvA=mvB=15. In the ordered case, campaign A has the highest Shapley Value and all campaigns have different Shapley Values.

This example illustrates the inherent differences between the set and sequences approach to Shapley values. Real life data is more likely to resemble the ordered case as conversion probabilities may for any given set of campaigns be influenced by the order through which the campaigns appear.

Advantages

Shapley value has become popular in allocation problems in cooperative games because it is the unique allocation which satisfies different axioms:

Efficiency: Shaple Values of all channels add up to the total gains (here, orders) observed.
Symmetry: if channels A and B bring the same contribution to any coalition of campaigns, then their Shapley Value i sthe same
Null player: if a channel brings no additionnal gains to all coalitions, then its Shapley Value is zero
Strong monotony: the Shapley Value of a player increases weakly if all its marginal contributions increase weakly

These properties make the Shapley Value close to what we intuitively define as a fair attribution.

Issues

The Shapley Value is based on combinatory mathematics, and the number of possible coalitions and ordered sequences becomes huge when the number of campaigns increases.
If unordered, the Shapley Value assumes the contribution of campaign A is the same if followed by campaign B or by C.
If ordered, the number of combinations for which data must be available and sufficient is huge.
Channels rarely present or present in long journeys will be played down.
Generally, gains are supposed to grow with the number of players in the game. However, it is plausible that in the marketing context a journey with a high number of channels will not necessarily bring more orders than a journey with less channels involved.

References:

R package: GameTheoryAllocation

Article:
Zhao & al, 2018 “Shapley Value Methods for Attribution Modeling in Online Advertising “
https://link.springer.com/content/pdf/10.1007/s13278-017-0480-z.pdf
Courses: https://www.lamsade.dauphine.fr/~airiau/Teaching/CoopGames/2011/coopgames-7%5b8up%5d.pdf
Blogs: https://towardsdatascience.com/one-feature-attribution-method-to-supposedly-rule-them-all-shapley-values-f3e04534983d

B) Markov Chains

Markov Chains are used to model random processes, i.e events that occur in a sequential manner and in such a way that the probability to move to a certain state only depends on the past steps. The number of previous steps that are taken into account to model the transition probability is called the memory parameter of the sequence, and for the model to have a solution must be comprised between 0 and 4. A Markov Chain process is thus defined entirely by its Transition Matrix and its initial vector (i.e the starting point of the process).

Markov Chains are applied in many scientific fields. Typically, they are used in weather forecasting, with the sequence of Sunny and Rainy days following a Markov Process of memory parameter 0, so that for each given day the probability that the next day will be rainy or sunny only depends on the weather of the current day. Other applications can be found in sociology to understand the dynamics of social classes intergenerational reproduction. To get more both mathematical and applied illustration, I recommend the reading of this course.

In the marketing context, Markov Chains are an interesting way to model the conversion funnel. To go from the from the Markov Model to the Attribution logic, we calculate the Removal Effect of each channel, i.e the difference in conversions that happen if the channel is removed. Please read below for an introduction to the methodology.

The first step in a Markov Chains Attribution Model is to build the transition matrix that captures the transition probabilities between the campaigns accross existing customer journeys. This Matrix is to be read as a “From state A to state B” table, from the left to the right. A first difficulty is finding the right memory parameter to use. A large memory parameter would allow to take more into account interraction effects within the conversion funnel but would lead to increased computationnal time, a non-readable transition matrix, and be more sensitive to noisy data. Please note that this transition matrix provides useful information on the conversion funnel and on the relationships between campaigns and can be used as such as an analytical tool. I suggest the clear and easily R code which can be found here or here.

Here is an illustration of a Markov Chain with memory Parameter of 0: the probability to go to a certain campaign B in the next step only depend on the campaign we are currently at:

The associated Transition Matrix is then (with null probabilities left as Blank):

The second step is to compute the actual responsibility of a channel in total conversions. As mentionned above, the main philosophy to do so is to calculate the Removal Effect of each channel, i.e the changes in the number of conversions when a channel is entirely removed. All customer journeys which went through this channel are settled out to be unsuccessful. This calculation is done by applying the transition matrix with and without the removed channels to an initial vector that contains the number of desired simulations.

Building on our current example, we can then settle an initial vector with the desired number of simulations, e.g 10 000:

It is possible at this stage to add a constraint on the maximum number of times the matrix is applied to the data, i.e on the maximal number of campaigns a simulated journey is allowed to have.

Advantages

The dynamic journey is taken into account, as well as the transition between two states. The funnel is not assumed to be linear.
It is possile to build a conversion graph that maps the customer journey provides valuable insights.
It is possible to evaluate partly the accuracy of the Attribution Model based on Markov Chains. It is for example possible to see how well the transition matrix help predict the future by analysing the number of correct predictions at any given step over all sequences.

Disadvantages

It can be somewhat difficult to set the memory parameter. Complementarity effects between channels are not well taken into account if the memory is low, but a parameter too high will lead to over-sensitivity to noise in the data and be difficult to implement if customer journeys tend to have a number of campaigns below this memory parameter.
Long journeys with different channels involved will be overweighted, as they will count many times in the Removal Effect. For example, if there are n-1 channels in the customer journey, this journey will be considered as failure for the n-1 channel-RE. If the volume effects (i.e the impact of the overall number of channels in a journey, irrelevant from their type° are important then results may be biased.

References:

R package: ChannelAttribution

Git:

https://github.com/MatCyt/Markov-Chain/blob/master/README.md

Course:

https://www.ssc.wisc.edu/~jmontgom/markovchains.pdf

Article:

“Mapping the Customer Journey: A Graph-Based Framework for Online Attribution Modeling”; Anderl, Eva and Becker, Ingo and Wangenheim, Florian V. and Schumann, Jan Hendrik, 2014. Available at SSRN: https://ssrn.com/abstract=2343077 or http://dx.doi.org/10.2139/ssrn.2343077

“Media Exposure through the Funnel: A Model of Multi-Stage Attribution”, Abhishek & al, 2012

“Multichannel Marketing Attribution Using Markov Chains”, Kakalejčík, L., Bucko, J., Resende, P.A.A. and Ferencova, M. Journal of Applied Management and Investments, Vol. 7 No. 1, pp. 49-60. 2018

Blogs:

https://analyzecore.com/2016/08/03/attribution-model-r-part-1

https://analyzecore.com/2016/08/03/attribution-model-r-part-2

3.3 To go further: Tackling selection biases with Quasi-Experiments

Exposure to certain types of advertisement is usually highly correlated to non-observable variables. Differences in the behaviour of users exposed to different campaigns may thus only be driven by core differences in converison probabilities between groups whether than by the campaign effect. These potential selection effects may bias the results obtained using historical data.

Quasi-Experiments can help correct this selection effect while still using available observationnal data. These methods recreate the settings on a randomized setting. The goal is to come as close as possible to the ideal of comparing two populations that are identical in all respects except for the advertising exposure. However, populations might still diﬀer with respect to some unobserved characteristics.

Common quasi-experimental methods used for instance in Public Policy Evaluation are:

Discontinuity Regressions
Matching Methods, such as Exact Matching, Propensity-score matching or k-nearest neighbourghs.

References:

Article:

“Towards a digital Attribution Model: Measuring the impact of display advertising on online consumer behaviour”, Anindya Ghose & al, MIS Quarterly Vol. 40 No. 4, pp. 1-XX, 2016

https://pdfs.semanticscholar.org/4fa6/1c53f281fa63a9f0617fbd794d54911a2f84.pdf

4. First Steps towards a Practical Implementation

Identify key points of interests

Identify the nature of touchpoints available: is the data based on clicks? If so, is there a way to complement the data with A/B tests to measure the influence of ads without clicks (display, video) ? For example, what happens to sales when display campaign is removed? Analysing this multiplier effect would give the overall responsibility of display on sales, to be deduced from current attribution values given to click-based channels. More interestingly, what is the impact of the removal of display campaign on the occurences of click-based campaigns ? This would give us an idea of the impact of display ads on the exposure to each other campaigns, which would help correct the attribution values more precisely at the campaign level.
Define the KPI to track. From a pure Marketing perspective, looking at purchases may be sufficient, but from a financial perspective looking at profits, though a bit more difficult to compute, may drive more interesting results.
Define a customer journey. It may seem obvious, but the notion needs to be clarified at first. Would it be defined by a time limit? If so, which one? Does it end when a conversion is observed? For example, if a customer makes 2 purchases, would the campaigns he’s been exposed to before the first order still be accounted for in the second order? If so, with a time decay?
Define the research framework: are we interested only in customer journeys which have led to conversions or in all journeys? Keep in mind that successful customer journeys are a non-representative sample of customer journeys. Models built on the analysis of biased samples may be conservative. Take an extreme example: 80% of customers who see campaign A buy the product, VS 1% for campaign B. However, campaign B exposure is great and 100 Million people see it VS only 1M for campaign A. An Attribution Model based on successful journeys will give higher credit to campaign B which is an auguable conclusion. Taking into account costs per campaign (in the case where costs are calculated by clicks) may of course tackle this issue partly, as campaign A could then exhibit higher returns, but a serious fallacious reasonning is at stake here.

Analyse the typical customer journey

Performing a duration analysis on the data may help you improve the definition of the customer journey to be used by your organization. After which days are converison probabilities null? Should we consider the effect of campaigns disappears after x days without orders? For example, if 99% of orders are placed in the 30 days following a first click, it might be interesting to define the customer journey as a 30 days time frame following the first oder.
Look at the distribution of the number of campaigns in a typical journey. If you choose to calculate the effect of campaigns interraction in your Attribution Model, it may indeed help you determine the maximum number of campaigns to be included in a combination. Indeed, you may not need to assess the impact of channel combinations with above than 4 different channels if 95% of orders are placed after less then 4 campaigns.
Transition matrixes: what if a campaign A systematically leads to a campaign B? What happens if we remove A or B? These insights would give clues to ask precise questions for a latter AB test, for example to find out if there is complementarity between channels A and B – (implying none should be removed) or mere substitution (implying one can be given up).
If conversion rates are available: it can be interesting to perform a survival analysis i.e to analyse the likelihood of conversion based on duration since first click. This could help us excluse potential outliers or individuals who have very low conversion probabilities.

Summary

Attribution is a complex topic which will probably never be definitively solved. Indeed, a main issue is the difficulty, or even impossibility, to evaluate precisely the accuracy of the attribution model that we’ve built. Attribution Models should be seen as a good yet always improvable approximation of the incremental values of campaigns, and be presented with their intrinsinc limits and biases.

A common trap when it comes to sampling from a population that intrinsically includes outliers

February 21, 2019/0 Comments/in Business Analytics, Business Intelligence, Data Mining, Data Science, Data Science Hack, Main Category, R Statistics, Statistics /by Moussa El-Zaarah

I will discuss a common fallacy concerning the conclusions drawn from calculating a sample mean and a sample standard deviation and more importantly how to avoid it.

Suppose you draw a random sample $x_1$ , $x_2$ , … $x_N$ of size $N$ and compute the ordinary (arithmetic) sample mean $x_m$ and a sample standard deviation $sd$ from it. Now if (and only if) the (true) population mean µ (first moment) and population variance (second moment) obtained from the actual underlying PDF are finite, the numbers $x_m$ and $sd$ make the usual sense otherwise they are misleading as will be shown by an example.

By the way: The common correlation coefficient will also be undefined (or in practice always point to zero) in the presence of infinite population variances. Hopefully I will create an article discussing this related fallacy in the near future where a suitable generalization to Lévy-stable variables will be proposed.

Drawing a random sample from a heavy tailed distribution and discussing certain measures

As an example suppose you have a one dimensional random walker whose step length is distributed by a symmetric standard Cauchy distribution (Lorentz-profile) with heavy tails, i.e. an alpha-stable distribution with alpha being equal to one. The PDF of an individual independent step is given by $p(x) = \frac{\pi^{-1}}{(1 + x^2)}$ , thus neither the first nor the second moment exist whereby the first exists and vanishes at least in the sense of a principal value due to symmetry.

Still let us generate $N = 3000$ (pseudo) standard Cauchy random numbers in R* to analyze the behavior of their sample mean and standard deviation $sd$ as a function of the reduced sample size $n \leq N$ .

*The R-code is shown at the end of the article.

Here are the piecewise sample mean (in blue) and standard deviation (in red) for the mentioned Cauchy sampling. We see that both the sample mean and $sd$ include jumps and do not converge.

Especially the mean deviates relatively largely from zero even after 3000 observations. The sample $sd$ has no target due to the population variance being infinite.

If the data is new and no prior distribution is known, computing the sample mean and $sd$ will be misleading. Astonishingly enough the sample mean itself will have the (formally exact) same distribution as the single step length $p(x)$ . This means that the sample mean is also standard Cauchy distributed implying that with a different Cauchy sample one could have easily observed different sample means far of the presented values in blue.

What sense does it make to present the usual interval $x_m \pm sd / \sqrt{N}$ in such a case? What to do?

The sample median, median absolute difference (mad) and Inter-Quantile-Range (IQR) are more appropriate to describe such a data set including outliers intrinsically. To make this plausible I present the following plot, whereby the median is shown in black, the mad in green and the IQR in orange.

This example shows that the median, mad and IQR converge quickly against their assumed values and contain no major jumps. These quantities do an obviously better job in describing the sample. Even in the presence of outliers they remain robust, whereby the mad converges more quickly than the IQR. Note that a standard Cauchy sample will contain half of its sample in the interval median $\pm$ mad meaning that the IQR is twice the mad.

Drawing a random sample from a PDF that has finite moments

Just for comparison I also show the above quantities for a standard normal (pseudo) sample labeled with the same color as before as a counter example. In this case not only do both the sample mean and median but also the $sd$ and mad converge towards their expected values (see plot below). Here all the quantities describe the data set properly and there is no trap since there are no intrinsic outliers. The sample mean itself follows a standard normal, so that the $sd$ in deed makes sense and one could calculate a standard error $\frac{sd}{\sqrt{N}}$ from it to present the usual stochastic confidence intervals for the sample mean.

A careful observation shows that in contrast to the Cauchy case here the sampled mean and $sd$ converge more quickly than the sample median and the IQR. However still the sampled mad performs about as well as the $sd$ . Again the mad is twice the IQR.

And here are the graphs of the prementioned quantities for a pseudo normal sample:

The take-home-message:

Just be careful when you observe outliers and calculate sample quantities right away, you might miss something. At best one carefully observes how the relevant quantities change with sample size as demonstrated in this article.

Such curves should become of broader interest in order to improve transparency in the Data Science process and reduce fallacies as well.

Thank you for reading.

P.S.: Feel free to play with the set random seed in the R-code below and observe how other quantities behave with rising sample size. Of course you can also try different PDFs at the beginning of the code. You can employ a Cauchy, Gaussian, uniform, exponential or Holtsmark (pseudo) random sample.

QUIZ: Which one of the recently mentioned random samples contains a trap** and why?

**in the context of this article

R-code used to generate the data and for producing plots:

#R-script for emphasizing convergence and divergence of sample means

####install and load relevant packages ####

#uncomment these lines if necessary

#install.packages(c('ggplot2',’stabledist’))

#library(ggplot2)

#library(stabledist)

#####drawing random samples #####

#Setting a random seed for being able to reproduce results

set.seed(1234567)

N= 2000 #sample size

#Choose a PDF from which a sample shall be drawn

#To do so (un)comment the respective lines of following code

data <- rcauchy(N) # option1(default): standard Cauchy sampling

#data <- rnorm(N) #option2: standard Gaussian sampling

#data <- rexp(N) # option3: standard exponential sampling

#data <- rstable(N,alpha=1.5,beta=0) # option4: standard symmetric Holtsmark sampling

#data <- runif(N) #option5: standard uniform sample

#####descriptive statistics####

#preparations/declarations

SUM = vector()

sd =vector()

mean = vector()

SQ =vector()

SQUARES = vector()

median = vector()

mad =vector()

quantiles = data.frame()

sem =vector()

#piecewise calculaion of descrptive quantities

for (k in 1:length(data)){ #mainloop

SUM[k] <- sum(data[1:k]) # sum of sample

mean[k] <- mean(data[1:k]) # arithmetic mean

sd[k] <- sd(data[1:k]) # standard deviation

sem[k] <- sd[k]/(sqrt(k)) #standard error of the sample mean (for finite variances)

mad[k] <- mad(data[1:k],const=1) # median absolute deviation

for (j in 1:5){

qq <- quantile(data[1:k],na.rm = T)

quantiles[k,j] <- qq[j] #quantiles of sample

}

colnames(quantiles) <- c('min','Q1','median','Q3','max')

for (i in 1:length(data[1:k])){

SQUARES[i] <- data[i]*data[i]

}

SQ[k] <- sum(SQUARES[1:k]) #sum of squares of random sample

} #end of mainloop

#create table containing all relevant data

TABLE <- as.data.frame(cbind(quantiles,mean,sd,SQ,SUM,sem))

#####plotting results###

x11()

print(ggplot(TABLE,aes(1:N,median))+

geom_point(size=.5)+xlab('sample size n')+ylab('sample median'))

x11()

print(ggplot(TABLE,aes(1:N,mad))+geom_point(size=.5,color ='green')+

xlab('sample size n')+ylab('sample median absolute difference'))

x11()

print(ggplot(TABLE,aes(1:N,sd))+geom_point(size=.5,color ='red')+

xlab('sample size n')+ylab('sample standard deviation'))

x11()

print(ggplot(TABLE,aes(1:N,mean))+geom_point(size=.5, color ='blue')+

xlab('sample size n')+ylab('sample mean'))

x11()

print(ggplot(TABLE,aes(1:N,Q3-Q1))+geom_point(size=.5, color ='blue')+

xlab('sample size n')+ylab('IQR'))

#uncomment the following lines of code to see further plots

#x11()

#print(ggplot(TABLE,aes(1:N,sem))+geom_point(size=.5)+

#xlab('sample size n')+ylab('sample sum of r.v.'))

#x11()

#print(ggplot(TABLE,aes(1:N,SUM))+geom_point(size=.5)+

#xlab('sample size n')+ylab('sample sum of r.v.'))

#x11()

#print(ggplot(TABLE,aes(1:N,SQ))+geom_point(size=.5)+

#xlab('sample size n')+ylab('sample sum of squares'))

Modelling Data – Case Study: Importance of domain knowledge

August 29, 2018/2 Comments/in Artificial Intelligence, Audit Analytics, Big Data, Data Science, Insights, Machine Learning, Realtime Analytics, Statistics, Use Case /by Marius Schulz

What´s the relation between earnings and happiness? I saw this chart and was strongly irritated – why is there a linear regression, it´s clearly a logarithmic relationship.

So I got angry and wanted to know, which model is the better fit. I started to work immediatly, because it´s a huge difference for man kind. Think about it: you give a poor person money and he gets as happy as a rich person with the same amount added – that´s against common sense and propaganda to get rich. Like an cultural desease.

So I gathered the data and did a first comparation, and this logarithmic model was the better fit:

I was right and seriously willing to clear the mess up – so posted the “correct” model on facebook, to explain things to my friends.

Once I came down…

I asked myself: “What´s the model that fits the data best – that would be more correct?”

So I started to write an algorithm to check polynominal regression levels for fit using a random train and test data split. Finally, I got to this result and was amazed:

This seriously hit me: “What the f***! There seems to be maximum happiness reachable with a certain amount of income / GDP.” Can you understand, what this result would mean for our world and economy? Think about all economies growing continiously, but well happiest was there or will come there. What would you do? Send income to less developed countries, because you don´t need it? Stop invention and progress, because it´s of no use? Seriously, I felt like a socialist: Stop progress at this point and share.

So I thought a while and concluded: “F***ing statistics, we need a profound econometric model.”

I started modelling: Well, the first amount of money in a market based on money leverages a huge amount of happiness, because you can participate and feed yourself. We can approximate that by infinit marginal utility. Then the more you have, the less utility should be provided by the additional same amount added. Finally, more income is more options, so more should be always better. I concluded, that this is catched by a Cobb Douglas production function. Here´s the graph:

That´s it, that´s the final model. Here I feel home, this looks like a normal world – for an economist.

The Relevance of Domain Knowledge

As this short case study shows, we get completly wrong information and conclusions, if we don´t do it right. If you were the most important decision making algorithm in global economic politics, imagine what desasterous outcomes it would have produced to automatically find an optimum of income.

This is a serious border of AI. If you want to analyse Big Data with algorithms, you may produce seriously wrong information and conclusions. Statistical analysis is allways about using the right model. And modelling is about the assumptions of the model. As long as you can not create the right assumtions for the statistical model automatically, Big Data analysis is near to crazy. So out of this point of view, Big Data analysis is either about very simplistic tendencies (like linear trends) or it´s bound to Data Scientists with domain knowledge checking each model – that´s slow.

Discussion

I´m quite new to the field of Data Science, but this case study shows very though limitations, clearly. It´s not about flexible fitting of data, it´s about right models. And right models don´t scale into the Big Data domain. What do you think is the solution for this issue?

Countries of Happiness – the Full Article

If you are interested in my final article on my personal blog, explaining the final results: Please feel welcome to read the article here. There is a translation widget in the menu, to read in your favorite language. The original article is german.

ID3-Algorithmus: Ein Rechenbeispiel

December 8, 2017/4 Comments/in Artificial Intelligence, Business Analytics, Customer Analytics, Data Mining, Data Science, Machine Learning, Main Category, Predictive Analytics, Statistics, Visualization /by Benjamin Aunkofer

Dieser Artikel ist Teil 3 von 4 der Artikelserie Maschinelles Lernen mit Entscheidungsbaumverfahren und nun wollen wir einen Entscheidungsbaum aus Daten herleiten, jedoch ohne Programmierung, sondern direkt auf Papier (bzw. HTML :-).

Folgender Datensatz sei gegeben:

Zeile	Kundenart	Zahlungsgeschwindigkeit	Kauffrequenz	Herkunft	Zahlungsmittel: Rechnung?
1	Neukunde	niedrig	niedrig	Inland	false
2	Neukunde	niedrig	niedrig	Ausland	false
3	Stammkunde	niedrig	niedrig	Inland	true
4	Normalkunde	mittel	niedrig	Inland	true
5	Normalkunde	hoch	hoch	Inland	true
6	Normalkunde	hoch	hoch	Ausland	false
7	Stammkunde	hoch	hoch	Ausland	true
8	Neukunde	mittel	niedrig	Inland	false
9	Neukunde	hoch	hoch	Inland	true
10	Normalkunde	mittel	hoch	Inland	true
11	Neukunde	mittel	hoch	Ausland	true
12	Stammkunde	mittel	niedrig	Ausland	true
13	Stammkunde	niedrig	hoch	Inland	true
14	Normalkunde	mittel	niedrig	Ausland	false

Gleich vorweg ein Disclaimer: Der Datensatz ist natürlich überaus klein, ja gerade zu winzig. Dafür würden wir in der Praxis niemals einen Machine Learning Algorithmus einsetzen. Dennoch bleiben wir besser übersichtlich und nachvollziehbar mit diesen 14 Zeilen. Das Lernziel dieser Übung ist es, ein Gefühl für die Erstellung von Entscheidungsbäumen zu erhalten.
Zu beachten ist ferner, dass dieser Datensatz bereits aggregiert ist, denn eigentlich nummerisch abbildbare Daten wurden in Klassen zusammengefasst.

Das Ziel:

Der Datensatz spielt wieder, welchem Kunden (ID) bisher die Zahlung per Rechnung erlaubt und nicht widerrufen wurde. Das Ziel soll sein, eine Vorhersage darüber zu machen zu können, wann ein Kunde per Rechnung zahlen darf und wann nicht (dann per Vorkasse).

Der Algorithmus:

Wir verwenden den ID3-Algorithmus in seiner Reinform. Der ID3-Algorithmus ist der gängigste Algorithmus zum Aufbau datengetriebener Entscheidungsbäume und es gibt mehrere Abwandlungen. Die Vorgehensweise des Algorithmus wird in dem Teil 2 der Artikelserie Entscheidungsbaum-Algorithmus ID3 erläutert.

1. Schritt: Auswählen des Attributes mit dem höchsten Informationsgewinn

Der Informationsgewinn eines Attributes ( $A$ ) im Sinne des ID3-Algorithmus ist die Differenz aus der Entropie ( $E(S)$ ) (siehe Teil 1 der Artikelserie Entropie, ein Maß für die Unreinheit in Daten) des gesamten Datensatzes ( $S$ ) und der Summe aus den gewichteten Entropien des Attributes für jeden einzelnen Wert (Value $i$ ), der im Attribut vorkommt:
$IG(S, A) = H(S) - \sum_{i=1}^n \frac{\bigl|S_i\bigl|}{\bigl|S\bigl|} \cdot H(S_i)$

1.1 Gesamt-Entropie des Datensatzes berechnen

Erstmal schauen wir uns die Entropie des gesamten Datensatzes an. Die Entropie bezieht sich dabei auf das gewünschte Klassifikationsergebnis, also ist die Zahlung via Rechnung erlaubt oder nicht? Diese Frage wird entweder mit true oder false beantwortet.

$H(S) = - \frac{9}{14} \cdot \log_2(\frac{9}{14}) - \frac{5}{14} \cdot \log_2(\frac{5}{14}) = 0.94$

1.2 Berechnung der Informationsgewinne aller Attribute

Berechnen wir nun also die Informationsgewinne über alle Spalten.

Attribut	Subset	Count(true)	Count(false)
Kundenart	“Neukunde”	2	3
	“Stammkunde”	4	0
	“Normalkunde”	3	2

Wir zerlegen den gesamten Datensatz gedanklich in drei Kategorien der Kundenart und berechnen die Entropie bezogen auf das Klassifikationsziel:

$H(S_{Neukunde}) = - \frac{2}{5} \cdot \log_2(\frac{2}{5}) - \frac{3}{5} \cdot \log_2(\frac{3}{5}) = 0.97$

$H(S_{Stammkunde}) = - \frac{4}{4} \cdot \log_2(\frac{4}{4}) - \frac{0}{4} \cdot \log_2(\frac{0}{4}) = 0.00$

$H(S_{Normalkunde}) = - \frac{3}{5} \cdot \log_2(\frac{3}{5}) - \frac{2}{5} \cdot \log_2(\frac{2}{5}) = 0.97$

Zur Erinnerung, der Informationsgewinn (Information Gain) wird wie folgt berechnet:

$IG(S, A_{Kundenart}) = - \sum_{i=1}^n \frac{\bigl|S_i\bigl|}{\bigl|S\bigl|} \cdot H(S_i)$

Angewendet auf das Attribut “Kundenart”…

$IG(S, A_{Kundenart}) = H(S) - \frac{\bigl|S_{Neukunde}\bigl|}{\bigl|S\bigl|} \cdot H(S_{Neukunde}) - \frac{\bigl|S_{Stammkunde}\bigl|}{\bigl|S\bigl|} \cdot H(S_{Stammkunde}) - \frac{\bigl|S_{Normalkunde}\bigl|}{\bigl|S\bigl|} \cdot H(S_{Normalkunde})$

… erhalten wir der Formal nach folgenden Informationsgewinn:

$IG(S, A_{Kundenart}) = 0.94 - \frac{5}{14} \cdot 0.97 - \frac{4}{14} \cdot 0.00 - \frac{5}{14} \cdot 0.97 = 0.247$

Nun für die weiteren Spalten:

Attribut	Subset	Count(true)	Count(false)
Zahlungsgeschwindigkeit	“niedrig”	2	2
	“mittel”	4	2
	“schnell”	3	1

Entropien für die “Zahlungsgeschwindigkeit”:

$H(S_{niedrig}) = - \frac{2}{4} \cdot \log_2(\frac{2}{4}) - \frac{2}{4} \cdot \log_2(\frac{2}{4}) = 1.00$

$H(S_{mittel}) = - \frac{4}{6} \cdot \log_2(\frac{4}{6}) - \frac{2}{6} \cdot \log_2(\frac{2}{6}) = 0.92$

$H(S_{schnell}) = - \frac{3}{4} \cdot \log_2(\frac{3}{4}) - \frac{1}{4} \cdot \log_2(\frac{1}{4}) = 0.81$

So berechnen wir wieder den Informationsgewinn:

$IG(S, A_{Zahlungsgeschwindigkeit}) = H(S) - \frac{\bigl|S_{niedrig}\bigl|}{\bigl|S\bigl|} \cdot H(S_{niedrig}) - \frac{\bigl|S_{mittel}\bigl|}{\bigl|S\bigl|} \cdot H(S_{mittel}) - \frac{\bigl|S_{schnell}\bigl|}{\bigl|S\bigl|} \cdot H(S_{schnell})$

Einsatzen und ausrechnen:

$IG(S, A_{Zahlungsgeschwindigkeit}) = 0.94 - \frac{4}{14} \cdot 1.00 - \frac{6}{14} \cdot 0.92 - \frac{4}{14} \cdot 0.81 = 0.029$

Und nun für die Spalte “Kauffrequenz”:

Attribut	Subset	Count(true)	Count(false)
Kauffrequenz	“niedrig”	3	4
Kauffrequenz	“hoch”	6	1

Entropien:

$H(S_{niedrig}) = - \frac{3}{7} \cdot \log_2(\frac{3}{7}) - \frac{4}{7} \cdot \log_2(\frac{4}{7}) = 0.99$

$H(S_{hoch}) = - \frac{6}{7} \cdot \log_2(\frac{6}{7}) - \frac{1}{7} \cdot \log_2(\frac{1}{7}) = 0.59$

Informationsgewinn:

$IG(S, A_{Kauffrequenz}) = H(S) - \frac{\bigl|S_{niedrig}\bigl|}{\bigl|S\bigl|} \cdot H(S_{niedrig}) - \frac{\bigl|S_{hoch}\bigl|}{\bigl|S\bigl|} \cdot H(S_{hoch})$

Einsetzen und Ausrechnen:

$IG(S, A_{Kauffrequenz}) = 0.94 - \frac{7}{14} \cdot 1.00 - \frac{7}{14} \cdot 0.59 = 0.150$

Und last but not least die Spalte “Herkunft”:

Attribut	Subset	Count(true)	Count(false)
Herkunft	“Inland”	6	2
Herkunft	“Ausland”	3	3

Entropien:

$H(S_{Inland}) = - \frac{6}{8} \cdot \log_2(\frac{6}{8}) - \frac{2}{8} \cdot \log_2(\frac{2}{8}) = 0.81$

$H(S_{Ausland}) = - \frac{3}{6} \cdot \log_2(\frac{3}{6}) - \frac{3}{6} \cdot \log_2(\frac{3}{6}) = 1.00$

Informationsgewinn:

$IG(S, A_{Herkunft}) = H(S) - \frac{\bigl|S_{Inland}\bigl|}{\bigl|S\bigl|} \cdot H(S_{Inland}) - \frac{\bigl|S_{Ausland}\bigl|}{\bigl|S\bigl|} \cdot H(S_{Ausland})$

Einsetzen und Ausrechnen:

$IG(S, A_{Herkunft}) = 0.94 - \frac{8}{14} \cdot 0.81 - \frac{6}{14} \cdot 1.00 = 0.05$

2. Schritt: Anlegen des Wurzel-Knotens

Der Informationsgewinn ist für das Attribut “Kundenart” am größten, daher entscheiden wir uns im Sinne des ID3-Algorithmus für dieses Attribut als Wurzel-Knoten.

3. Schritt: Rekursive Wiederholung (!!!)

Nun stellt sich natürlich die Frage: Wie geht es weiter?

Der Algorithmus kann eigentlich nur eines: Einen Wurzelknoten finden. Diesen Vorgang müssen wir nun nur noch rekursiv wiederholen, und das tun wir wie folgt.

Der Datensatz wurde bereits aufgeteilt in die drei Kundenarten. Für jede Kundenart ergibt sich jeweils ein Subset mit den verbleibenden Attributen. Für alle drei Subsets erstellen wir dann wieder einen Wurzelknoten, so dass ein neuer Ast entsteht.

3.1 Erster Rekursionsschritt

Machen wir also weiter und bestimmen wir das nächste Attribut nach der Kundenart, für die Fälle Kundenart = “Neukunde”:

Zeile	Kundenart	Zahlungsgeschwindigkeit	Kauffrequenz	Herkunft	Zahlungsmittel: Rechnung?
1	Neukunde	niedrig	niedrig	Inland	false
2	Neukunde	niedrig	niedrig	Ausland	false
8	Neukunde	mittel	niedrig	Inland	false
9	Neukunde	hoch	hoch	Inland	true
11	Neukunde	mittel	hoch	Ausland	true

Die Entropie des Gesamtdatensatzes (ja, es ist für diesen Schritt betrachtet der gesamte Datensatz!) ist wie folgt:

$H(S_{Neukunde}) = - \frac{2}{5} \cdot \log_2(\frac{2}{5}) - \frac{3}{5} \cdot \log_2(\frac{3}{5}) = 0.97$

Die Entropie ist weit weg von einer bestimmten Wahrscheinlichkeit (nahe der Gleichverteilung). Daher müssen wir hier nochmal ansetzen und losrechnen:

Entropien für “Zahlungsgeschwindigkeit” bei Neukunden:

$H(S_{niedrig}) = 0.00$

$H(S_{mittel}) = 1.00$

$H(S_{hoch}) = 0.00$

Informationsgewinn des Attributes “Zahlungsgeschwindigkeit” bei Neukunden:

$IG(S_{Neukunde},A_{Zahlungsgeschwindigkeit}) = 0.97 - \frac{3}{5} \cdot 0.00 - \frac{2}{5} \cdot 1.00 - \frac{1}{5} \cdot 0.00 = 0.57$

Betrachtung der Spalte “Kauffrequenz” bei Neukunden:

Entropien für “Kauffrequenz” bei Neukunden:

$H(S_{niedrig}) = 0.00$

$H(S_{hoch}) = 0.00$

Informationsgewinn des Attributes “Kauffrequenz” bei Neukunden:

$IG(S_{Neukunde},A_{Kauffrequenz}) = 0.97 - \frac{3}{5} \cdot 0.00 - \frac{2}{5} \cdot 0.00 = 0.97$

Betrachtung der Spalte “Herkunft” bei Neukunden:

Entropien für “Herkunft” bei Neukunden:

$H(S_{Inland}) = 0.92$

$H(S_{hoch}) = 1.00$

Informationsgewinn des Attributes “Herkunft” bei Neukunden:

$IG(S_{Neukunde},A_{Herkunft}) = 0.97 - \frac{3}{5} \cdot 0.92 - \frac{2}{5} \cdot 1.00 = 0.018$

Wir entscheiden uns also für das Attribut “Kauffrequenz” als Ast nach der Entscheidung “Neukunde”, denn dieses Attribut bring uns den größten Informationsgewinn und trennt uns die Unterscheidung für oder gegen das Zahlungsmittel “Rechnung” eindeutig auf.

3.1 Zweiter Rekursionsschritt

Was passiert mit der Kundenart “Stammkunde”?

Zeile	Kundenart	Zahlungsgeschwindigkeit	Kauffrequenz	Herkunft	Zahlungsmittel: Rechnung?
3	Stammkunde	niedrig	niedrig	Inland	true
7	Stammkunde	hoch	hoch	Ausland	true
12	Stammkunde	mittel	niedrig	Ausland	true
13	Stammkunde	niedrig	hoch	Inland	true

Die Antwort ist einfach: Nichts!
Wer ein Stammkunde ist, dem wurde stets die Zahlung per Rechnung erlaubt.

$H(S_{Stammkunde}) = 0.0$

3.1 Dritter Rekursionsschritt

Fehlt nun nur noch die Frage nach der Unterscheidung von Normalkunden.

Zeile	Kundenart	Zahlungsgeschwindigkeit	Kauffrequenz	Herkunft	Zahlungsmittel: Rechnung?
4	Normalkunde	mittel	niedrig	Inland	true
5	Normalkunde	hoch	hoch	Inland	true
6	Normalkunde	hoch	hoch	Ausland	false
14	Normalkunde	mittel	niedrig	Ausland	false

Zwar ist die Entropie des Subsets der Normalkunden…

$H(S_{Normalkunde}) = 1.0$

… denkbar schlecht, da maximal. Aber wir können genauso vorgehen, wie wir es bei dem Subset der Neukunden getan haben. Ich nehme es nun aber vorweg: Wenn wir uns den Datensatz näher ansehen, erkennen wir, dass wir diese Gesamtentropie von 1.0 für das Subset “Normalkunde” nicht mit den Attributen “Kauffrequenz” oder “Zahlungsgeschwindigkeit” reduzieren können, da dieses auch für sich betrachtet in Entropien der Größe 1.0 erhalten werden. Das Attribut “Herkunft” hingegen teilt den Datensatz sauber in true und false auf:

Somit ist der Informationsgewinn für das Attribut “Herkunft” am größten und wir haben unseren Baum komplett und – glücklicherweise – eindeutig bestimmen können!

Ergebnis: Der Entscheidungsbaum

Somit haben wir den Entscheidungsbaum über den ID3-Algorithmus erstellt, der eine Auskunft darüber macht, ob einem Kunden die Zahlung über Rechnung (statt Vorkasse) erlaubt wird:

true = Rechnung als Zahlungsmittel erlaubt
false = Rechnung als Zahlungsmittel nicht erlaubt

Lineare Regression in Python mit Scitkit-Learn

October 17, 2017/4 Comments/in Data Mining, Data Science, Data Science Hack, Python, Statistics, Visualization /by Benjamin Aunkofer

Die lineare Regressionsanalyse ist ein häufiger Einstieg ins maschinelle Lernen um stetige Werte vorherzusagen (Prediction bzw. Prädiktion). Hinter der Regression steht oftmals die Methode der kleinsten Fehlerquadrate und die hat mehr als eine mathematische Methode zur Lösungsfindung (Gradientenverfahren und Normalengleichung). Alternativ kann auch die Maximum Likelihood-Methode zur Regression verwendet werden. Wir wollen uns in diesem Artikel nicht auf die Mathematik konzentrieren, sondern uns direkt an die Anwendung mit Python Scikit-Learn machen:

Haupt-Lernziele:

Einführung in Machine Learning mit Scikit-Learn
Lineare Regression mit Scikit-Learn

Neben-Lernziele:

Datenvorbereitung (Data Preparation) mit Pandas und Scikit-Learn
Datenvisualisierung mit der Matplotlib direkt und indirekt (über Pandas)

Was wir inhaltlich tun:

Der Versuch einer Vorhersage eines Fahrzeugpreises auf Basis einer quantitativ-messbaren Eigenschaft eines Fahrzeuges.

Die Daten als Download

Für dieses Beispiel verwende ich die Datei “Automobil_data.txt” von Kaggle.com. Die Daten lassen sich über folgenden Link downloaden, nur leider wird ein (kostenloser) Account benötigt:
https://www.kaggle.com/toramky/automobile-dataset/downloads/automobile-dataset.zip
Sollte der Download-Link unerwartet mal nicht mehr funktionieren, freue ich mich über einen Hinweis als Kommentar 🙂

Die Entwicklungsumgebung

Ich verwende hier die Python-Distribution Anaconda 3 und als Entwicklungs-Umgebung Spyder (in Anaconda enthalten). Genauso gut funktionieren jedoch auch Jupyter Notebook, Eclipse mit PyDev oder direkt die IPython QT-Console.

Zuerst einmal müssen wir die Daten in unsere Python-Session laden und werden einige Transformationen durchführen müssen. Wir starten zunächst mit dem Importieren von drei Bibliotheken NumPy und Pandas, deren Bedeutung ich nicht weiter erläutern werde, somit voraussetze.

import matplotlib.pyplot as plt # Die Nr.1 der Bibliotheken zur Datenvisualisierung

import numpy as np # Bibliothek "Nummerisches Python"

import pandas as pd # Bibliothek "Panel Data"

Wir nutzen die Pandas-Bibliothek, um die “Automobile_data.txt” in ein pd.DataFrame zu laden.

dataSet = pd.read_csv("Automobile_data.txt", # Hier liegt die Datei im selben Verzeichnis wie das Python-Skript!

delimiter = ',',

thousands = None,

decimal = '.')

Schauen wir uns dann die ersten fünf Zeilen in IPython via dataSet.head().

In : dataSet.head()

Out:

symboling normalized-losses make fuel-type aspiration num-of-doors \

0 3 ? alfa-romero gas std two

1 3 ? alfa-romero gas std two

2 1 ? alfa-romero gas std two

3 2 164 audi gas std four

4 2 164 audi gas std four

body-style drive-wheels engine-location wheel-base ... engine-size \

0 convertible rwd front 88.6 ... 130

1 convertible rwd front 88.6 ... 130

2 hatchback rwd front 94.5 ... 152

3 sedan fwd front 99.8 ... 109

4 sedan 4wd front 99.4 ... 136

fuel-system bore stroke compression-ratio horsepower peak-rpm city-mpg \

0 mpfi 3.47 2.68 9.0 111 5000 21

1 mpfi 3.47 2.68 9.0 111 5000 21

2 mpfi 2.68 3.47 9.0 154 5000 19

3 mpfi 3.19 3.4 10.0 102 5500 24

4 mpfi 3.19 3.4 8.0 115 5500 18

highway-mpg price

0 27 13495

1 27 16500

2 26 16500

3 30 13950

4 22 17450

[5 rows x 26 columns]

Hinweis: Der Datensatz hat viele Spalten, so dass diese in der Darstellung mit einem Backslash \ umgebrochen werden.

Gleich noch eine weitere Ausgabe dataSet.info(), die uns etwas über die Beschaffenheit der importierten Daten verrät:

In : dataSet.info()

RangeIndex: 205 entries, 0 to 204

Data columns (total 26 columns):

symboling 205 non-null int64

normalized-losses 205 non-null object

make 205 non-null object

fuel-type 205 non-null object

aspiration 205 non-null object

num-of-doors 205 non-null object

body-style 205 non-null object

drive-wheels 205 non-null object

engine-location 205 non-null object

wheel-base 205 non-null float64

length 205 non-null float64

width 205 non-null float64

height 205 non-null float64

curb-weight 205 non-null int64

engine-type 205 non-null object

num-of-cylinders 205 non-null object

engine-size 205 non-null int64

fuel-system 205 non-null object

bore 205 non-null object

stroke 205 non-null object

compression-ratio 205 non-null float64

horsepower 205 non-null object

peak-rpm 205 non-null object

city-mpg 205 non-null int64

highway-mpg 205 non-null int64

price 205 non-null object

dtypes: float64(5), int64(5), object(16)

memory usage: 41.7+ KB

Einige Spalten entsprechen hinsichtlich des Datentypes nicht der Erwartung. Für die Spalten ‘horsepower’ und ‘peak-rpm’ würde ich eine Ganzzahl (Integer) erwarten, für ‘price’ hingegen eine Fließkommazahl (Float), allerdings sind die drei Spalten als Object deklariert. Mit Trick 17 im Data Science, der Anzeige der Minimum- und Maximum-Werte einer zu untersuchenden Datenreihe, kommen wir dem Übeltäter schnell auf die Schliche:

dataSet['horsepower'].min()

Out: '100'

dataSet['horsepower'].max()

Out: '?'

Datenbereinigung

Für eine Regressionsanalyse benötigen wir nummerische Werte (intervall- oder ratioskaliert), diese möchten wir auch durch richtige Datentypen-Deklaration herstellen. Nun wird eine Konvertierung in den gewünschten Datentyp jedoch an den (mit ‘?’ aufgefüllten) Datenlücken scheitern.

Schauen wir uns doch einmal die Datenreihen an, in denen in der Spalte ‘peak-rpm’ Fragezeichen stehen:

dataSet[dataSet['peak-rpm'] == '?'][['engine-type', 'num-of-cylinders']]

Out:

engine-type num-of-cylinders

130 ohc four

131 ohc four

Zwei Datenreihen sind vorhanden, bei denen ‘peak-rpm’ mit einem ‘?’ aufgefüllt wurde. Nun könnten wir diese Datenreihen einfach rauslöschen. Oder mit sinnvollen (im Sinne von wahrscheinlichen) Werten auffüllen. Vermutlichen haben beide Einträge – beide sind OHC-Motoren mit 4 Zylindern – eine ähnliche Drehzahl-Angabe wie vergleichbare Motoren. Mit folgendem Quellcode, gruppieren wir die Spalten ‘engine-type’ und ‘num-of-cylinders’ und bilden für diese Klassen den arithmetischen Mittelwert (.mean()) für die ‘peak-rpm’.

dataSet_rpm = dataSet[dataSet['peak-rpm'] != '?'][['engine-type', 'num-of-cylinders','peak-rpm']]

dataSet_rpm['peak-rpm'] = dataSet_rpm['peak-rpm'].astype(float)

dataSet_rpm_grouped = dataSet_rpm.groupby(['engine-type', 'num-of-cylinders'])

dataSet_rpm_grouped['peak-rpm'].mean()

Und schauen wir uns das Ergebnis an:

dataSet_rpm_grouped['peak-rpm'].mean()

Out:

engine-type num-of-cylinders

dohc four 5700.000000 --

six 5050.000000

dohcv eight 5750.000000

l four 4668.181818

three 5100.000000

ohc five 5081.818182

four 5155.468750

six 4821.428571

ohcf four 4775.000000

six 5900.000000

ohcv eight 4625.000000

six 5212.500000

twelve 5000.000000

rotor two 6000.000000

Name: peak-rpm, dtype: float64

Ein Vier-Zylinder-OHC-Motor hat demnach durchschnittlich einen Drehzahl-Peak von 5155 Umdrehungen pro Minute. Ohne nun (fahrlässigerweise) auf die Verteilung in dieser Klasse zu achten, nehmen wir einfach diesen Schätzwert, um die zwei fehlende Datenpunkte zu ersetzen.

Wir möchten jedoch die Original-Daten erhalten und legen ein neues DataSet (dataSet_c) an, in welches wir die Korrekturen vornehmen:

1 2	dataSet_c = dataSet.copy() # das "c"-Anhängsel steht für "corrected"

Nun können wir die fehlenden Peak-RPM-Einträge mit unserem Schätzwert ersetzen:

1 2	dataSet_c.loc[dataSet_c['peak-rpm'] == '?', 'peak-rpm'] = 5155

Was bei einer Drehzahl-Angabe noch funktionieren mag, ist für anderen Spalten bereits etwas schwieriger: Die beiden Spalten ‘price’ und ‘horsepower’ sind ebenfalls vom Typ Object, da sie ‘?’ enthalten. Verzichten wir einfach auf die betroffenen Zeilen:

dataSet_c = dataSet_c[dataSet_c['price'] != '?'] # entsprechende Zeilen herausfiltern

dataSet_c['price'] = dataSet_c['price'].astype(float) # Typ-Konvertierung zu Float

dataSet_c = dataSet_c[dataSet_c.horsepower != '?'] # entsprechende Zeilen herausfiltern

dataSet_c['horsepower'] = dataSet_c['horsepower'].astype(float) # Typ-Konvertierung in Int

Datenvisualisierung mit Pandas

Wir wollen uns nicht lange vom eigentlichen Ziel ablenken, dennoch nutzen wir die Visualisierungsfähigkeiten der Pandas-Library (welche die Matplotlib inkludiert), um uns dann die Anzahlen an Einträgen nach Hersteller der Fahrzeuge (Spalte ‘make’) anzeigen zu lassen:

dataSet_grouped_make = dataSet_c.groupby('make')

dataSet_grouped_make['make'].count().plot(kind = 'bar', figsize = (10, 10))

plt.show() # Besser jedes Plot abschließen! Auch wenn es in Pandas entstanden ist.

Oder die durchschnittliche PS-Zahl nach Hersteller:

(dataSet_c.groupby('make'))['horsepower'].mean().plot(kind = 'barh',

title = 'Mean Horsepower',

figsize = (10, 10))

plt.show()

Vorbereitung der Regressionsanalyse

Nun kommen wir endlich zur Regressionsanalyse, die wir mit Scikit-Learn umsetzen möchten. Die Regressionsanalyse können wir nur mit intervall- oder ratioskalierten Datenspalten betreiben, daher beschränken wir uns auf diese. Die “price”-Spalte nehmen wir jedoch heraus und setzen sie als unsere Zielgröße fest.

""" ----- Vorbereitung für die Regressionsanalyse ----- """

cols_ratio = ['horsepower', 'wheel-base', 'length', 'width', 'height', 'curb-weight', 'engine-size', 'compression-ratio', 'city-mpg', 'highway-mpg']

cols_target = ['price']

dataSet_ratio = dataSet_c.loc[:, cols_ratio]

dataSet_target = dataSet_c[cols_target]

Interessant ist zudem die Betrachtung vorab, wie die einzelnen nummerischen Attribute untereinander korrelieren. Dafür nehmen wir auch die ‘price’-Spalte wieder in die Betrachtung hinein und hinterlegen auch eine Farbskala mit dem Preis (höhere Preise, hellere Farben).

grr = pd.plotting.scatter_matrix(dataSet_c[cols_target + cols_ratio]

,c = dataSet_target

,figsize=(15, 15)

,marker = 'o'

,hist_kwds={'bins' : 20}

,s = 60

,alpha = 0.8)

plt.show()

Die lineare Korrelation ist hier sehr interessant, da wir auch nur eine lineare Regression beabsichtigen.

Wie man in dieser Scatter-Matrix recht gut erkennen kann, scheinen einige Größen-Paare nahezu perfekt zu korrelieren, andere nicht.

Korrelation…

…nahezu perfekt linear: highway-mpg vs city-mpg (mpg = Miles per Gallon)
… eher nicht gegeben: highway-mpg vs height
… nicht linear, dafür aber nicht-linear: highway-mpg vs price

Nun, wir wollen den Preis eines Fahrzeuges vorhersagen, wenn wir eine andere quantitative Größe gegeben haben. Auf den Preis bezogen, erscheint mir die Motorleistung (Horsepower) einigermaßen linear zu korrelieren. Versuchen wir hier die lineare Regression und setzen somit die Spalte ‘horsepower’ als X und ‘price’ als y fest.

X = dataSet_ratio[['horsepower']] # doppelte [], da eine Liste von Spalten zu übergeben ist

y = dataSet_c[cols_target]

Die gängige Konvention ist übrigens, X groß zu schreiben, weil hier auch mehrere x-Dimensionen enthalten sein dürfen (multivariate Regression). y hingegen, ist stets nur eine Zielgröße (eine Dimension).

Die lineare Regression ist ein überwachtes Verfahren des maschinellen Lernens, somit müssen wir unsere Prädiktionsergebnisse mit Test-Daten testen, die nicht für das Training verwendet werden dürfen. Scitkit-Learn (oder kurz: sklearn) bietet hierfür eine Funktion an, die uns das Aufteilen der Daten abnimmt:

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y,

test_size = 0.3, # 70% der Daten für das Training

random_state = None) # bei Bedarf kann hier "dem Zufall auf die Sprünge geholfen" werden

Zu beachten ist dabei, dass die Daten vor dem Aufteilen in Trainings- und Testdaten gut zu durchmischen sind. Auch dies übernimmt die train_test_split-Funktion für uns, nur sollte man im Hinterkopf behalten, dass die Ergebnisse (auf Grund der Zufallsauswahl) nach jedem Durchlauf immer wieder etwas anders aussehen.

Lineare Regression mit Scikit-Learn

Nun kommen wir zur Durchführung der linearen Regression mit Scitkit-Learn, die sich in drei Zeilen trainieren lässt:

""" ----- Lineare Regressionsanalyse ------- """

from sklearn.linear_model import LinearRegression # importieren der Klasse

lr = LinearRegression() # instanziieren der Klasse

lr.fit(X_train, y_train) # trainieren

Aber Vorsicht! Bevor wir eine Prädiktion durchführen, wollen wir festlegen, wie wir die Güte der Prädiktion bewerten wollen. Die gängigsten Messungen für eine lineare Regression sind der MSE und R².

$MSE = \frac{\sum_{i=1}^n (y_i - \hat{y_i})^2}{n}$

Ein großer MSE ist schlecht, ein kleiner gut.

$R^2 = 1 - \frac{MSE}{Var(y)}= \frac{\frac{1}{n} \cdot \sum_{i=1}^n (y_i - \hat{y_i})^2}{\frac{1}{n} \cdot \sum_{i=1}^n (y_i - \hat{\mu_y})^2}$

Ein kleines R² ist schlecht, ein großes R² gut. Ein R² = 1.0 wäre theoretisch perfekt (da der Fehler = 0.00 wäre), jedoch in der Praxis unmöglich, da dieser nur bei absolut perfekter Korrelation auftreten würde. Die Klasse LinearRegression hat eine R²-Messmethode implementiert (score(x, y)).

print('------ Lineare Regression -----')

print('Funktion via sklearn: y = %.3f * x + %.3f' % (lr.coef_[0], lr.intercept_))

print("Alpha: {}".format(lr.intercept_))

print("Beta: {}".format(lr.coef_[0]))

print("Training Set R² Score: {:.2f}".format(lr.score(X_train, y_train)))

print("Test Set R² Score: {:.2f}".format(lr.score(X_test, y_test)))

print("\n")

Die Ausgabe (ein Beispiel!):

------ Lineare Regression -----

Funktion via sklearn: y = 170.919 * x + -4254.701 # Die Funktion ist als y = 171 * x - 4254.7

Alpha: [-4254.70114803] # y-Achsenschnitt bei x = 0

Beta: [ 170.91919086] # Steigung der Gerade

Training Set R² Score: 0.62

Test Set R² Score: 0.73

Nach jedem Durchlauf ändert sich mit der Datenaufteilung (train_test_split()) das Modell etwas und auch R² schwankt um eine gewisse Bandbreite. Berauschend sind die Ergebnisse dabei nicht, und wenn wir uns die Regressionsgerade einmal ansehen, wird auch klar, warum:

plt.figure(figsize=(10,10))

plt.scatter(X_train, y_train, color = 'blue') # Blaue Punkte sind Trainingsdaten

plt.scatter(X_test, y_test, color = 'green') # Grüne Punkte sind Testdaten

plt.plot(X_train, lr.predict(X_train), color = 'red') # Hier ensteht die Gerade (x, y) = (x, lr.predict(x)

plt.xlabel(X_train.columns[0])

plt.ylabel(cols_target[0])

plt.show()

Bei kleineren Leistungsbereichen, etwa bis 100 PS, ist die Preis-Varianz noch annehmbar gering, doch bei höheren Leistungsbereichen ist die Spannweite deutlich größer. (Nachträgliche Anmerkung vom 06.05.2018: relativ betrachtet, bleibt der Fehler über alle Wertebereiche ungefähr gleich [relativer Fehler]. Die absoluten Fehlerwerte haben jedoch bei größeren x-Werten so eine Varianz der möglichen y-Werte, dass keine befriedigenden Prädiktionen zu erwarten sind.)

Egal wie wir eine Gerade in diese Punktwolke legen, wir werden keine befriedigende Fehlergröße erhalten.

Nehmen wir einmal eine andere Spalte für X, bei der wir vor allem eine nicht-lineare Korrelation erkannt haben: “highway-mpg”

X = dataSet_ratio[['highway-mpg']]

y = dataSet_c[cols_target]

Wenn wir dann das Training wiederholen:

------ Lineare Regression -----

Funktion via sklearn: y = -868.787 * x + 40575.036

Alpha: [ 40575.03556055]

Beta: [-868.7869183]

Training Set R² Score: 0.49

Test Set R² Score: 0.40

Die R²-Werte sind nicht gerade berauschend, und das erklärt sich auch leicht, wenn wir die Trainings- und Testdaten sowie die gelernte Funktionsgerade visualisieren:

Die Gerade lässt sich nicht wirklich gut durch diese Punktwolke legen, da letztere eher eine Kurve als eine Gerade bildet. Im Grunde könnte eine Gerade noch einigermaßen gut in den Bereich von 22 bis 43 mpg passen und vermutlich annehmbare Ergebnisse liefern. Die Wertebereiche darunter und darüber jedoch verzerren zu sehr und sorgen zudem dafür, dass die Gerade auch innerhalb des mittleren Bereiches zu weit nach oben verschoben ist (ggf. könnte hier eine Ridge-/Lasso-Regression helfen).

Richtig gute Vorhersagen über nicht-lineare Verhältnisse können jedoch nur mit einer nicht-linearen Regression erreicht werden.

Nicht-lineare Regression mit Scikit-Learn

Nicht-lineare Regressionsanalysen erlauben es uns, nicht-lineare korrelierende Werte-Paare als Funktion zu erlernen. Im folgenden Scatter-Plot sehen wir zum einen die gewohnte lineare Regressionsgerade (y = a * x + b) in rot, eine polinominale Regressionskurve dritten Grades (y = a * x³ + b * x² + c * x + d) in violet sowie einen Entscheidungsweg einer Entscheidungsbaum-Regression in gelb.

Nicht-lineare Regressionsanalysen passen sich dem Verlauf der Punktwolke sehr viel besser an und können somit in der Regel auch sehr gute Vorhersageergebnisse liefern. Ich ziehe hier nun jedoch einen Gedankenstrich, liefere aber den Quellcode für die lineare Regression als auch für die beiden nicht-linearen Regressionen mit:

Python Script Regression via Scikit-Learn

Weitere Anmerkungen

Bibliotheken wie Scitkit-Learn erlauben es, machinelle Lernverfahren schnell und unkompliziert anwenden zu können. Allerdings sollte man auch verstehen, wei diese Verfahren im Hintergrund mathematisch arbeiten. Diese Bibliotheken befreien uns also nicht gänzlich von der grauen Theorie.
Statt der “reinen” lineare Regression (LinearRegression()) können auch eine Ridge-Regression (Ridge()), Lasso-Regression (Lasso()) oder eine Kombination aus beiden als sogenannte ElasticNet-Regression (ElasticNet()). Bei diesen kann über Parametern gesteuert werden, wie stark Ausreißer in den Daten berücksichtigt werden sollen.
Vor einer Regression sollten die Werte skaliert werden, idealerweise durch Standardisierung der Werte (sklearn.preprocessing.StandardScaler()) oder durch Normierung (sklearn.preprocessing.Normalizer()).
Wir haben hier nur zwei-dimensional betrachtet. In der Praxis ist das jedoch selten ausreichend, auch der Fahrzeug-Preis ist weder von der Motor-Leistung, noch von dem Kraftstoffverbrauch alleine abhängig – Es nehmen viele Größen auf den Preis Einfluss, somit benötigen wir multivariate Regressionsanalysen.

Unsupervised Learning in R: K-Means Clustering

May 15, 2017/0 Comments/in Data Mining, Data Science, Machine Learning, R Statistics, Statistics, Tutorial /by Markus Lang

Die Clusteranalyse ist ein gruppenbildendes Verfahren, mit dem Objekte Gruppen – sogenannten Clustern zuordnet werden. Die dem Cluster zugeordneten Objekte sollen möglichst homogen sein, wohingegen die Objekte, die unterschiedlichen Clustern zugeordnet werden möglichst heterogen sein sollen. Dieses Verfahren wird z.B. im Marketing bei der Zielgruppensegmentierung, um Angebote entsprechend anzupassen oder im User Experience Bereich zur Identifikation sog. Personas.

Es gibt in der Praxis eine Vielzahl von Cluster-Verfahren, eine der bekanntesten und gebräuchlichsten Verfahren ist das K-Means Clustering, ein sog. Partitionierendes Clusterverfahren. Das Ziel dabei ist es, den Datensatz in K Cluster zu unterteilen. Dabei werden zunächst K beliebige Punkte als Anfangszentren (sog. Zentroiden) ausgewählt und jedem dieser Punkte der Punkt zugeordnet, zu dessen Zentrum er die geringste Distanz hat. K-Means ist ein „harter“ Clusteralgorithmus, d.h. jede Beobachtung wird genau einem Cluster zugeordnet. Zur Berechnung existieren verschiedene Distanzmaße. Das gebräuchlichste Distanzmaß ist die quadrierte euklidische Distanz:

$D^2 = \sum_{i=1}^{v}(x_i - y_i)^2$

Nachdem jede Beobachtung einem Cluster zugeordnet wurde, wird das Clusterzentrum neu berechnet und die Punkte werden den neuen Clusterzentren erneut zugeordnet. Dieser Vorgang wird so lange durchgeführt bis die Clusterzentren stabil sind oder eine vorher bestimmte Anzahl an Iterationen durchlaufen sind.
Das komplette Vorgehen wird im Folgenden anhand eines künstlich erzeugten Testdatensatzes erläutert.

set.seed(123)

Alter <- c(24, 22, 28, 25, 41, 39, 35, 40, 62, 57, 60, 55)

Einkommen <- c(20000, 22000, 25000, 24000, 55000, 65000, 75000, 60000, 30000, 34000, 30000, 34000)

Daten <- as.data.frame(cbind(Alter, Einkommen))

Zunächst wird ein Testdatensatz mit den Variablen „Alter“ und „Einkommen“ erzeugt, der 12 Fälle enthält. Als Schritt des „Data preprocessing“ müssen zunächst beide Variablen standardisiert werden, da ansonsten die Variable „Alter“ die Clusterbildung zu stark beeinflusst.

1 2	Daten$Alter <- scale(Daten$Alter) Daten$Einkommen <- scale(Daten$Einkommen)

Das Ganze geplottet:

plot(Daten$Alter, Daten$Einkommen, col = "blue", pch = 19,

xlab = "Alter (scaled)",

ylab = "Einkommen (scaled)",

main = "Alter vs. Einkommen (scaled)")

Wie bereits eingangs erwähnt müssen Cluster innerhalb möglichst homogen und zu Objekten anderer Cluster möglichst heterogen sein. Ein Maß für die Homogenität die „Within Cluster Sums of Squares“ (WSS), ein Maß für die Heterogenität „Between Cluster Sums of Squares“ (BSS).

Diese sind beispielsweise für eine 3-Cluster-Lösung wie folgt:

KmeansObj <- kmeans(Daten, 3, nstart = 20)

KmeansObj$withinss # Within Cluster Sums of Squares (WSS)

KmeansObj$totss # Between Cluster Sums of Suqares (BSS)

> KmeansObj$withinss # Within Cluster Sums of Squares (WSS)

[1] 0.7056937 0.1281607 0.1792432

> KmeansObj$totss # Between Cluster Sums of Suqares (BSS)

[1] 22

Sollte man die Anzahl der Cluster nicht bereits kennen oder sind diese extern nicht vorgegeben, dann bietet es sich an, anhand des Verhältnisses von WSS und BSS die „optimale“ Clusteranzahl zu berechnen. Dafür wird zunächst ein leerer Vektor initialisiert, dessen Werte nachfolgend über die Schleife mit dem Verhältnis von WSS und WSS gefüllt werden. Dies lässt sich anschließend per „Screeplot“ visualisieren.

ratio <- vector()

for (k in 1:6) {

KMeansObj <- kmeans(Daten, k, nstart = 20)

ratio[k] <- KMeansObj$tot.withinss / KMeansObj$totss

}

plot(ratio, type = "b",

xlab = "Anzahl der Cluster",

ylab = "Ratio WSS/BSS",

main = "Screeplot für verschiedene Clusterlösungen",

col = "blue", pch = 19)

Die „optimale“ Anzahl der Cluster zählt sich am Knick der Linie ablesen (auch Ellbow-Kriterium genannt). Alternativ kann man sich an dem Richtwert von 0.2 orientieren. Unterschreitet das Verhältnis von WSS und BSS diesen Wert, so hat man die beste Lösung gefunden. In diesem Beispiel ist sehr deutlich, dass eine 3-Cluster-Lösung am besten ist.

KmeansObj <- kmeans(Daten, centers = 3, nstart = 20)

plot(Daten$Alter, Daten$Einkommen, col = KmeansObj$cluster, pch = 19, cex = 4,

xlab = "Alter (scaled)",

ylab = "Einkommen (scaled)",

main = "3-Cluster-Lösung")

points(KmeansObj$centers, col = 1:3, pch = 3, cex = 5, lwd = 5)

Fazit: Mit K-Means Clustering lassen sich schnell und einfach Muster in Datensätzen erkennen, die, gerade wenn mehr als zwei Variablen geclustert werden, sonst verborgen blieben. K-Means ist allerdings anfällig gegenüber Ausreißern, da Ausreißer gerne als separate Cluster betrachtet werden. Ebenfalls problematisch sind Cluster, deren Struktur nicht kugelförmig ist. Dies ist vor der Durchführung der Clusteranalyse mittels explorativer Datenanalyse zu überprüfen.

R als Tool im Process Mining

March 15, 2017/2 Comments/in Business Analytics, Business Intelligence, Data Mining, Data Science, Data Science Hack, Industrie 4.0, Main Category, Process Mining, R Statistics, Statistics /by Erwin Enzenhofer

Die Open Source Sprache R ermöglicht eine Vielzahl von Analysemöglichkeiten, die von einer einfachen beschreibenden Darstellung eines Prozesses bis zur umfassenden statistischen Analyse reicht. Dabei können Daten aus einem Manufacturing Execution System, kurz MES, als Basis der Prozessanalyse herangezogen werden. R ist ein Open Source Programm, welches sich für die Lösung von statischen Aufgaben im Bereich der Prozessoptimierung sehr gut eignet, erfordert jedoch auf Grund des Bedienungskonzepts als Scriptsprache, grundlegende Kenntnisse der Programmierung. Aber auch eine interaktive Bedienung lässt sich mit einer Einbindung der Statistikfunktionen in ein Dashboard erreichen. Damit können entsprechend den Anforderungen, automatisierte Analysen ohne Programmierkenntnisse realisiert werden.

Der Prozess als Spagetti Diagramm

Um einen Überblick zu erhalten, wird der Prozess in einem „process value flowchart“, ähnlich einem Spagetti‐ Diagramm dargestellt und je nach Anforderung mit Angaben zu den Key Performance Indicators ergänzt. Im konkreten Fall werden die absolute Anzahl und der relative Anteil der bearbeiteten Teile angegeben. Werden Teile wie nachfolgend dargestellt, aufgrund von festgestellten Mängel bei der Qualitätskontrolle automatisiert ausgeschleust, können darüber Kennzahlen für den Ausschuss ermittelt werden.

Der Prozess in Tabellen und Diagrammen

Im folgenden Chart sind grundlegende Angaben zu den ausgeführten Prozessschritten, sowie deren Varianten dargestellt. Die Statistikansicht bietet eine Übersicht zu den Fällen, den sogenannte „Cases“, sowie zur Dauer und Taktzeit der einzelnen Aktivitäten. Dabei handelt es sich um eine Fertigungsline mit hohem Automatisierungsgrad, bei der jeder Fertigungsschritt im MES dokumentiert wird. Die Tabelle enthält statistische Angaben zur Zykluszeit, sowie der Prozessdauer zu den einzelnen Aktivitäten. In diesem Fall waren keine Timestamps für das Ende der Aktivität vorhanden, somit konnte die Prozessdauer nicht berechnet werden.

Die Anwendung von Six Sigma Tools

R verfügt über eine umfangreiche Sammlung von Bibliotheken zur Datendarstellung, sowie der Prozessanalyse. Darin sind auch Tools aus Six Sigma enthalten, die für die weitere Analyse der Prozesse eingesetzt werden können. In den folgenden Darstellungen wird die Möglichkeit aufgezeigt, zwei Produktionszeiträume, welche über eine einfache Datumseingabe im Dashboard abgegrenzt werden, gegenüber zu stellen. Dabei handelt es sich um die Ausbringung der Fertigung in Stundenwerten, die für jeden Prozessschritt errechnet wird. Das xbar und r Chart findet im Bereich der Qualitätssicherung häufig Anwendung zur ersten Beurteilung des Prozessoutputs.

Zwei weitere Six Sigma typische Kennzahlen zur Beurteilung der Prozessfähigkeit sind der Cp und Cpk Wert und deren Ermittlung ein Bestandteil der R Bibliotheken ist. Bei der Berechnung wird von einer Normalverteilung der Daten ausgegangen, wobei das Ergebnis aus der Überprüfung dieser Annahme im Chart durch Zahlen, als auch grafisch dargestellt wird.

Von Interesse ist auch die Antwort auf die Frage, welchem Trend folgt der Prozess? Bereits aus der Darstellung der beiden Produktionszeiträume im Box‐Whiskers‐Plot könnte man anhand der Mediane auf einen Trend zu einer Verschlechterung der Ausbringung schließen, den der Interquartilsabstand nicht widerspiegelt. Eine weitere Absicherung einer Aussage über den Trend, kann über einen statistischen Vergleichs der Mittelwerte erfolgen.

Der Modellvergleich

Besteht die Anforderung einer direkten Gegenüberstellung des geplanten, mit dem vorgefundenen, sogenannten „Discovered Model“, ist aufgrund der Komplexität beim Modellvergleich, dieser in R mit hohem Programmieraufwand verbunden. Besser geeignet sind dafür spezielle Process Miningtools. Diese ermöglichen den direkten Vergleich und unterstützen bei der Analyse der Ursachen zu den dargestellten Abweichungen. Bei Produktionsprozessen handelt es sich meist um sogenannte „Milestone Events“, die bei jedem Fertigungsschritt durch das MES dokumentiert werden und eine einfache Modellierung des Target Process ermöglichen. Weiterführende Analysen der Prozessdaten in R sind durch einen direkten Zugriff über ein API realisierbar oder es wurde vollständig integriert. Damit eröffnen sich wiederum die umfangreichen Möglichkeiten bei der statistischen Prozessanalyse, sowie der Einsatz von Six Sigma Tools aus dem Qualitätsmanagement. Die Analyse kann durch eine, den Kundenanforderungen entsprechende Darstellung in einem Dashboard vereinfacht werden, ermöglicht somit eine zeitnahe, weitgehend automatisierte Prozessanalyse auf Basis der Produktionsdaten.

Resümee

Process Mining in R ermöglicht zeitnahe Ergebnisse, die bis zur automatisierten Analyse in Echtzeit reicht. Der Einsatz beschleunigt erheblich das Process Controlling und hilft den Ressourceneinsatz bei der Datenerhebung, sowie deren Analyse zu reduzieren. Es kann als stand‐alone Lösung zur Untersuchung des „Discovered Process“ oder als Erweiterung für nachfolgende statistische Analysen eingesetzt werden. Als stand‐alone Lösung eignet es sich für Prozesse mit geringer Komplexität, wie in der automatisierten Fertigung. Besteht eine hohe Diversifikation oder sollen standortübergreifende Prozessanalysen durchgeführt werden, übersteigt der Ressourcenaufwand rasch die Kosten für den Einsatz einer Enterprise Software, von denen mittlerweile einige angeboten werden.

Einstieg in das Maschinelle Lernen mit Python(x,y)

January 30, 2017/1 Comment/in Artificial Intelligence, Business Analytics, Business Intelligence, Data Mining, Data Science, Data Science Hack, Data Warehousing, Machine Learning, Predictive Analytics, Python, Statistics /by Christoph Gresch

Python(x,y) ist eine Python-Distribution, die speziell für wissenschaftliche Arbeiten entwickelt wurde. Es umfasst neben der Programmiersprache auch die Entwicklungsumgebung Spyder und eine Reihe integrierter Python-Bibliotheken. Mithilfe von Python(x,y) kann eine Vielzahl von Interessensbereichen bearbeitet werden. Dazu zählen unter anderem Bildverarbeitung oder auch das maschinelle Lernen. Das All-in-One-Setup für Python(x,y) ist für alle gängigen Betriebssysteme online erhältlich. Read more

R Data Frames meistern mit dplyr – Teil 2

January 24, 2017/1 Comment/in Business Analytics, Business Intelligence, Data Migration, Data Mining, Data Science, Data Science Hack, Data Warehousing, Database, R Statistics, Statistics, Tool Introduction /by Dr. Uwe Schneider

Dieser Artikel ist Teil 2 von 2 aus der Artikelserie R Data Frames meistern mit dplyr.

Noch mehr Datenbank-Features

Im ersten Teil dieser Artikel-Serie habe ich die Parallelen zwischen Data Frames in R und Relationen in SQL herausgearbeitet und gezeigt, wie das Paket dplyr eine Reihe von SQL-analogen Operationen auf Data Frames standardisiert und optimiert. In diesem Teil möchte ich nun drei weitere Analogien aufzeigen. Es handelt sich um die

Window Functions in dplyr als Entsprechung zu analytischen Funktionen in SQL,
Joins zwischen Data Frames als Pendant zu Tabellen-Joins
Delegation von Data Frame-Operationen zu einer bestehenden SQL-Datenbank

Window Functions

Im letzten Teil habe ich gezeigt, wie durch die Kombination von group_by() und summarise() im Handumdrehen Aggregate entstehen. Das Verb group_by() schafft dabei, wie der Name schon sagt, eine Gruppierung der Zeilen des Data Frame anhand benannter Schlüssel, die oft ordinaler oder kategorialer Natur sind (z.B. Datum, Produkt oder Mitarbeiter).

Ersetzt man die Aggregation mit summarise() durch die Funktion mutate(), um neue Spalten zu bilden, so ist der Effekt des group_by() weiterhin nutzbar, erzeugt aber „Windows“, also Gruppen von Datensätzen des Data Frames mit gleichen Werten der Gruppierungskriterien. Auf diesen Gruppen können nun mittels mutate() beliebige R-Funktionen angewendet werden. Das Ergebnis ist im Gegensatz zu summarise() keine Verdichtung auf einen Datensatz pro Gruppe, sondern eine Erweiterung jeder einzelnen Zeile um neue Werte. Das soll folgendes Beispiel verdeutlichen:

library(dplyr)

set.seed(42)

df <- data.frame(id = 1:20,

a=sample(c("Hund","Katze","Maus","Tiger"),20,replace=T),

b=sample(1:10,20, replace = T))

id a b

1 1 Maus 7

2 2 Hund 3

3 3 Katze 3

4 4 Maus 4

5 5 Tiger 10

6 6 Maus 10

7 7 Hund 8

8 8 Hund 8

9 9 Hund 6

10 10 Katze 1

11 11 Maus 7

12 12 Hund 9

13 13 Hund 8

14 14 Tiger 5

15 15 Tiger 6

16 16 Maus 6

17 17 Katze 1

18 18 Maus 4

19 19 Maus 7

20 20 Maus 9

df %>%

group_by(a) %>%

mutate(r = row_number(), # aus dplyr

n_memb = n(), # aus dplyr

n_dist = n_distinct(b), # aus dplyr

ra=rank(desc(b)), # aus base und dplyr

last_b = lag(b), # aus dplyr

next_b = lead(b), # aus dplyr

mb = mean(b), # aus base

cs = cumsum(b) ) # aus base

Source: local data frame [20 x 11]

Groups: a [4]

id a b r n_memb n_dist ra last_b next_b mb cs

1 1 Maus 7 1 8 5 4.0 NA 4 6.750000 7

2 2 Hund 3 1 6 4 6.0 NA 8 7.000000 3

3 3 Katze 3 1 3 2 1.0 NA 1 1.666667 3

4 4 Maus 4 2 8 5 7.5 7 10 6.750000 11

5 5 Tiger 10 1 3 3 1.0 NA 5 7.000000 10

6 6 Maus 10 3 8 5 1.0 4 7 6.750000 21

7 7 Hund 8 2 6 4 3.0 3 8 7.000000 11

8 8 Hund 8 3 6 4 3.0 8 6 7.000000 19

9 9 Hund 6 4 6 4 5.0 8 9 7.000000 25

10 10 Katze 1 2 3 2 2.5 3 1 1.666667 4

11 11 Maus 7 4 8 5 4.0 10 6 6.750000 28

12 12 Hund 9 5 6 4 1.0 6 8 7.000000 34

13 13 Hund 8 6 6 4 3.0 9 NA 7.000000 42

14 14 Tiger 5 2 3 3 3.0 10 6 7.000000 15

15 15 Tiger 6 3 3 3 2.0 5 NA 7.000000 21

16 16 Maus 6 5 8 5 6.0 7 4 6.750000 34

17 17 Katze 1 3 3 2 2.5 1 NA 1.666667 5

18 18 Maus 4 6 8 5 7.5 6 7 6.750000 38

19 19 Maus 7 7 8 5 4.0 4 9 6.750000 45

20 20 Maus 9 8 8 5 2.0 7 NA 6.750000 54

Das group_by() unterteilt den Data Frame nach den 4 gleichen Werten von a. Innerhalb dieser Gruppen berechnen die beispielsweise eingesetzten Funktionen

row_number(): Die laufende Nummer in dieser Gruppe
n(): Die Gesamtgröße dieser Gruppe
n_distinct(b): Die Anzahl verschiedener Werte von b innerhalb der Gruppe
rank(desc(b)): Den Rang innerhalb der selben Gruppe, absteigend nach b geordnet
lag(b): Den Wert von b der vorherigen Zeile innerhalb derselben Gruppe
lead(b): Analog den Wert von b der folgenden Zeile innerhalb derselben Gruppe
mean(b): Den Mittelwert von b innerhalb der Gruppe
cumsum(b): Die kumulierte Summe der b-Werte innerhalb der Gruppe.

Wichtig ist hierbei, dass die Anwendung dieser Funktionen nicht dazu führt, dass die ursprüngliche Reihenfolge der Datensätze im Data Frame geändert wird. Hier erweist sich ein wesentlicher Unterschied zwischen Data Frames und Datenbank-Relationen von Vorteil: Die Reihenfolge von Datensätzen in Data Frames ist stabil und definiert. Sie resultiert aus der Abfolge der Elemente auf den Vektoren, die die Data Frames bilden. Im Gegensatz dazu haben Tabellen und Views keine Reihenfolge, auf die man sich beim SELECT verlassen kann. Nur mit der ORDER BY-Klausel über eindeutige Schlüsselwerte erreicht man eine definierte, stabile Reihenfolge der resultierenden Datensätze.

Die Wirkungsweise von Window Functions wird noch besser verständlich, wenn in obiger Abfrage das group_by(a) entfernt wird. Dann wirken alle genannten Funktionen auf der einzigen Gruppe, die existiert, nämlich dem gesamten Data Frame:

df %>%

mutate(r = row_number(), # aus dplyr

n_memb = n(), # aus dplyr

n_dist = n_distinct(b), # aus dplyr

ra=rank(desc(b)), # aus base und dplyr

last_b = lag(b), # aus dplyr

next_b = lead(b), # aus dplyr

mb = mean(b), # aus base

cs = cumsum(b) ) # aus base

id a b r n_memb n_dist ra last_b next_b mb cs

1 1 Maus 7 1 20 9 9.0 NA 3 6.1 7

2 2 Hund 3 2 20 9 17.5 7 3 6.1 10

3 3 Katze 3 3 20 9 17.5 3 4 6.1 13

4 4 Maus 4 4 20 9 15.5 3 10 6.1 17

5 5 Tiger 10 5 20 9 1.5 4 10 6.1 27

6 6 Maus 10 6 20 9 1.5 10 8 6.1 37

7 7 Hund 8 7 20 9 6.0 10 8 6.1 45

8 8 Hund 8 8 20 9 6.0 8 6 6.1 53

9 9 Hund 6 9 20 9 12.0 8 1 6.1 59

10 10 Katze 1 10 20 9 19.5 6 7 6.1 60

11 11 Maus 7 11 20 9 9.0 1 9 6.1 67

12 12 Hund 9 12 20 9 3.5 7 8 6.1 76

13 13 Hund 8 13 20 9 6.0 9 5 6.1 84

14 14 Tiger 5 14 20 9 14.0 8 6 6.1 89

15 15 Tiger 6 15 20 9 12.0 5 6 6.1 95

16 16 Maus 6 16 20 9 12.0 6 1 6.1 101

17 17 Katze 1 17 20 9 19.5 6 4 6.1 102

18 18 Maus 4 18 20 9 15.5 1 7 6.1 106

19 19 Maus 7 19 20 9 9.0 4 9 6.1 113

20 20 Maus 9 20 20 9 3.5 7 NA 6.1 122

Anwendbar sind hierbei sämtliche Funktionen, die auf Vektoren wirken. Diese müssen also wie in unserem Beispiel nicht unbedingt aus dplyr stammen. Allerdings komplettiert das Package die Menge der sinnvoll anwendbaren Funktionen um einige wichtige Elemente wie cumany() oder n_distinct().

Data Frames Hand in Hand…

In relationalen Datenbanken wird häufig angestrebt, das Datenmodell zu normalisieren. Dadurch bekommt man die negativen Folgen von Datenredundanz, wie Inkonsistenzen bei Datenmanipulationen und unnötig große Datenvolumina, in den Griff. Dies geschieht unter anderem dadurch, dass tabellarische Datenbestände aufgetrennt werden Stammdaten- und Faktentabellen. Letztere beziehen sich über Fremdschlüsselspalten auf die Primärschlüssel der Stammdatentabellen. Durch Joins, also Abfragen über mehrere Tabellen und Ausnutzen der Fremdschlüsselbeziehungen, werden die normalisierten Tabellen wieder zu einem fachlich kompletten Resultat denormalisiert.

In den Data Frames von R trifft man dieses Modellierungsmuster aus verschiedenen Gründen weit seltener an als in RDBMS. Dennoch gibt es neben der Normalisierung/Denormalisierung andere Fragestellungen, die sich gut durch Joins beantworten lassen. Neben der Zusammenführung von Beobachtungen unterschiedlicher Quellen anhand charakteristischer Schlüssel sind dies bestimmte Mengenoperationen wie Schnitt- und Differenzmengenbildung.

Die traditionelle R-Funktion für den Join zweier Data Frames lautet merge(). dplyr erweitert den Funktionsumfang dieser Funktion und sorgt für sprechendere Funktionsnamen und Konsistenz mit den anderen Operationen.

Hier ein synthetisches Beispiel:

products <- data.frame(

id = 1:5,

name = c("Desktop", "Laptop", "Maus", "Tablet", "Smartphone"),

preis = c(500, 700, 10, 300, 500)

)

set.seed(1)

(salesfacts <- data.frame(

prod_id = sample(1:5,size = 8,replace = T),

date = as.Date('2017-01-01') + sample(1:5,size = 8,replace = T)

) )

prod_id date

1 2 2017-01-05

2 2 2017-01-02

3 3 2017-01-03

4 5 2017-01-02

5 2 2017-01-05

6 5 2017-01-03

7 5 2017-01-05

8 4 2017-01-04

Nun gilt es, die Verkäufe aus dem Data Frame sales mit den Produkten in products zusammenzuführen und auf Basis von Produkten Bilanzen zu erstellen. Diese Denormalisierung geschieht durch das Verb inner_join() auf zweierlei Art und Weise:

salesfacts %>%

inner_join(products, by = c("prod_id" = "id"))

prod_id date name preis

1 2 2017-01-05 Laptop 700

2 2 2017-01-02 Laptop 700

3 3 2017-01-03 Maus 10

4 5 2017-01-02 Smartphone 500

5 2 2017-01-05 Laptop 700

6 5 2017-01-03 Smartphone 500

7 5 2017-01-05 Smartphone 500

8 4 2017-01-04 Tablet 300

products %>%

inner_join(salesfacts, by = c("id" = "prod_id"))

id name preis date

1 2 Laptop 700 2017-01-05

2 2 Laptop 700 2017-01-02

3 2 Laptop 700 2017-01-05

4 3 Maus 10 2017-01-03

5 4 Tablet 300 2017-01-04

6 5 Smartphone 500 2017-01-02

7 5 Smartphone 500 2017-01-03

8 5 Smartphone 500 2017-01-05

Die Ergebnisse sind bis auf die Reihenfolge der Spalten und der Zeilen identisch. Außerdem ist im einen Fall der gemeinsame Schlüssel der Produkt-Id als prod_id, im anderen Fall als id enthalten. dplyr entfernt also die Spalten-Duplikate der Join-Bedingungen. Letzere wird bei Bedarf im by-Argument der Join-Funktion angegeben. R-Experten erkennen hier einen „Named Vector“, also einen Vektor, bei dem jedes Element einen Namen hat. Diese Syntax verwendet dplyr, um elegant die äquivalenten Spalten zu kennzeichnen. Wird das Argument by weggelassen, so verwendet dplyr im Sinne eines „Natural Join“ automatisch alle Spalten, deren Namen in beiden Data Frames vorkommen.

Natürlich können wir dieses Beispiel mit den anderen Verben erweitern, um z.B. eine Umsatzbilanz pro Produkt zu erreichen:

salesfacts %>%

inner_join(products, by = c("prod_id" = "id")) %>%

group_by(prod_id) %>%

summarise(n_verk = n(), sum_preis = sum(preis), letzt_dat = max(date))

# A tibble: 4 × 4

prod_id n_verk sum_preis letzt_dat

1 2 3 2100 2017-01-05

2 3 1 10 2017-01-03

3 4 1 300 2017-01-04

4 5 3 1500 2017-01-05

dplyr bringt insgesamt 6 verschiedene Join-Funktionen mit: Neben dem bereits verwendeten Inner Join gibt es die linksseitigen und rechtsseitigen Outer Joins und den Full Join. Diese entsprechen genau der Funktionalität von SQL-Datenbanken. Daneben gibt es die Funktion semi_join(), die in SQL etwa folgendermaßen ausgedrückt würde:

SELECT ...

FROM a

WHERE EXISTS (SELECT * FROM b WHERE b.a_id = a.id)

Das Gegenteil, also ein NOT EXISTS, realisiert die sechste Join-Funktion: anti_join(). Im folgenden Beispiel sollen alle Produkte ausgegeben werden, die noch nie verkauft wurden:

products %>% anti_join(salesfacts,c("id" = "prod_id"))

id name preis

1 1 Desktop 500

… und in der Datenbank

Wir schon mehrfach betont, hat dplyr eine Reihe von Analogien zu SQL-Operationen auf relationalen Datenbanken. R Data Frames entsprechen Tabellen und Views und die dplyr-Operationen den Bausteinen von SELECT-Statements. Daraus ergibt sich die Möglichkeit, dplyr-Funktionen ohne viel Zutun auf eine bestehende Datenbank und deren Relationen zu deligieren.

Mir fallen folgende Szenarien ein, wo dies sinnvoll erscheint:

Die zu verarbeitende Datenmenge ist zu groß für das Memory des Rechners, auf dem R läuft.
Die interessierenden Daten liegen bereits als Tabellen und Views auf einer Datenbank vor.
Die Datenbank hat Features, wie z.B. Parallelverarbeitung oder Bitmap Indexe, die R nicht hat.

In der aktuellen Version 0.5.0 kann dplyr nativ vier Datenbank-Backends ansprechen: SQLite, MySQL, PostgreSQL und Google BigQuery. Ich vermute, unter der Leserschaft des Data Science Blogs dürfte MySQL (oder der Fork MariaDB) die weiteste Verbreitung haben, weshalb ich die folgenden Beispiele darauf zeige. Allerdings muss man beachten, dass MySQL keine Window Funktionen kennt, was sich 1:1 auf die Funktionalität von dplyr auswirkt.

Im folgenden möchte ich zeigen, wie dplyr sich gegen eine bestehende MySQL-Datenbank verbindet und danach einen bestehenden R Data Frame in eine neue Datenbanktabelle wegspeichert:

mysql_db <- src_mysql(host = "localhost", user = "testuser",

password = "********", dbname = "test")

library(ggplot2)

str(diamonds)

Classes ‘tbl_df’, ‘tbl’ and 'data.frame': 53940 obs. of 10 variables:

$ carat : num 0.23 0.21 0.23 0.29 0.31 0.24 0.24 0.26 0.22 0.23 ...

$ cut : chr "Ideal" "Premium" "Good" "Premium" ...

$ color : chr "E" "E" "E" "I" ...

$ clarity: chr "SI2" "SI1" "VS1" "VS2" ...

$ depth : num 61.5 59.8 56.9 62.4 63.3 62.8 62.3 61.9 65.1 59.4 ...

$ table : num 55 61 65 58 58 57 57 55 61 61 ...

$ price : int 326 326 327 334 335 336 336 337 337 338 ...

$ x : num 3.95 3.89 4.05 4.2 4.34 3.94 3.95 4.07 3.87 4 ...

$ y : num 3.98 3.84 4.07 4.23 4.35 3.96 3.98 4.11 3.78 4.05 ...

$ z : num 2.43 2.31 2.31 2.63 2.75 2.48 2.47 2.53 2.49 2.39 ...

diamonds %>% mutate(cut = as.character(cut),

color = as.character(color),

clarity = as.character(clarity)) -> diamonds

diamonds_mysql <- copy_to(mysql_db, diamonds, name="diamonds",

temporary = FALSE, indexes = list(

c("cut", "color", "clarity"), "carat", "price"))

diamonds_mysql %>% summarise(count = n())

Source: query [?? x 1]

Database: mysql 5.5.54-0ubuntu0.14.04.1 [testuser@localhost:/test]

count

1 53940

Die erste Anweisung verbindet R mit einer bestehenden MySQL-Datenbank. Danach lade ich den Data Frame diamonds aus dem Paket ggplot2. Mit str() wird deutlich, dass drei darin enthaltene Variablen vom Typ Factor sind. Damit dplyr damit arbeiten kann, werden sie mit mutate() in Character-Vektoren gewandelt. Dann erzeugt die Funktion copy_to() auf der MySQL-Datenbank eine leere Tabelle namens diamonds, in die die Datensätze kopiert werden. Danach erhält die Tabelle noch drei Indexe (von dem der erste aus drei Segmenten besteht), und zum Schluß führt dplyr noch ein ANALYSE der Tabelle durch, um die Werteverteilungen auf den Spalten für kostenbasierte Optimierung zu bestimmen.

Meistens aber wird bereits eine bestehende Datenbanktabelle die interessierenden Daten enthalten. In diesem Fall lautet die Funktion zum Erstellen des Delegats tbl():

diamonds_mysql2 <- tbl(mysql_db,"diamonds")

identical(diamonds_mysql,diamonds_mysql2)

[1] TRUE

Die Rückgabewerte von copy_to() und von tbl() sind natürlich keine reinrassigen Data Frames, sondern Objekte, auf die die Operationen von dplyr wirken können, indem sie auf die Datenbank deligiert werden. Im folgenden Beispiel sollen alle Diamanten, die ein Gewicht von mindestens 1 Karat haben, pro Cut, Color und Clarity nach Anzahl und mittlerem Preis bilanziert werden:

bilanz <- diamonds_mysql2 %>%

filter(carat >= 1) %>%

group_by(cut,color,clarity) %>%

summarise(count = n(), mean_price = mean(price))

bilanz

Source: query [?? x 5]

Database: mysql 5.5.54-0ubuntu0.14.04.1 [testuser@localhost:/test]

Groups: cut, color

cut color clarity count mean_price

1 Fair D I1 3 9013.667

2 Fair D SI1 26 6398.192

3 Fair D SI2 29 6138.552

4 Fair D VS1 1 7083.000

5 Fair D VS2 7 8553.429

6 Fair D VVS1 1 10752.000

7 Fair D VVS2 2 9639.000

8 Fair E I1 5 2469.800

9 Fair E SI1 28 6407.464

10 Fair E SI2 45 5627.489

# ... with more rows

explain(bilanz)

SELECT `cut`, `color`, `clarity`, count(*) AS `count`, AVG(`price`) AS `mean_price`

FROM (SELECT *

FROM `diamonds`

WHERE (`carat` >= 1.0)) `cttxnwlelz`

GROUP BY `cut`, `color`, `clarity`

id select_type table type possible_keys key key_len ref rows

1 1 PRIMARY ALL 19060

2 2 DERIVED diamonds ALL diamonds_carat 50681

Extra

1 Using temporary; Using filesort

2 Using where

Die Definition der Variablen bilanz geschieht dabei komplett ohne Interaktion mit der Datenbank. Erst beim Anzeigen von Daten wird das notwendige SQL ermittelt und auf der DB ausgeführt. Die ersten 10 resultierenden Datensätze werden angezeigt. Mittels der mächtigen Funktion explain() erhalten wir das erzeugte SQL-Kommando und sogar den Ausführungsplan auf der Datenbank. SQL-Kundige werden erkennen, dass die verketteten dplyr-Operationen in verschachtelte SELECT-Statements umgesetzt werden.

Zu guter Letzt sollen aber meistens die Ergebnisse der dplyr-Operationen irgendwie gesichert werden. Hier hat der Benutzer die Wahl, ob die Daten auf der Datenbank in einer neuen Tabelle gespeichert werden sollen oder ob sie komplett nach R transferiert werden sollen. Dies erfolgt mit den Funktionen compute() bzw. collect():

compute(bilanz, name = "t_bilanz", temporary = F)

df <- collect(bilanz)

str(df)

Classes ‘grouped_df’, ‘tbl_df’, ‘tbl’ and 'data.frame': 265 obs. of 5 variables:

$ cut : chr "Fair" "Fair" "Fair" "Fair" ...

$ color : chr "D" "D" "D" "D" ...

$ clarity : chr "I1" "SI1" "SI2" "VS1" ...

$ count : num 3 26 29 1 7 1 2 5 28 45 ...

$ mean_price: num 9014 6398 6139 7083 8553 ...

...

Durch diese beiden Operationen wurde eine neue Datenbanktabelle „t_bilanz“ erzeugt und danach der Inhalt der Bilanz als Data Frame zurück in den R-Interpreter geholt. Damit schließt sich der Kreis.

Fazit

Mit dem Paket dplyr von Hadley Wickham wird die Arbeit mit R Data Frames auf eine neue Ebene gehoben. Die Operationen sind konsistent, vollständig und performant. Durch den Verkettungs-Operator %>% erhalten sie auch bei hoher Komplexität eine intuitive Syntax. Viele Aspekte der Funktionalität lehnen sich an Relationale Datenbanken an, sodass Analysten mit SQL-Kenntnissen rasch viele Operationen auf R Data Frames übertragen können.

Zurück zu R Data Frames meistern mit dplyr – Teil 1.

Numerical Python – Einführung in wissenschaftliches Rechnen mit NumPy

January 23, 2017/2 Comments/in Data Science, Python, Python, Statistics /by Benjamin Aunkofer

NumPy steht für Numerical Python und ist eines der bekanntesten Pakete für alle Python-Programmierer mit wissenschaftlichen Hintergrund. Von persönlichen Kontakten erfuhr ich, dass NumPy heute in der Astrophysik fast genauso verwendet wird wie auch von sogenannten Quants im Investment-Banking. Das NumPy-Paket ist sicherlich ein Grundstein des Erfolges für Python in der Wissenschaft und für den häufigen Einsatz für die Implementierung von Algorihtmen des maschinellen Lernens in Python.

Die zentrale Datenstruktur in NumPy ist das mehrdimensionale Array. Dieses n-dimensionale Array (ndarray) ist eine sehr mächtige Datenstruktur und verwende ich beispielsweise in meinem Artikel über den k-Nächste-Nachbarn-Algorithmus. Die Besonderheit des NumPy-Arrays ist, dass es ein mehrdimensionaler Container für homogene Daten ist. Ein Datentyp gilt also für das gesamte Array, nicht nur für bestimmte Zeilen oder Spalten!

1 2	import numpy as np

Attribution Models

A Business and Statistical Case

INTRODUCTION

1. A/B-Tests

2. Attribution models

2.1 General Issues

2.2 Today’s main practices

3. Data-Driven Attribution Models in practice

3.1 Issues

3.2 Main models

A) Logistic Regression and Classification models

B) Shapley Value

Theory

Advantages

B) Markov Chains

3.3 To go further: Tackling selection biases with Quasi-Experiments

4. First Steps towards a Practical Implementation

Identify key points of interests

Analyse the typical customer journey

Summary

The Relevance of Domain Knowledge

Discussion

Countries of Happiness – the Full Article

1. Schritt: Auswählen des Attributes mit dem höchsten Informationsgewinn

1.1 Gesamt-Entropie des Datensatzes berechnen

1.2 Berechnung der Informationsgewinne aller Attribute

2. Schritt: Anlegen des Wurzel-Knotens

3. Schritt: Rekursive Wiederholung (!!!)

3.1 Erster Rekursionsschritt

3.1 Zweiter Rekursionsschritt

3.1 Dritter Rekursionsschritt

Ergebnis: Der Entscheidungsbaum

Die Daten als Download

Die Entwicklungsumgebung

Datenbereinigung

Datenvisualisierung mit Pandas

Vorbereitung der Regressionsanalyse

Lineare Regression mit Scikit-Learn

Nicht-lineare Regression mit Scikit-Learn

Weitere Anmerkungen

Der Prozess als Spagetti Diagramm

Der Prozess in Tabellen und Diagrammen

Die Anwendung von Six Sigma Tools

Der Modellvergleich

Resümee

Noch mehr Datenbank-Features

Window Functions

Data Frames Hand in Hand…

… und in der Datenbank

Fazit

Interesting links

Pages

Categories

Archive