I would like to suggest a (standard) preliminary analysis to remove the principal effects of (a) variation among users, (b) the typical response among all users to the change, and (c) typical variation from one time period to the next.
A simple (but by no means the best) way to do this is to perform a few iterations of "median polish" on the data to sweep out user medians and time period medians, then smooth the residuals over time. Identify the smooths that change a lot: they are the users you want to emphasize in the graphic.
Because these are count data, it's a good idea to re-express them using a square root.
As an example of what can result, here is a simulated 60-week dataset of 240 users who typically undertake 10 to 20 actions per week. A change in all users occurred after week 40. Three of these were "told" to respond negatively to the change. The left plot shows the raw data: counts of action by user (with users distinguished by color) over time. As asserted in the question, it's a mess. The right plot shows the results of this EDA--in the same colors as before--with the unusually responsive users automatically identified and highlighted. The identification--although it is somewhat ad hoc--is complete and correct (in this example).
Here is the R
code that produced these data and carried out the analysis. It could be improved in several ways, including
Using a full median polish to find the residuals, rather than just one iteration.
Smoothing the residuals separately before and after the change point.
Perhaps using a more sophisticated outlier detection algorithm. The current one merely flags all users whose range of residuals is more than twice the median range. Albeit simple, it is robust and appears to work well. (A user-settable value, threshold
, can be adjusted to make this identification more or less stringent.)
Testing nevertheless suggests this solution works well for a wide range of user counts, 12 - 240 or more.
n.users <- 240 # Number of users (here limited to 657, the number of colors)
n.periods <- 60 # Number of time periods
i.break <- 40 # Period after which change occurs
n.outliers <- 3 # Number of greatly changed users
window <- 1/5 # Temporal smoothing window, fraction of total period
response.all <- 1.1 # Overall response to the change
threshold <- 2 # Outlier detection threshold
# Create a simulated dataset
set.seed(17)
base <- exp(rnorm(n.users, log(10), 1/2))
response <- c(rbeta(n.users - n.outliers, 9, 1),
rbeta(n.outliers, 5, 45)) * response.all
actual <- cbind(base %o% rep(1, i.break),
base * response %o% rep(response.all, n.periods-i.break))
observed <- matrix(rpois(n.users * n.periods, actual), nrow=n.users)
# ---------------------------- The analysis begins here ----------------------------#
# Plot the raw data as lines
set.seed(17)
colors = sample(colors(), n.users) # (Use a different method when n.users > 657)
par(mfrow=c(1,2))
plot(c(1,n.periods), c(min(observed), max(observed)), type="n",
xlab="Time period", ylab="Number of actions", main="Raw data")
i <- 0
apply(observed, 1, function(a) {i <<- i+1; lines(a, col=colors[i])})
abline(v = i.break, col="Gray") # Mark the last period before a change
# Analyze the data by time period and user by sweeping out medians and smoothing
x <- sqrt(observed + 1/6) # Re-express the counts
mean.per.period <- apply(x, 2, median)
residuals <- sweep(x, 2, mean.per.period)
mean.per.user <- apply(residuals, 1, median)
residuals <- sweep(residuals, 1, mean.per.user)
smooth <- apply(residuals, 1, lowess, f=window) # Smooth the residuals
smooth.y <- sapply(smooth, function(s) s$y) # Extract the smoothed values
ends <- ceiling(window * n.periods / 4) # Prepare to drop near-end values
range <- apply(smooth.y[-(1:ends), ], 2, function(x) max(x) - min(x))
# Mark the apparent outlying users
thick <- rep(1, n.users)
thick[outliers <- which(range >= threshold * median(range))] <- 3
type <- ifelse(thick==1, 3, 1)
cat(outliers) # Print the outlier identifiers (ideally, the last `n.outliers`)
# Plot the residuals
plot(c(1,n.periods), c(min(smooth.y), max(smooth.y)), type="n",
xlab="Time period", ylab="Smoothed residual root", main="Residuals")
i <- 0
tmp <- lapply(smooth,
function(a) {i <<- i+1; lines(a, lwd=thick[i], lty=type[i], col=colors[i])})
abline(v = i.break, col="Gray")
facet_wrap
function of ggplot2 to create a block of 4 x 5 charts (4 rows, 5 columns - adjust depending on the desired aspect ratio) with ~5 users per chart. That should be clear enough and you could scale it up to about 10 users per chart, giving room for 200 with a 4x5 plot or 360 with a 6x6 plot. $\endgroup$