Advanced Attribution: Markov Chain

Why rule-based attribution is not good enough?

Continuing our discussion about attribution modeling, We’ve already convered basic rule-based attribution, which can by easily produced in excel. One major limitation of rule-based attribution is that it ignores customers that don’t convert. This means that we throw out the majority of data given that the ratio of non-converting paths tends to be higher than that of converting paths, which in return can produce biased result.

In order to take advantage of non-converting paths, Two advanced attribution methodologies are introduced. I will cover Markov Chain in this post, and the data-mining approach in the next post.

What is markov chain?

The mathmatical definition of Markov Chain is “a stochastic model describing a sequence of possible events in which the probability of each event depends only on the state attained in the previous event.”

Markov chain describes a probabilistic relationship between states. Here is a sample example:

Given we have two channels, a and b; and two possible outcomes, convert or churn. A series of probabilities are given as follows:

P(Customer start(first interact) with A channel): 0.2

P(Customer start(first interact) with B channel): 0.8

P(Customer interact with A|Given previously interact with B): 0.3

P(Customer interact with B|Given previously interact with A): 0.4

P(Customer Convert|Given previously interact with A): 0.3

P(Customer Churn|Given previously interact with A): 0.3

P(Customer Convert|Given previously interact with B): 0.25

P(Customer Churn|Given previously interact with B): 0.45

We can constract a matrix out from these relationship.

example_matrix = matrix(c(0,0,0,0,0,0.2,0,0.3,0,0,0.8,0.4,0,0,0,0,0.3,0.25,0,0,0,0.3,0.45,1,1),nrow = 5,ncol = 5)
colnames(example_matrix) <- c("start","a","b","convert","churn")
rownames(example_matrix) <- c("start","a","b","convert","churn")
print(example_matrix)

        start   a   b convert churn
start       0 0.2 0.8    0.00  0.00
a           0 0.0 0.4    0.30  0.30
b           0 0.3 0.0    0.25  0.45
convert     0 0.0 0.0    0.00  1.00
churn       0 0.0 0.0    0.00  1.00

How do we construct a matrix like this? If you take a close look, you will find out that the matrix entry [x,y] annotates the probability of transiting from state x to state y, namingly, the probability of x occurs given y occurs.

This matrix is called transition matrix.

What problem does transition matrix solve?

Say if we are given the transition matrix, what is the probability of conversion given one interaction, given 2 interactions, given n interactions?

Let’s forget about the transition matrix and solve it using conditional probability.

Q: What is the probability of conversion given one interaction? A:since there is no direct connection between start to convert, the probability of conversion after one interaction will be 0.

Q:What is the probablity of conversion given two interaction? A:P(Customer start(first interact) with A channel) x P(Customer Convert|Given previously interact with A) + P(Customer start(first interact) with B channel) x P(Customer Convert|Given previously interact with B): 0.7

0.2 x 0.3 + 0.8 x 0.25 = 0.15 + 0.14 = 0.26

As you can see, it is still pretty managible to compute the first couple interactions, but it will get more cumbersome as number of interaction increases. Also if the number of channels increases, the complexity of computation will also increase. In short, this solution of not scalable.

Let’s look at the althernatives. (For you to understand this, a little bit of brushing off on simple linear algebra calculation may help you a long way:https://www.khanacademy.org/math/linear-algebra)

inital_state <- c(1,0,0,0,0)
inital_state%*%example_matrix

     start   a   b convert churn
[1,]     0 0.2 0.8       0     0

Look at the convert column, We get the same result!

initial_state <- c(1,0,0,0,0)
initial_state%*%example_matrix%*%example_matrix

     start    a    b convert churn
[1,]     0 0.24 0.08    0.26  0.42

Again, the same result is produced.

And there comes the trick when the number of iteration scales up.

logic <- TRUE
initial_state <- c(1,0,0,0,0)
new_state <- initial_state
iter = 0
conv.proba <- new_state[4]
precision <- 0.0001
while(logic) {
  old_state <- new_state
  new_state <- new_state%*%example_matrix
  conv.proba <- conv.proba+new_state[4]
  #As iteraction approach infinity, all customers will eventually migrate to "Churn". therefore it is safe to say that our program approaches equilibrium when probability for churn approaches 1 and other channels approachs 0.
  logic <- c(new_state[1:4] < precision,(abs(new_state[5] - 1) < precision))
  logic <- sum(logic==FALSE) > 0
  iter <- iter + 1
}
print(paste("after",iter,"iterations","the probability of conversion is",round(conv.proba,2)))

[1] "after 11 iterations the probability of conversion is 0.4"

Scaling this program up will give us the probability of conversion for all numbers of interactions, which is the global conversion rate of this program. From the program we just ran, we know that the global conversion rate is 0.45, given there are two channels a,b, and the transition matrix.

Now say, we think that a is not performing and we decide to stop investment in channel a. but we’d like to gauge the impact of this decision before we implement it. What can we do?

We can simplely remove channel a from the transition matrix:

(Here I need to point out that this operation is based on premise that the remaining probabilities will be independent of channel a. which may not be the case in real life scenarios. It is possible that without channel a, a subset of customers who previously interact with channel a gets transmitted to channel b - You will see later in an example in which I used a different assumption).

example_matrix_b <- example_matrix[-2,-2]

Error: object 'example_matrix' not found

Let’s calculate the conversion rate when there is only channel b.

logic <- TRUE
initial_state <- c(1,0,0,0)
new_state <- initial_state
iter = 0
conv.proba.no.a <- new_state[3]
while(logic) {
  old_state <- new_state
  new_state <- new_state%*%example_matrix_b
  conv.proba.no.a <- conv.proba.no.a+new_state[3]
  logic <- c(new_state[1:3] < precision,(abs(new_state[4] - 1) < precision))
  logic <- sum(logic==FALSE) > 0
  iter <- iter + 1
}
print(paste("after",iter,"iterations","the probability of conversion is",round(conv.proba.no.a,2)))

[1] "after 3 iterations the probability of conversion is 0.2"

We think that channel a is not efficient but without channel a the conversion rate is reduced from 0.45 to 0.2! this means that channel a increases the probability of conversion by 200%! it will be a bad decision to scale down investment in channel a.

Let’s do the same thing to channel b.

example_matrix_a <- example_matrix[-3,-3]
print(example_matrix_a)

        start   a convert churn
start       0 0.2     0.0   0.0
a           0 0.0     0.3   0.3
convert     0 0.0     0.0   1.0
churn       0 0.0     0.0   1.0

example_matrix_a[1,4] <- 1 - sum(example_matrix_a[1,]) + example_matrix_a[1,4]
example_matrix_a[2,4] <- 1 - sum(example_matrix_a[2,]) + example_matrix_a[2,4]
example_matrix_a

        start   a convert churn
start       0 0.2     0.0   0.8
a           0 0.0     0.3   0.7
convert     0 0.0     0.0   1.0
churn       0 0.0     0.0   1.0

logic <- TRUE
initial_state <- c(1,0,0,0)
new_state <- initial_state
iter = 0
conv.proba.no.b <- new_state[3]
while(iter<1000) {
  old_state <- new_state
  new_state <- new_state%*%example_matrix_a
  conv.proba.no.b <- conv.proba.no.b+new_state[3]
  #logic <- 1-new_state[4] > 0.0000001 #As number of iteration approach infinity, the probablities in each state will be approaching 0. Here we set the cutting point at 0.0000001 instead of 0 to avoid program running too long or crashing the computer. Essentially I'm saying I will not care the probability loss as long as it is less 0.000001.
  iter <- iter + 1
}
print(paste("after",iter,"iterations","the probability of conversion is",round(conv.proba.no.b,2)))

[1] "after 1000 iterations the probability of conversion is 0.06"

After removing channel b, the conversion rate is reduced from 0.4 to 0.06. Channel b increases the probability of conversion by 667%.

From here we can easily induce attribution model from the results we have got.

a.lift <- conv.proba/conv.proba.no.a
b.lift <- conv.proba/conv.proba.no.b
a.attr <- a.lift/(a.lift+b.lift)
b.attr <- b.lift/(a.lift+b.lift)
print(paste("attributed value of channel a:",round(a.attr,2)))

[1] "attributed value of channel a: 0.23"

print(paste("attributed value of channel b:",round(b.attr,2)))

[1] "attributed value of channel b: 0.77"

Shortcut: R packages

I built this calculation pretty much from scratch with minimum dependencies on existing R packages. The intention of this is solely pedagogical, as I believe only by doing the calculation step-by-step, taking no shortcut could we truely grasp the mechanism of how things work. But it is a loss of efficiency if you do not leverage existing packages that can acomplish the same task faster and better. As far as I know there are two packages: markovchain(https://cran.r-project.org/web/packages/markovchain/index.html) and ChannelAttribution(https://cran.r-project.org/web/packages/ChannelAttribution/ChannelAttribution.pdf). I will briefly show you want to use markovchain. And I will use ChannelAttribution in a seperate post for comparasion purposes.

#Import library 
library(markovchain)
states <- c("start","a","b","convert","churn")
example.mc <- new("markovchain", states = states,transitionMatrix = example_matrix,name = "attribution")

plot(example.mc)

initialState <- c(1,0,0,0,0)
after2Iter <- initialState * (example.mc * example.mc)
print(after2Iter)

     start    a    b convert churn
[1,]     0 0.24 0.08    0.26  0.42

after7Iter <- initialState * (example.mc ^ 7)
print(after7Iter)

     start         a         b   convert     churn
[1,]     0 0.0003456 0.0013824 0.0013248 0.9969472

Production Code on Stimulated Results

Now having grasped how Markov Chain works, let’s apply this method on a stimuated dataset. Aslo I’m going to modularize the code by creating to functions, ConvRateEst and newTransitionMatrix to aid reproduction.

data <- read.csv("C:/Users/rinlin/Desktop/Attribution/Blogpost/attribution_raw.csv",stringsAsFactors = FALSE)

ConvRateEst <- function(TransitionMatrix) {
  iter <- 0
  logic <- TRUE
  precision <- 0.0001
  
  new_state <- matrix(rep(0,length(states(TransitionMatrix))), nrow = 1)
  colnames(new_state) <- states(TransitionMatrix)
  new_state[,"start"] <- 1
  
  ConvRate <- new_state[,"converted"]
  
  
  
  while(logic) {
    
    new_state <- new_state*TransitionMatrix
    ConvRate <- ConvRate + new_state[,"converted"]
    logic <- ((abs(new_state[,"churned"] - 1) < precision))
    logic <- sum(logic==FALSE) > 0
    iter <- iter + 1
  }
  
  print(paste("after",iter,"iterations","the probability of conversion is",round(ConvRate,2)))
  return(ConvRate)
  
}
newTransitionMatrix <- function(mcMatrix,removeChannel) {
  newMatrix <- mcMatrix[-which(states(mcMatrix) == removeChannel),-which(states(mcMatrix) == removeChannel)]
  for(state in rownames(newMatrix)) {
    total <- sum(newMatrix[state,])
    for(state.2 in rownames(newMatrix)) {
      newMatrix[state,state.2] <- newMatrix[state,state.2]/total 
    }
    #newMatrix[state,"churned"] <- 1- sum(newMatrix[state,]) + newMatrix[state,"churned"]
  }
  newMatrix <- new("markovchain",states = rownames(newMatrix),transitionMatrix = newMatrix,name = "attribution")
  return(newMatrix)
}

data.1 <- data
data.1$status <- "start"
data.2 <- data
data.2$status <- "end"
df <- rbind(data.1,data.2)
temp <- df %>%
  group_by(Path.ID,Conversion) %>%
  summarise_at("Time.Stamp",funs(max)) %>%
  ungroup()
end <- temp 
end$Time.Stamp <- end$Time.Stamp + 1
end$Channel <- ifelse(end$Conversion == 1,"converted","churned")
end$status <- "end"
start <- temp
start$Time.Stamp <- 0
start$Channel <- "start"
start$status <- "start"
df <- rbind(df,start,end)
df <- df %>%
  arrange(Path.ID,Time.Stamp,status) %>%
  group_by(Path.ID) %>%
  mutate(new.id = lag(Time.Stamp)) %>%
  ungroup() 
df[is.na(df$new.id),]$new.id <- 0
output <- df %>%
  mutate(id = paste0(Path.ID,new.id)) %>%
  select(id,status,Channel) %>%
  tidyr::spread(status,Channel) %>%
  group_by(start,end) %>%
  summarise_at("id",funs(length)) %>%
  ungroup() %>%
  mutate(count = id) %>%
  select(-id) 
stop()

Error: 

ConvRate <- ConvRateEst(mc.attr)

[1] "after 20 iterations the probability of conversion is 0.37"

ConRateNoSearch <- ConvRateEst(newTransitionMatrix(mc.attr,"Search"))

[1] "after 15 iterations the probability of conversion is 0.19"

ConRateNoDisplay <- ConvRateEst(newTransitionMatrix(mc.attr,"Display"))

[1] "after 13 iterations the probability of conversion is 0.45"

ConRateNoEmail <- ConvRateEst(newTransitionMatrix(mc.attr,"Email"))

[1] "after 17 iterations the probability of conversion is 0.37"

SearchLift <- ConvRate/ConRateNoSearch
DisplayLift <- ConvRate/ConRateNoDisplay
EmailLift <- ConvRate/ConRateNoEmail
data.frame(channel = c("Search","Display","Email"),attribution = c(round(SearchLift/(SearchLift + DisplayLift + EmailLift),2),round(DisplayLift/(SearchLift + DisplayLift + EmailLift),2),round(EmailLift/(SearchLift + DisplayLift + EmailLift),2)))

Other use cases on Markov Chain

As you might have already found out, often time one statistical methods can be applied to multiple disciplines.Another use case of Markov Chain in marketing is in CRM analytics(Which I will be talking about in a different post on migration model). It is also possible to apply markov chain into other fields as long as transitions of states exist. Can you think of any other scenarios where Markov Chain can be applied?