Adding archer terms

He was, of course, stopped by his brother Tyrion, who said “now Jaime, can’t you see that you chose archers of wildly different statures? Only a fool would not account for this in his analysis.” Grumpily, Jaime re-fit his model with additional terms for the archers:

\[ Y_{ijk} = \mu + \tau\mathbb{1}(i = 2) + \alpha_j + \epsilon_{ik} \]

Here, \(Y_{ijk}\) is the \(k\)-th shot (now \(k\) goes from 1 to 10) from the \(j\)-th archer using the \(i\)-th bow. The \(\alpha_j\) parameters represent archer \(j\)’s baseline skill. Here were Jaime’s results with the new model:

lm1 <- lm(Distance_ ~ Bow_ + Name_, data = JaimeExp1$Data)
summary(lm1)

## 
## Call:
## lm(formula = Distance_ ~ Bow_ + Name_, data = JaimeExp1$Data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.1652  -3.4554  -0.5064   2.5754  13.1611 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          156.427      1.238 126.386  < 2e-16 ***
## Bow_Southwood         -9.767      1.238  -7.891 1.17e-10 ***
## Name_Fabiar Darklyn   38.585      1.516  25.454  < 2e-16 ***
## Name_Margan Drumm     -4.212      1.516  -2.779  0.00741 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.794 on 56 degrees of freedom
## Multiple R-squared:  0.9485, Adjusted R-squared:  0.9457 
## F-statistic: 343.6 on 3 and 56 DF,  p-value: < 2.2e-16

After accounting for the different archer, the estimated effect of the Southwood bow became quite statistically significant. These results weren’t as pleasing to Jaime, who stormed past a smirking Tyrion. “Well? What read your results?” he heard as he slammed the training field gate behind him.

So what happened? There were two separable sources of error around the bow effect: per-shot random error, and variance among archer shooting strengths. At Tyrion’s suggestion, Jaime controlled for this in the second model by including archer effects.

“I drink, and I account for variance in experimental units”

Here’s an illustration of why it’s important to include archer effects, or in general, “experimental unit” effects. Without knowledge of the particular archer, Jaime’s data look like this:

The red shots look somewhat above the blue shots, but the shot points seem to be layered. Here’s the plot with points shaped differently for different archers:

It’s now apparent that the red shots are usually higher than the green shots, and the shots cluster around distinct locations that depend on the marksman. Including marksman terms in the model allows for estimation of the noise variance around the marksman centers, rather than the overall center. This is beneficial from (at least) two standpoints:

1. Interpretation: the variance due to archer skill is clearly not random error, so a model that does not account for the marksman will give an inaccurate estimate of the true per-shot variance.
2. Estimation of \(\tau\): we are better able to “see” the effect of the Southwood bow when we separate the shots by marksman.

After cooling down a little, Jaime realized that, despite not achieving the result he wanted, his findings could be useful to the Lannister army. He decided to run a larger experiment with 50 archers to solidify the results. Here was the outcome:

lm2 <- lm(Distance_ ~ Bow_ + Name_, data = JaimeExp2$Data)
summary(lm2)

## 
## Call:
## lm(formula = Distance_ ~ Bow_ + Name_, data = JaimeExp2$Data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.4424  -3.3524  -0.1689   3.2859  14.8538 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            168.8198     1.1263 149.883  < 2e-16 ***
## Bow_Southwood          -10.4792     0.3154 -33.221  < 2e-16 ***
## Name_Aldo Apperford    -12.2736     1.5772  -7.782 1.86e-14 ***
## Name_Antorn Feller      34.5155     1.5772  21.884  < 2e-16 ***
## Name_Arren Durrandon    29.3949     1.5772  18.637  < 2e-16 ***
## Name_Ashtin Orme        11.5839     1.5772   7.345 4.44e-13 ***
## Name_Barret Foral      -16.6079     1.5772 -10.530  < 2e-16 ***
## Name_Barrock Ryswell     6.4531     1.5772   4.091 4.65e-05 ***
## Name_Brack Mallister   -53.4039     1.5772 -33.860  < 2e-16 ***
## Name_Brandeth Erenford   6.1860     1.5772   3.922 9.41e-05 ***
## Name_Brayan Trante     -35.2405     1.5772 -22.344  < 2e-16 ***
## Name_Charres Apperford   0.3317     1.5772   0.210 0.833494    
## Name_Clarrik Parge       8.2105     1.5772   5.206 2.37e-07 ***
## Name_Codin Sadelyn      37.8102     1.5772  23.973  < 2e-16 ***
## Name_Culler Vyrwel       6.8461     1.5772   4.341 1.57e-05 ***
## Name_Curtass Shiphard  -33.6791     1.5772 -21.354  < 2e-16 ***
## Name_Denzin Garner     -21.6985     1.5772 -13.758  < 2e-16 ***
## Name_Deran Swygert       5.3668     1.5772   3.403 0.000695 ***
## Name_Donnal Brackwell  -18.5525     1.5772 -11.763  < 2e-16 ***
## Name_Dontar Brask        3.0880     1.5772   1.958 0.050534 .  
## Name_Dran Kenning       -4.3976     1.5772  -2.788 0.005406 ** 
## Name_Dwan Plumm        -14.3560     1.5772  -9.102  < 2e-16 ***
## Name_Eddard Meadows    -24.9352     1.5772 -15.810  < 2e-16 ***
## Name_Erock Conklyn      11.5451     1.5772   7.320 5.28e-13 ***
## Name_Evin Rane         -33.1439     1.5772 -21.014  < 2e-16 ***
## Name_Fabiar Darklyn    -33.2426     1.5772 -21.077  < 2e-16 ***
## Name_Harden Blest      -30.0566     1.5772 -19.057  < 2e-16 ***
## Name_Howar Clifton      -2.6749     1.5772  -1.696 0.090221 .  
## Name_Jaesse Cantell    -15.7270     1.5772  -9.971  < 2e-16 ***
## Name_Jares Grell         3.6472     1.5772   2.312 0.020966 *  
## Name_Jarvas Bywater      8.6005     1.5772   5.453 6.32e-08 ***
## Name_Jeffary Mudd       -4.4122     1.5772  -2.797 0.005254 ** 
## Name_Jeran Parsin       -8.9081     1.5772  -5.648 2.14e-08 ***
## Name_Kober Maeson       -5.4609     1.5772  -3.462 0.000559 ***
## Name_Koryn Staunton     32.5297     1.5772  20.625  < 2e-16 ***
## Name_Marak Herston     -51.0331     1.5772 -32.357  < 2e-16 ***
## Name_Mikal Netley      -11.2779     1.5772  -7.151 1.72e-12 ***
## Name_Mykal Ryser       -21.1308     1.5772 -13.398  < 2e-16 ***
## Name_Nelsor Baxter     -19.2805     1.5772 -12.225  < 2e-16 ***
## Name_Portar Parge       36.7979     1.5772  23.331  < 2e-16 ***
## Name_Raman Brakker      36.6719     1.5772  23.251  < 2e-16 ***
## Name_Randar Condon      -4.3076     1.5772  -2.731 0.006428 ** 
## Name_Roberd Spyre        3.6156     1.5772   2.292 0.022098 *  
## Name_Rolan Morrass      25.9665     1.5772  16.464  < 2e-16 ***
## Name_Samn Whitehill    -18.9911     1.5772 -12.041  < 2e-16 ***
## Name_Sarrac Ridman      -9.2904     1.5772  -5.890 5.34e-09 ***
## Name_Seamas Knigh        3.9833     1.5772   2.526 0.011713 *  
## Name_Seldan Harker     -21.1080     1.5772 -13.383  < 2e-16 ***
## Name_Sharun Mudd       -27.9034     1.5772 -17.692  < 2e-16 ***
## Name_Zandren Kneight   -53.3039     1.5772 -33.797  < 2e-16 ***
## Name_Zarin Goodbrother -21.9110     1.5772 -13.892  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.988 on 949 degrees of freedom
## Multiple R-squared:  0.9585, Adjusted R-squared:  0.9564 
## F-statistic: 438.8 on 50 and 949 DF,  p-value: < 2.2e-16

Having disregarded much of his scientific education, Jaime felt quite overwhelmed with this massive list of parameter estimates. In particular, he was confused by the few archers without a statistically significant estimate. Should he re-run the experiment without them? (Note, the standard errors for the archer effect estimates are the same because they all fired the same number of shots.)

This time, Cersei was hanging around the training grounds. Overhearing Jaime, she decided to step in. “Jaime, don’t trouble yourself with those worrisome parameter estimates for the archers. I looked at your charts from the first experiment. It is certainly proper to account for archers in the model, but the only value we care about is the variance of their skills.” Later that afternoon she developed Westeros’s first Linear Mixed Effects model. In this model, archer effects are considered to be random draws from a larger population.

“With these equations I shall secure the reputation of the Lannister name in both might and maths”

The mixed-effects model and sources of randomness

Cersei’s key observation is that Jaime chose 50 archers out of the large population of archers in and around King’s Landing. Furthermore, he did so in a random fashion, distributing the bow test across a representative population of archers. This randomness can be incorporated in the model by treating archer effects as Normal random variables with a common, fixed variance: \[ Y_{ijk} = \mu + \tau\mathbb{1}(i = 2) + \alpha_j + \epsilon_{ik}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\alpha_j\sim\mathcal{N}(0, \sigma^2_a) \] The \(\alpha_j\)’s are now not of interest to estimate; they exist to separate out the contribution of marksman skill to per-shot variance. This “variance component”, \(\sigma^2_a\), is now the archer parameter of interest.

Note that, though the \(\alpha_j\)’s as modeled as random variables, they do not change shot-to-shot. Mechanistically, they are still fixed features of the archer. We regard them as random only to focus on the underlying distribution of archer skill, and its effect on the predictive model.

Fitting the mixed-effects model

Let’s fit the model to Jaime’s 2nd set of experimental data:

lme2 <- lmer(Distance_ ~ Bow_ + (1 | Name_), data = JaimeExp2$Data)

The switch to mixed effects modeling shifts the focus from inference on the individual archer to inference on the population of archers. This focus is important in applications from many fields, like biology (where observations may come from many organisms), economics (where observations may come from many regions or agents), or clinical trials (where observations may come from many different labs). In our example, it makes sense for Cersei to want to learn about the variance in city-wide marksman skill along with the effect of Southbows. The estimate \(\hat\sigma^2_a\) gives her this:

VarCorr(lme2)

##  Groups   Name        Std.Dev.
##  Name_    (Intercept) 22.9698 
##  Residual              4.9875

confint(lme2)[1:2, ]

##            2.5 %    97.5 %
## .sig01 18.907078 28.054556
## .sigma  4.768896  5.217741

The output above gives the estimates for \(\sigma_a\) and \(\sigma_e\), and the assocated confidence intervals. These parameters don’t exist in a fixed-effects model.

Accuracy and speed of the mixed effects model

This last section compares the performance of the Binomial GLMM to that of standard fixed-effects models, with respect to both predictive potential and compuation time.

As mentioned just before the CTR example, a common addition to large-scale fixed-effect regression modeling is a LASSO penalty. Roughly, the LASSO penalty zeros out user coefficients for users with effects not significantly different from the mean user effect. Fitting a LASSO is a little tricky in R if you have multi-level factor variables, especially if one of the variables has many, many levels (like our “user” variable). In case you’re interested in the best way to do this, I’ll print this chunk of code:

suppressPackageStartupMessages(library(glmnet))
id.matrix <- sparse.model.matrix(~ CTRdata$id - 1)
Sparse.Data <- cbind(id.matrix, CTRdata$display)
weights <- CTRdata$expos
nclicks <- CTRdata$clicks * weights
y <- unname(cbind(weights - nclicks, nclicks))
set.seed(12345)
CTR_lasso <- cv.glmnet(Sparse.Data, y, family = "binomial")

I’ll also fit two other fixed-effects models: a Ridge regression, which uses another form of parameter penalization, and a straight-up fixed-effect logistic model.

set.seed(12345)
CTR_ridge <- cv.glmnet(Sparse.Data, y, family = "binomial", alpha = 0)
CTR_glm <- glm(y[ , 2:1] ~ display + id, data = CTRdata, family = "binomial")

The \(\lambda\) regularization parameter for the Ridge and LASSO models were chosen by 10-fold cross-validation, as is default for the cv function. Here are the effect estimates for each of these models:

##     LASSO     Ridge       GLM 
## 0.1247575 0.1371491 0.5138967

The true display effect is 0.5, so the effect of Ridge and LASSO regularization is apparent. Though Ridge and LASSO have their uses, at least in this case they’re preventing an interpretable estimate of display effect. Since the logistic is a full model fit, we should not be surprised that its point estimate for \(\tau\) is close to the truth.

How do the predictive powers of these approaches compare? Let’s generate 50 “future” data sets from the same population of users. For each new data set, we’ll predict the logit means using the fitted model from the original data. To compare the prediction accuracies of the approaches, we’ll calculate the mean squared prediction error¹. Note that, for the future data sets, it is possible for users to be assigned to Treatment when in the original data set they were Control, and vice-versa. This provides some nice heterogenity in this comparative study. The following are the means and standard errors of the methods’ MSPEs:

## Warning: 'cBind' is deprecated.
##  Since R version 3.2.0, base's cbind() should work fine with S4 objects

# MSPE means
apply(MSEP_mat, 1, mean)
# MSPE standard errors
apply(MSEP_mat, 1, sd)

At least for this number of users, the difference between methods appears vast. Let’s explore the MSPE for varying numbers of users. I ran the following experiment for a range of user sample sizes:

1. Sample the desired number of users from the population and generate data
2. Fit the GLMM, GLM, LASSO, and Ridge models to the data
3. Generate 50 new data sets from the same originally chosen users
4. Calculate and store MSPE for each data set, for each method

These were the results:

Some conclusions:

1. The LASSO and Ridge models are worse than the GLM and GLMM models. Possibly, the model parsimony provided by the former methods is actually harmful in this setting. This result possibly reflects the same phenomenon at play in the Google blog post.

2. The MSPEs of the GLM and GLMM fits are nearly identical. This should not be surprising, since the model matrix equations are essentially the same.

3. While GLM and GLMM fits are identical, we seem to be gaining something computationally by fitting a mixed-effects model. This is probably only due to the particular R packages involved. Specifically, I think the fact that the GLM estimates come with standard errors is the main bottleneck. If we just cared about prediction, there is probably an equally fast GLM implementation. That being said, the choice to disregard user effects as parameters of interest provided the additional benefit of a faster off-the-shelf solution.

Caveats:

For the fixed-effect models, there might be better rules-of-thumb or different regularization techniques that would make more sense for situations like the CTR example. Maybe these would improve their performance, but I’m not sure. The results shown are from the most straightforward off-the-shelf implementations for some of the most common regularization approaches.

Wait, you work at Google, why are you responding to one of their blog posts on a personal blog? I wrote this post before I got the job (almost a year before my first day), and published it on my old webpage. I recently completely revamped my website and blogging workflow (everything uses blogdown now), so I’m using this as an easy “first” post. As mentioned at beginning, I just tweaked the writing style; the analysis is unchanged. I think it’s still a good reference for these very general methods, and I’ve gone back to it a few times myself!

Further reading and references:

1. http://glmm.wikidot.com/faq
2. http://datascienceplus.com/introduction-to-bootstrap-with-applications-to-mixed-effect-models
3. http://www.ats.ucla.edu/stat/r/dae/melogit.htm
4. Variance Components (Searle, Casella, and McCulloch, 1992)
5. Generalized Linear Models, 2nd Edition (McCullagh and Nelder, 1989)

Code examples:

To make the CTR simulation data:

make_click_data <- function (p, # Number of people to choose 
                             mu, # Baseline logit-probability of click 
                             tau, # Effect of new display
                             sig_a, # Variance of random effects
                             pois_mean, # Mean number of impressions
                             effects = NULL) { # In case you actually want to fix effects
  
  if (is.null(effects))
    effects <- rnorm(p, 0, sig_a)
  
  # Preparing design
  treat <- rbinom(p, 1, 0.5)
  treat_expos <- rpois(p, pois_mean * 0.5) + 1
  null_expos <- rpois(p, pois_mean * (1 - 0.5)) + 1
  id <- 1:p
  expos <- treat_expos + null_expos
  expos <- c(expos[!treat == 1], null_expos[treat == 1], treat_expos[treat == 1])
  id <- c(id[!treat == 1], rep(id[treat == 1], 2))
  treat <- c(treat[!treat == 1], rep(0, sum(treat)), rep(1, sum(treat)))
  means <- mu + treat * tau + effects[id]
  probs <- exp(means) / (1 + exp(means))
  
  # Getting response and making data matrix
  Y <- rbinom(length(probs), expos, probs)
  dataMat <- data.frame("clicks" = Y,
                        "expos" = expos,
                        "display" = treat,
                        "id" = as.factor(id))
  
  return(list("Data" = dataMat,
              "Effects" = effects))
  
}

# Make simulated data
set.seed(1)
CTR_Exp1 <- make_click_data(p = 1e3,
                            mu = -4.5,
                            tau = 0.5,
                            sig_a = 5,
                            pois_mean = 20)

To calculate MSPE of the methods over many repetitions of the data from the original simulation model:

user_effects <- CTR_Exp1$Effects
nsims <- 50
MSEP_mat <- matrix(0, 4, nsims)
rownames(MSEP_mat) <- c("GLMM", "LASSO", "Ridge", "GLM")
for (i in 1:nsims) {
  
  # Generating data
  new_CTR <- make_click_data(p = 1e3, 
                             mu = -4.5,
                             tau = 0.5,
                             sig_a = 5,
                             pois_mean = 20,
                             effects = user_effects)
  
  # Prepping LASSO data
  new_id.matrix <- sparse.model.matrix(~ new_CTR$Data$id - 1)
  new_Sparse.Data <- cBind(new_id.matrix, new_CTR$Data$display)
  
  # Getting y data
  new_weights <- new_CTR$Data$expos
  new_nclicks <- new_CTR$Data$clicks
  y_new <- unname(cbind(new_weights - new_nclicks, new_nclicks))
  
  calc_MSEP <- function (pred) {
    err0 <- pred^2 * y_new[ , 1]
    err1 <- (1 - pred)^2 * y_new[ , 2]
    return((sum(err0) + sum(err1)) / sum(new_weights))
  }
  
  # Calculating prediction logit means
  lme_pred <- predict(CTR_lme1, newdata = new_CTR$Data)
  lasso_pred <- predict(CTR_lasso, newx = new_Sparse.Data)
  ridge_pred <- predict(CTR_ridge, newx = new_Sparse.Data)
  glm_pred <- predict(CTR_glm, newdata = new_CTR$Data)
  
  # Calculating prediction probs
  lme_pred_probs <- exp(lme_pred) / (1 + exp(lme_pred))
  lasso_pred_probs <- exp(lasso_pred) / (1 + exp(lasso_pred))
  ridge_pred_probs <- exp(ridge_pred) / (1 + exp(ridge_pred))
  glm_pred_probs <- exp(glm_pred) / (1 + exp(glm_pred))
  
  # Calculating prediction probs
  prob_mat <- cbind(exp(lme_pred) / (1 + exp(lme_pred)),
                    exp(lasso_pred) / (1 + exp(lasso_pred)),
                    exp(ridge_pred) / (1 + exp(ridge_pred)),
                    exp(glm_pred) / (1 + exp(glm_pred)))
  
  MSEP_mat[ , i] <- unname(apply(prob_mat, 2, calc_MSEP))
  
}