title: Fast and elegant way to calculate fold change between several groups for many variables?
body: <p>Out of curiosity I have been playing with several ways to calculate fold changes and I am trying to find the fastest and the most elegant way to do that (hoping that would also be the same solution).</p>
  
  <p><br></p>
  
  <p>The kind of matrix I am interested in would look like this:</p>
  
  <pre><code># Some data
  nvars &lt;- 10000
  nsamples &lt;- 500
  sample_groups &lt;- 5
  MAT &lt;- replicate(nvars, runif(n=nsamples))
  </code></pre>
  
  <p>And a grouping vector that look like this:</p>
  
  <pre><code>f &lt;- rep_len(1:sample_groups, nsamples)
  f &lt;- LETTERS[f]
  </code></pre>
  
  <p>The output I would like in the end from the above input is a 10 x 10,000 matrix with one row for each combination of the levels in <code>f</code>.</p>
  
  <p><br><br></p>
  
  <p>To do this the first task is to calculate the mean for each group for all column.
  I have come up with 4 possible methods to do that:</p>
  
  <hr>
  
  <pre><code># Settings
  aggr_FUN  &lt;- mean
  combi_FUN &lt;- function(x,y) "/"(x,y) 
  
  # helper function
  pasteC &lt;- function(x,y) paste(x,y,sep=" - ")
  </code></pre>
  
  <p><br></p>
  
  <pre><code>#A1. Loop
  system.time({
  f_un &lt;- unique(f)
  temp1 &lt;- matrix(NA,nrow = length(f_un),ncol=ncol(MAT))
  rownames(temp1) &lt;- f_un
  
  for(i in 1:length(f_un)){
    temp1[i,] &lt;- apply(MAT[f_un[i] == f,,drop=FALSE],2,aggr_FUN)
  }
  })
  
  
   user  system elapsed 
     0.41    0.00    0.41 
  </code></pre>
  
  <p><br><br></p>
  
  <pre><code>#A2. aggregate
  system.time({
  temp2 &lt;- aggregate(. ~ class, data = cbind.data.frame(class=f,MAT), aggr_FUN)
  })
  
  
   user  system elapsed 
     7.76    0.05    7.81 
  </code></pre>
  
  <p><br><br></p>
  
  <pre><code>#A3. reshape2
  library(reshape2)
  system.time({
  temp3 &lt;- recast(data.frame(class=f,MAT),class ~ variable,id.var="class",aggr_FUN)
  })
  
  
    user  system elapsed 
     1.82    0.30    2.12 
  </code></pre>
  
  <p><br><br></p>
  
  <pre><code>#A4. purrr
  library(purrr)
  system.time({
  temp4 &lt;- data.frame(class = f, MAT) %&gt;%
    slice_rows("class") %&gt;%
    by_slice(map, aggr_FUN)
  })
  
  
     user  system elapsed 
     0.47    0.00    0.47 
  </code></pre>
  
  <hr>
  
  <p><br><br></p>
  
  <p>As you can see the loop is actually the fastest solution with the <code>purrr</code> package being only slightly slower. <code>recast</code> is 5 times slower and <code>aggregate</code> is a clear looser.
  I also tried the dplyr package but that turned out to be very slow for some reason (<a href="https://github.com/hadley/dplyr/issues/1395" rel="nofollow">https://github.com/hadley/dplyr/issues/1395</a>).</p>
  
  <p><code>purrr</code> is both fast and elegant so for this part it is a pretty academic exercise looking for a better way.</p>
  
  <p>At this point the output we have is:</p>
  
  <pre><code>&gt; temp1[,1:6]
         [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
  A 0.4804964 0.4779168 0.5292458 0.4401357 0.4728515 0.5009800
  B 0.4819612 0.5260592 0.5291887 0.5095620 0.4840777 0.4792213
  C 0.4661714 0.4886010 0.5006018 0.5061170 0.5058892 0.5432819
  D 0.4566942 0.4519988 0.5334207 0.4912822 0.4542889 0.4898384
  E 0.4967948 0.5630683 0.4941777 0.5239327 0.5045152 0.5227140
  </code></pre>
  
  <p><br><br>
  So if you are still reading here comes the more challenging part. We need to calculate the fold change between all combinations of the groups/rows.</p>
  
  <p>I found two ways to do this:</p>
  
  <pre><code>#B1. by loop
  combs &lt;- t(combn(as.character(f_un),2))
  combi_FUN_vec &lt;- Vectorize(combi_FUN)
  
  out &lt;- matrix(NA,nrow = nrow(combs),ncol=ncol(temp1))
  rownames(out) &lt;- pasteC(combs[,1],combs[,2])
  colnames(out) &lt;- 1:ncol(temp1)
  
  system.time({
  
  for( i in 1:nrow(combs)){
    out[i,] &lt;- combi_FUN_vec(    temp1[combs[i,1],]      ,         temp1[combs[i,2],]   )
  }
  
  })
  
  
    user  system elapsed 
     0.13    0.00    0.13
  </code></pre>
  
  <p><br><br></p>
  
  <pre><code>#B2. by apply
  class_computed &lt;- as.character(temp2[,1])
  temp2 &lt;- as.matrix(temp2[,-1])
  combs &lt;- t(combn(class_computed,2))
  rownames(temp2) &lt;- class_computed
  
  combi_FUN_vec &lt;- Vectorize(combi_FUN)
  
  system.time({
  
  out &lt;- apply(temp2,2,function(x){
    v &lt;- combi_FUN_vec(    x[combs[,1]]      ,        x[combs[,2]]    )
    names(v) &lt;- pasteC(combs[,1],combs[,2])
    return(v)
  })
  
  })
  
  
     user  system elapsed 
     0.91    0.00    0.91 
  </code></pre>
  
  <p><br><br></p>
  
  <p>Not surprisingly the loop is a clear winner and the output is this:</p>
  
  <pre><code>&gt; out[,1:5]
                1         2         3         4         5
  A - B 1.2128952 1.0161608 0.9313115 0.9767619 1.0132362
  A - C 1.0946079 1.0524154 0.9882857 0.9546686 0.9604382
  A - D 1.1872958 0.9113349 0.9941437 0.8751611 0.9863873
  A - E 1.1457396 0.9669100 0.9146375 0.8760513 1.0604971
  B - C 0.9024753 1.0356780 1.0611763 0.9773810 0.9478918
  B - D 0.9788940 0.8968413 1.0674664 0.8959820 0.9735018
  B - E 0.9446320 0.9515325 0.9820962 0.8968933 1.0466435
  C - D 1.0846768 0.8659461 1.0059275 0.9167172 1.0270179
  C - E 1.0467123 0.9187532 0.9254788 0.9176496 1.1041804
  D - E 0.9649993 1.0609821 0.9200254 1.0010171 1.0751326
  </code></pre>
  
  <p><br><br></p>
  
  <p>Now here is what bugs me... These two last methods are extremely ugly.
  Is there a better/cleaner/faster way? Preferably in <code>dplyr</code>/<code>purrr</code> style syntax? Perhaps without having to go through <code>combn</code> even?</p>
  
  <p>Any hints are appreciated.</p>
  
  <p><br><br></p>
  
  <p>EDIT:</p>
  
  <p>I manage to make a much more compact version in dplyr style:</p>
  
  <pre><code>f_un &lt;- unique(f)
    combs &lt;- t(combn(as.character(f_un),2))
  
    out3 &lt;- data.frame(class = f, MAT) %&gt;%  slice_rows("class") %&gt;% by_slice(map, aggr_FUN) %&gt;% 
            do(combi_FUN( slice(.,match(combs[,1], class))[,-1]  ,slice(.,match(combs[,2], class))[,-1]      )) %&gt;% 
            as.data.frame(row.names = pasteC(combs[,1],combs[,2]))
  </code></pre>
  
  <p>Is there a way to simplify that and speed it up? It is 10x slower than the fastest above.</p>
  
  <p><br><br>
  EDIT2:
  Based on the suggestions so far the fastest and cleanest is the below function.</p>
  
  <pre><code>fold.change &lt;- function(MAT,f,aggr_FUN=mean,combi_FUN=function(x,y) "/"(x,y)   ){
  
    # mean using purrr
    x &lt;- data.frame(class = f, MAT) %&gt;%  slice_rows("class") %&gt;% by_slice(map, aggr_FUN)
    rownames &lt;- as.character(as.data.frame(x[,1])[,1])
    x &lt;- as.matrix(x[,-1])
    rownames(x) &lt;- rownames
  
    # calculate changes between all rows
    i &lt;- combn(unique(f), 2)
    ret &lt;- combi_FUN(x[i[1,],] , x[i[2,],])
    rownames(ret) &lt;- pasteC(i[1,], i[2,])
  
    # Put original colnames
    colnames(ret) &lt;- colnames(MAT)
  
    return(ret)
  }
  </code></pre>