## Measuring user retention using cohort analysis with R

Cohort analysis is super important if you want to know if your service is in fact a leaky bucket despite nice growth of absolute numbers. There’s a good write up on that subject “Cohorts, Retention, Churn, ARPU” by Matt Johnson.

So how to do it using R and how to visualize it. Inspired by examples described in “Retention, Cohorts, and Visualizations” I came up with the following solution.

First, get the data in a suitable format, like this:

```cohort  signed_up  active_m0  active_m1  active_m2
2011-10 12345      10432      8765       6754
2011-11 12345      10432      8765       6754
2011-12 12345      10432      8765       6754```

Cohort here is in “YYYY-MM” format, signed_up is the number of users who have created accounts in the given month, active_m0 – number of users who have been active in the same month as they registered, active_m1 – number of users who have been active in the following month, and so forth. For newest cohorts you’ll be getting zeroes in some of active_mN columns, since there’s no data on them yet. This is taken into account in processing scripts.

```require(plyr)

# Load SystematicInvestor's plot.table (https://github.com/systematicinvestor/SIT)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
source(con)
close(con)

# Read the data
cohorts
# Let's convert absolute values to percentages (% of the registered users remaining active)
cohort_p as.numeric(df\$active_m0/df\$signed_up), as.numeric(df\$active_m1/df\$signed_up), as.numeric(df\$active_m2/df\$signed_up),
as.numeric(df\$active_m3/df\$signed_up), as.numeric(df\$active_m4/df\$signed_up), as.numeric(df\$active_m5/df\$signed_up),
as.numeric(df\$active_m6/df\$signed_up), as.numeric(df\$active_m7/df\$signed_up), as.numeric(df\$active_m8/df\$signed_up) ))

# Create a matrix
temp = as.matrix(cohort_p[,3:(length(cohort_p[1,])-1)])
colnames(temp) = paste('Month', 0:(length(temp[1,])-1), sep=' ')
rownames(temp) = as.vector(cohort_p\$V1)

# Drop 0 values and format data
temp[] = plota.format(100 * as.numeric(temp), 0, '', '%')
temp[temp == " 0%"] # Plot cohort analysis table
plot.table(temp, smain='Cohort(users)', highlight = TRUE, colorbar = TRUE)
```

This code produces nice visualizations of the cohort analysis as a table: I used articles “Visualizing Tables with plot.table” and “Response to Flowingdata Challenge: Graphing obesity trends” as an inspiration for this R code.

If you want to get nice colours as in the example above, you’ll need to adjust rainbow interval for plot.table. I managed to do it by editing functions code directly from R environment:

`plot.table.helper.color <- edit(plot.table.helper.color)`
```function
(
temp # matrix to plot
){
# convert temp to numerical matrix
temp = matrix(as.double(gsub('[%,\$]', '', temp)), nrow(temp), ncol(temp))

highlight = as.vector(temp)
cols = rep(NA, len(highlight))
ncols = len(highlight[!is.na(highlight)])
cols[1:ncols] = rainbow(ncols, start = 0, end = 0.3)

o = sort.list(highlight, na.last = TRUE, decreasing = FALSE)
o1 = sort.list(o, na.last = TRUE, decreasing = FALSE)
highlight = matrix(cols[o1], nrow = nrow(temp))
highlight[is.na(temp)] = NA
return(highlight)
}
```

Adjust interval in line 11 to 0.5, 0.6 to get shades of blue.
plot.table.helper.colorbar <- edit(plot.table.helper.colorbar)

```function
(
plot.matrix # matrix to plot
)
{
nr = nrow(plot.matrix) + 1
nc = ncol(plot.matrix) + 1

c = nc
r1 = 1
r2 = nr

rect((2*(c - 1) + .5), -(r1 - .5), (2*c + .5), -(r2 + .5), col='white', border='white')
rect((2*(c - 1) + .5), -(r1 - .5), (2*(c - 1) + .5), -(r2 + .5), col='black', border='black')

y1= c( -(r2) : -(r1) )

graphics::image(x = c( (2*(c - 1) + 1.5) : (2*c + 0.5) ),
y = y1,
z = t(matrix( y1 , ncol = 1)),
col = t(matrix( rainbow(len( y1 ), start = 0.5, end = 0.6) , ncol = 1)),
add = T)
}
```

Adjust interval in line 21 to 0.5, 0.6 to get shades of blue.

Now if you want to draw the cycle-like graph:

```# make matrix shorter for the graph (limit to 0-6 months)
temp = as.matrix(cohort_p[,3:(length(cohort_p[1,])-1)])
temp temp[temp == "0"]
library(RColorBrewer)
colnames(temp) = paste('Month', 0:(length(temp[1,])-1), 'retention', sep=' ')
palplot(temp[,1],pch=19,xaxt="n",col=pal,type="o",ylim=c(0,as.numeric(max(temp[,-2],na.rm=T))),xlab="Cohort by Month",ylab="Retention",main="Retention by Cohort")

for(i in 2:length(colnames(temp))) {
points(temp[,i],pch=19,xaxt="n",col=pal[i])
lines(temp[,i],pch=19,xaxt="n",col=pal[i])
}

axis(1,at=1:length(cohort_p\$cohort),labels=as.vector(cohort_p\$cohort),cex.axis=0.75)
legend("bottomleft",legend=colnames(temp),col=pal,lty=1,pch=19,bty="n")
abline(h=(seq(0,1,0.1)), col="lightgray", lty="dotted")
```

This code produces nice visualizations of the cohort analysis as multicolour cycle graph: ## Heat map visualization of sick day trends in Finland with R, ggplot2 and Google Correlate

Inspired by Margintale’s post “ggplot2 Time Series Heatmaps” and Google Flu Trends I decided to use a heat map to visualize sick days logged by HeiaHeia.com Finnish users.

I got the data from our database, filtering results by country (Finnish users only) in a tab separated form with the first line as the header. Three columns contained date, count of sick days logged on that date and count of Finnish users in the service on that date.

```date count(*) user_cnt
2011-01-01 123 12345
2011-01-02 456 67890
...```

Below is R source code for plotting the heat map. I made some small changes to the original code:

• data normalization (line 9): this is specific to the data used in this example
• days of the week have to be 1..7, not 0..6 as returned by \$wday (line 19): dat\$weekday = as.numeric(format(as.POSIXlt(dat\$date),”%u”))
• date format (line 31): week of year calculation required date conversion to POSIX dat\$week <- as.numeric(format(as.POSIXlt(dat\$date),”%W”))
• custom header for the legend (line 39): adding + labs(fill=”per user per day”) allows you to customize legend header
```require(zoo)
require(ggplot2)
require(plyr)

dat<-read.csv("~/data/sick_days_per_day.txt",header=TRUE,sep="\t")</pre>
colnames(dat) <- c("date", "count", "user_cnt")

# normalize data by number of users on each date
dat\$norm_count <- dat\$count / dat\$user_cnt

# facet by year ~ month, and each subgraph will show week-of-month versus weekday the year is simple
dat\$year<-as.numeric(as.POSIXlt(dat\$date)\$year+1900)
dat\$month<-as.numeric(as.POSIXlt(dat\$date)\$mon+1)

# turn months into ordered facors to control the appearance/ordering in the presentation
dat\$monthf<-factor(dat\$month,levels=as.character(1:12),labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)

# the day of week is again easily found
dat\$weekday = as.numeric(format(as.POSIXlt(dat\$date),"%u"))

# again turn into factors to control appearance/abbreviation and ordering
# I use the reverse function rev here to order the week top down in the graph
# you can cut it out to reverse week order
dat\$weekdayf<-factor(dat\$weekday,levels=rev(1:7),labels=rev(c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")),ordered=TRUE)

# the monthweek part is a bit trickier - first a factor which cuts the data into month chunks
dat\$yearmonth<-as.yearmon(dat\$date)
dat\$yearmonthf<-factor(dat\$yearmonth)

# then find the "week of year" for each day
dat\$week <- as.numeric(format(as.POSIXlt(dat\$date),"%W"))

# and now for each monthblock we normalize the week to start at 1
dat<-ddply(dat,.(yearmonthf),transform,monthweek=1+week-min(week))

# Now for the plot
P<- ggplot(dat, aes(monthweek, weekdayf, fill = dat\$norm_count)) +
geom_tile(colour = "white") + facet_grid(year~monthf) + scale_fill_gradient(low="green", high="red") +
opts(title = "Time-Series Calendar Heatmap - HeiaHeia.com sick days logged") + xlab("Week of Month") + ylab("") + labs(fill="per user per day")
P
```

Here are the results. Green indicates the healthiest days with lowest values of sick days logged per user, red indicates the worst days with highest values of sick days logged per user. It’s quite clear that there are seasonal peaks around February, and 2011 was a lot worse than 2012 (one should note that January-February of 2011 were exceptionally cold in Finland). It matches quite well with the coverage in the national press: Flu season reaching peak (Feb’2012), Employers grapple with sick leaves brought by flu wave (Feb’2012).

It’s interesting that there are less sick days logged on the weekends than on the work days, and traditional holiday month of July is the healthiest month of all. (click to see full-sized image)

To get a more formal validation of the data logged by HeiaHeia users, I used Google Correlate lab tool to check that heat map results make sense. I uploaded sick days per user weekly time series and plotted a correlation with Google search queries for “kuumeen hoito” (treatment of fever in Finnish). (click to see full-sized image)

Pearson Correlation Coefficient r between HeiaHeia sick days time series and Google search activity σ (both normalized so that mean is 0 and standard deviation is 1) is 0.8257 – this is a pretty good match.