title: Cross-tabulating data with a function
body: <p>I have data in tree columns</p>
  
  <pre><code>set.seed(42)
  N = 1000
  XYp = as.data.frame(matrix(cbind(round(runif(N)*100),
                                   round(runif(N)*1000+1000),
                                   round(runif(N),2)),N,3))
  colnames(XYp) &lt;- c('X','Y','p')
  </code></pre>
  
  <p>Now I would like to cross-tabulate the data based on deciles in 2 dimension:</p>
  
  <pre><code>colX_deciles = quantile(data[,'X'], probs=seq(0,1,1/10))
  colY_deciles = quantile(data[,'Y'], probs=seq(0,1,1/10))
  XYp['X_decile'] &lt;- findInterval(XYp[,'X'],colX_deciles,all.inside = TRUE)
  XYp['Y_decile'] &lt;- findInterval(XYp[,'Y'],colY_deciles,all.inside = TRUE)
  &gt; colX_deciles
     0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
    0.0   9.9  18.0  29.0  39.0  48.0  57.0  69.0  79.2  91.0 100.0 
  &gt; colY_deciles
      0%    10%    20%    30%    40%    50%    60%    70%    80%    90%   100% 
  1000.0 1088.0 1180.0 1279.4 1392.0 1502.5 1602.4 1711.3 1805.2 1902.0 2000.0 
  </code></pre>
  
  <p>I have figured out that it is possible to calculate the sum of elements in column <code>p</code> using <code>xtabs</code>:</p>
  
  <pre><code>&gt; xtabs(p ~ X_decile + Y_decile, XYp)
          Y_decile
  X_decile    1    2    3    4    5    6    7    8    9   10
        1  2.57 8.74 5.51 5.74 4.40 1.77 5.79 3.43 4.66 3.80
        2  6.43 4.25 7.29 5.41 3.08 4.43 8.70 2.62 3.37 4.45
        3  1.99 2.80 7.54 2.56 5.02 4.30 7.99 2.03 4.91 6.28
        4  4.53 4.90 8.04 3.49 2.25 2.87 7.47 5.41 3.54 9.28
        5  2.32 5.82 7.18 4.58 5.39 2.26 0.59 9.61 5.91 5.37
        6  7.70 5.50 6.45 7.83 4.65 8.45 1.70 6.40 4.88 4.32
        7  7.05 3.87 3.54 3.79 6.15 5.55 6.31 2.31 3.42 6.14
        8  4.43 4.50 3.04 3.62 9.92 5.66 3.75 7.01 4.92 7.08
        9  3.67 5.56 3.56 7.92 5.05 5.00 3.64 6.74 5.85 3.26
        10 5.75 3.17 9.50 5.44 3.64 6.13 3.18 5.93 6.18 3.71
  </code></pre>
  
  <p>But how to elegantly apply any function to the cross-tabulated matrix element and get the results, for example <code>avg(p)</code> in the following manner? :</p>
  
  <pre><code>&gt; xtabs(mean(p) ~ X_decile + Y_decile, XYp)
  Error in model.frame.default(formula = mean(p) ~ X_decile + Y_decile,  : 
    variable lengths differ (found for 'X_decile')
  </code></pre>
  
  <p>As a bonus, the values of <code>colX_deciles[1:10]</code> and <code>colY_deciles[1:10]</code> could be set as row names and column names, respectively.</p>