Efficient implementation of average predictive comparison (APC, Gelman & Pardoe 2007)?

[Guido Biele]

Hello,

I’d appreciate some feedback on the implementation of average predictive comparisons (APC; Gelman & Pardoe 2007) in stan (for now only for non-categorical variables).

######### if you know APC, jump over this section #######

Background: In a nutshell, the APC provide an estimate about how much on average a response variable y changes for a one unit change in the predictor variable u, given all other predictor variables v are fixed.

For example, one could examine the effect of the age, number of publications, number of citations,  on the response-variable “has tenure” by running a logistic regression. Calculating APCs for each variable would then  provide an estimate by how many percent the probability to be tenured increases for each age year, number of publications, number of publications, and number of citations, respectively.

Equation 5 in Gelman & Pardoe (2007), describes an algorithm to estimate APCs and uses the following notation:

• u, the variable of interest

• v, a vector with the other variables

• theta, the model parameters

• yhat, the model prediction

• n, the number data columns

• S, the number of mcmc sample

to estimate the APC for a variable u as

here, E(yhat|u_j,v_i,theta_s) is a known function mapping data and parameters to predicted values (e.g. the inverse logit for logistic regressions). w_ij is a Mahalanobis-distance based weight that insuring comparisons u_i-u_j for which v_i and v_j are similar have a greater weight than comparisons u_i-u_j for which v_i and v_j are dissimilar. The term (u_j-u_i) in the denominator has the function to give smaller weight to comparisons with small differences of u, because these are potentially unstable due to noise or incomplete models. Finally, Gelman and Pardoe use the term sign(u_j-u_i) because their algorithm only considers cases where u_i<u_j. Otherwise the denominator and numerator would each sum to zero.

################### continue here ##################

My APC implementation deviates a bit from a literal implementation of the equation above, and I’d like to know if if my implementation makes sense.

In a nutshell, i prepare pairs for predictive comparisons and weights in the transformed parameters block calculate pair by pair predictive comparisons and their weighted average in the generated quantities block.

(I) Generate pairs ui,uj for predictive comparisons

Because I have too many data rows to make all pair comparisons in a reasonable amount of time, I submit 2 vectors u1s and u2s to stan, which contain row-numbers of the rows to be used to calculate the APCs. Pairs are created from indices in u1s and u2s by combining each u1 in u1s with each u2 in u2s. to avoid the problem of summation to zero, vectors u1s and u2s must be non-overlapping.

In R for example i would create non-overlapping indices as follows:

uidxs = sample(1:n,500,replace = F)  # sample 500 out of n rows
u1s = uidxs[1:250]
u2s = uidxs[251:500]

For smaller data sets one could define all pairs in stan by looping only through the lower triangle of the matrix of all possible pairs u_i,u_j:

for (u1 in 1:(n-1)) {
 for (u2 in (u1+1):n){
 ...
 }
}

(II) Calculate weights in the transformed data block.

I define weights as w_ij*|u_j-u_i|, where w_ij  is the  Mahalanobis distance based weight proposed by Gelman & Pardoe, and |u_j-u_i| downweights pairs with small differences in u.

While this also follows Gelman and Pardoe, it still seems a bit ad hoc to me, because one could also omit this second term (when the model is complete and there is not to much noise) or one could take the sqrt(|ui-uj|) of it to avoid overweighting of extreme predictor values. (see also the discussion here: http://www.davidchudzicki.com/predcomps/more-pairs-and-weights.html)

Because below I increment the APC comparison by comparison (i.e. I don’t first calculate separate sums for denominator and numerator), my weights have to sum to 1.

stan code:

(the full model code is attached)

transformed data {
  // matrix for weights
  // K = number of predictor variables u (where the intercept is the first predictor)
  // nus = number of indices first and second part of pair
  matrix[K-1,nus*nus] weights;  

  // matrix for delta_u = u1-u2 to avoid calculating it repeatedly later  
  matrix[K-1,nus*nus] delta_u;

  // matrices to store indices of nonzero u1-u2 differences  
  int<lower=0,upper=N> u1idx[K-1,nus*nus];
  int<lower=0,upper=N> u2idx[K-1,nus*nus]

  row_vector[K-2] delta_v;
  matrix[1,1] mhlbd_sq;

  // array to store number of nonzero differences u1-u2
  int<lower=1> n_pairs[K-1];

  n_pairs <- rep_array(0,K-1);
  weights <- rep_matrix(0.0, K-1, nus*nus);


  for (u in 2:K){ // loop through predictor variables
    for (u1 in 1:nus){
      for (u2 in 1:nus){
        if (x[u1s[u1],u] != x[u2s[u2],u]) {
          n_pairs[u-1] <- n_pairs[u-1]+1;

          u1idx[u-1,n_pairs[u-1]] <- u1s[u1];
          u2idx[u-1,n_pairs[u-1]] <- u2s[u2];

          delta_u[u-1,n_pairs[u-1]] <- x[u1s[u1],u] - x[u2s[u2],u];

          // calculate vector delta_v  for Mahalanobis distance
          // x is the design matrix for the regression, with the intercept in the first column
         if( u == 2) {
             delta_v <- sub_row(x,u1s[u1],u+1,K-u) - sub_row(x,u2s[u2],u+1,K-u);
         } else if (u < K) {
             delta_v <- append_col(sub_row(x,u1s[u1],2,u-2),sub_row(x,u1s[u1],u+1,K-u)) -
                        append_col(sub_row(x,u2s[u2],2,u-2),sub_row(x,u2s[u2],u+1,K-u));
         } else {
             delta_v <- sub_row(x,u1s[u1],2,u-2) - sub_row(x,u2s[u2],2,u-2);
         }


         // square of mahalabonis distance between vectors v for u1 and u2
         mhlbd_sq <- to_matrix(delta_v) * iCVs[u-1] * to_matrix(delta_v)';

         // calculate weights from 2 parts: 1) |u1-u2|, 2) square of mahalabonis distance
         weights[u-1,n_pairs[u-1]] <- fabs(delta_u[u-1,n_pairs[u-1]])*1/(1+mhlbd_sq[1,1]);
        }
      }
    }
    // normalize weights so that they sum to 0
    weights[u-1] <- weights[u-1]/sum(weights[u-1]);
  }
  print(n_pairs);
}

(III) Calculate APCs

The basic idea to calculate the APC for a variable u is simple:

• loop through all comparisons,

• calculate a pair’s predictive comparison as (E(y|u_j,v_i,theta_s)-(E_y|u_i,v_i,theta_s))/(u_j-u_i)

• multiply the pair's predictive comparison with its weight (which was prepared in the transformed parameters block)

• increment the APC with the weighted predictive comparison

stan code:

generated quantities {
  // array to store APCs
  real APC[K-1];
  APC <- rep_array(0,K-1);

  for (u in 2:K){ // loop through predictor variables
    for (k in 1:n_pairs[u-1]) { // loop through pairs with nonzero difference u1-u2
      // increment            eta=log(mu) of a neg. binomial  * weight
      APC[u-1] <- APC[u-1] + (exp(eta[u1idx[u-1,k]])-exp(eta[u2idx[u-1,k]]))/delta_u[u-1,k] * weights[u-1,k];
    }
  }
}

Finally, I am uncertain about the use of matrices vs arrays. I could glean from the stan-reference that arrays have the advantage that values can be returned by reference. On the other hand, I also read that "traversing matrices is more efficient than traversing arrays" (though I haven't checked yet if I am traversing them in the right order). I am not sure what is more relevant here.

Speed is certainly an issue, because with the current code the sampling phase is slower than the warmup phase.

Best - Guido

[Andrew Gelman]

Hi, Guido.  I’m cc-ing David Chudzicki because he’s been thinking a lot about this recently and he could have some useful thoughts.

Andrew

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<HZAxNB2RegressionVapc.stan>

[Bob Carpenter]

>> On Jan 13, 2015, at 3:57 PM, Guido Biele <guido...@gmail.com> wrote:

>> ...

>> Finally, I am uncertain about the use of matrices vs arrays. I could glean from the stan-reference that arrays have the advantage that values can be returned by reference. On the other hand, I also read that "traversing matrices is more efficient than traversing arrays" (though I haven't checked yet if I am traversing them in the right order). I am not sure what is more relevant here.
>> Speed is certainly an issue, because with the current code the sampling phase is slower than the warmup phase.
>>

The matrix/vector decision is a fairly minor one in terms of speed
in most cases.

Matrices are stored column major, so you want this order
for speed:

for (n in 1:N)
for (m in 1:M)
... do something with x[m,n] ...

If you need to pull a row out of a matrix and use it for something,
it's better to use an array of row vectors (or just an array of vectors
if you'd otherwise need to transform).

- Bob

[David]

Seems reasonable, but I'm surprised by the decision to create the pairs/weights inside Stan rather than do that step in another language and provide it to Stan as input data. I don't know if that step is slow, but if it is there'd be an advantage to doing it outside Stan and saving it once and for all. Then as you tweak the model, you don't need to do that step again each time you run Stan.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

David J. Chudzicki
http://davidchudzicki.com
dch...@gmail.com
(518) 366-7303

[Guido Biele]

Hi David,

it is a bit a hassle to program the generation of pairs in stan. But I am still playing around with how many pairs to generate and find it convenient to be able to simply submit the data set and let stan to the rest. Creating the pairs & weights in stan takes a bit time, but this is negligible compared to sampling times. (another issue is that matrices with pairs can become pretty large if one has thousands of data points and 20 or more regressors).

I'm still finalizing the code and will send around an update early this week.

Cheers - Guido

 

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/NYodG--ecSI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.