<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/compare-tube-types-with-r-repeated-measures-anova/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://labrtorian.com/2017/08/21/compare-tube-types-with-r-repeated-measures-anova/"> The Lab-R-torian</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Background</h2>
<p>Sometimes we might want to compare three or four tube types for a particular analyte on a group of patients or we might want to see if a particular analyte is stable over time in aliqioted samples. In these experiments are essentially doing the multivariable analogue of the paired t-test. In the tube-type experiment, the factor that is differing between the ('paired') groups is the container: serum separator tubes (SST), EDTA plasma tubes, plasma separator tubes (PST) etc. In a stability experiment, the factor that is differing is storage duration.</p>
<p>Since this is a fairly common clinical lab experiment, I thought I would just jot down how this is accomplished in R – though I must confess I know just about \(\lim_{x\to0}x\) about statistics. In any case, the statistical test is a repeated-measures ANOVA and this is one way to do it (there are many) including an approach to the post-hoc testing. </p>
<h2>Some Fake Data to Work With</h2>
<p>I'm going to make some fake data. I tried to dig up the data from an <a rel="nofollow" href="http://clinchem.aaccjnls.org/content/51/5/915.short" target="_blank">experiment I did as a resident</a> but alas, I think the raw data died on an old laptop. But fake data will do for demonstration purposes. Let's suppose we are looking at parathyroid hormone (PTH) in three different <a rel="nofollow" href="https://en.wikipedia.org/wiki/Vacutainer" target="_blank">vacutainer tubes</a>: SST, EDTA and PST. For the sake of argument, let's say that we collect samples from 20 patients simultaneously and we anlayze them all as per our usual process. This means that each patient has three samples of material that should be otherwise identical outside of the effects of the collection contained. </p>

<pre>library(magrittr)
set.seed(100) #to force the same pseudo-random each time
#data in pmol/L
#induce some heteroscedastic error
SST &lt;- runif(20,3,50)  
PST &lt;- 1.03*SST + rnorm(20,0,0.1)*SST #set the data up to show no difference
EDTA &lt;- 1.15*SST + rnorm(20,0,0.1)*SST  #set the data up to show a difference
tube.data &lt;- data.frame(SST,PST,EDTA) %>% round(.,1)
tube.data &lt;- data.frame(Subject = factor(1:20), tube.data)</pre>

<p>This is the way we usually express (and receive) data like this in an Excel spreadsheet:</p>
<table>
<thead>
<tr>
<th align="left">Subject</th>
<th align="right">SST</th>
<th align="right">PST</th>
<th align="right">EDTA</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left">1</td>
<td align="right">17.5</td>
<td align="right">18.1</td>
<td align="right">19.9</td>
</tr>
<tr>
<td align="left">2</td>
<td align="right">15.1</td>
<td align="right">15.7</td>
<td align="right">20.0</td>
</tr>
<tr>
<td align="left">3</td>
<td align="right">29.0</td>
<td align="right">29.2</td>
<td align="right">32.9</td>
</tr>
<tr>
<td align="left">4</td>
<td align="right">5.7</td>
<td align="right">6.2</td>
<td align="right">6.4</td>
</tr>
<tr>
<td align="left">5</td>
<td align="right">25.0</td>
<td align="right">26.1</td>
<td align="right">27.0</td>
</tr>
<tr>
<td align="left">6</td>
<td align="right">25.7</td>
<td align="right">26.4</td>
<td align="right">29.0</td>
</tr>
<tr>
<td align="left">7</td>
<td align="right">41.2</td>
<td align="right">40.8</td>
<td align="right">48.1</td>
</tr>
<tr>
<td align="left">8</td>
<td align="right">20.4</td>
<td align="right">22.1</td>
<td align="right">24.3</td>
</tr>
<tr>
<td align="left">9</td>
<td align="right">28.7</td>
<td align="right">26.9</td>
<td align="right">36.0</td>
</tr>
<tr>
<td align="left">10</td>
<td align="right">11.0</td>
<td align="right">13.9</td>
<td align="right">13.7</td>
</tr>
<tr>
<td align="left">11</td>
<td align="right">32.4</td>
<td align="right">31.9</td>
<td align="right">36.9</td>
</tr>
<tr>
<td align="left">12</td>
<td align="right">44.5</td>
<td align="right">49.2</td>
<td align="right">57.4</td>
</tr>
<tr>
<td align="left">13</td>
<td align="right">16.2</td>
<td align="right">17.1</td>
<td align="right">15.7</td>
</tr>
<tr>
<td align="left">14</td>
<td align="right">21.7</td>
<td align="right">24.1</td>
<td align="right">26.3</td>
</tr>
<tr>
<td align="left">15</td>
<td align="right">38.8</td>
<td align="right">36.8</td>
<td align="right">42.6</td>
</tr>
<tr>
<td align="left">16</td>
<td align="right">34.4</td>
<td align="right">34.0</td>
<td align="right">44.2</td>
</tr>
<tr>
<td align="left">17</td>
<td align="right">12.6</td>
<td align="right">12.1</td>
<td align="right">14.1</td>
</tr>
<tr>
<td align="left">18</td>
<td align="right">19.8</td>
<td align="right">20.9</td>
<td align="right">25.4</td>
</tr>
<tr>
<td align="left">19</td>
<td align="right">19.9</td>
<td align="right">18.2</td>
<td align="right">23.0</td>
</tr>
<tr>
<td align="left">20</td>
<td align="right">35.4</td>
<td align="right">37.4</td>
<td align="right">34.1</td>
</tr>
</tbody>
</table>
<p>This Excel-ish way of storing the data is referred to as the “datawide” format for obvious reasons.</p>
<h2>Gather the Grain</h2>
<p>As it turns out this is not the way that we want to store data to do the statistical analyses of interest. What we want to do is have the tube type in a single column because this is the factor that is different within the subjects. We want to <code>gather()</code> or <code>melt()</code> the data (depending on your package of choice) to be like so:</p>

<pre>library(tidyr)
tube.data.2 &lt;- gather(tube.data, key = &quot;Subject&quot;)
tube.data.2 %>% kable()</pre>

<table>
<thead>
<tr>
<th align="left">Subject</th>
<th align="left">Subject</th>
<th align="right">value</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left">1</td>
<td align="left">SST</td>
<td align="right">17.5</td>
</tr>
<tr>
<td align="left">2</td>
<td align="left">SST</td>
<td align="right">15.1</td>
</tr>
<tr>
<td align="left">3</td>
<td align="left">SST</td>
<td align="right">29.0</td>
</tr>
<tr>
<td align="left">4</td>
<td align="left">SST</td>
<td align="right">5.7</td>
</tr>
<tr>
<td align="left">5</td>
<td align="left">SST</td>
<td align="right">25.0</td>
</tr>
<tr>
<td align="left">6</td>
<td align="left">SST</td>
<td align="right">25.7</td>
</tr>
<tr>
<td align="left">7</td>
<td align="left">SST</td>
<td align="right">41.2</td>
</tr>
<tr>
<td align="left">8</td>
<td align="left">SST</td>
<td align="right">20.4</td>
</tr>
<tr>
<td align="left">9</td>
<td align="left">SST</td>
<td align="right">28.7</td>
</tr>
<tr>
<td align="left">10</td>
<td align="left">SST</td>
<td align="right">11.0</td>
</tr>
<tr>
<td align="left">11</td>
<td align="left">SST</td>
<td align="right">32.4</td>
</tr>
<tr>
<td align="left">12</td>
<td align="left">SST</td>
<td align="right">44.5</td>
</tr>
<tr>
<td align="left">13</td>
<td align="left">SST</td>
<td align="right">16.2</td>
</tr>
<tr>
<td align="left">14</td>
<td align="left">SST</td>
<td align="right">21.7</td>
</tr>
<tr>
<td align="left">15</td>
<td align="left">SST</td>
<td align="right">38.8</td>
</tr>
<tr>
<td align="left">16</td>
<td align="left">SST</td>
<td align="right">34.4</td>
</tr>
<tr>
<td align="left">17</td>
<td align="left">SST</td>
<td align="right">12.6</td>
</tr>
<tr>
<td align="left">18</td>
<td align="left">SST</td>
<td align="right">19.8</td>
</tr>
<tr>
<td align="left">19</td>
<td align="left">SST</td>
<td align="right">19.9</td>
</tr>
<tr>
<td align="left">20</td>
<td align="left">SST</td>
<td align="right">35.4</td>
</tr>
<tr>
<td align="left">1</td>
<td align="left">PST</td>
<td align="right">18.1</td>
</tr>
<tr>
<td align="left">2</td>
<td align="left">PST</td>
<td align="right">15.7</td>
</tr>
<tr>
<td align="left">3</td>
<td align="left">PST</td>
<td align="right">29.2</td>
</tr>
<tr>
<td align="left">4</td>
<td align="left">PST</td>
<td align="right">6.2</td>
</tr>
<tr>
<td align="left">5</td>
<td align="left">PST</td>
<td align="right">26.1</td>
</tr>
<tr>
<td align="left">6</td>
<td align="left">PST</td>
<td align="right">26.4</td>
</tr>
<tr>
<td align="left">7</td>
<td align="left">PST</td>
<td align="right">40.8</td>
</tr>
<tr>
<td align="left">8</td>
<td align="left">PST</td>
<td align="right">22.1</td>
</tr>
<tr>
<td align="left">9</td>
<td align="left">PST</td>
<td align="right">26.9</td>
</tr>
<tr>
<td align="left">10</td>
<td align="left">PST</td>
<td align="right">13.9</td>
</tr>
<tr>
<td align="left">11</td>
<td align="left">PST</td>
<td align="right">31.9</td>
</tr>
<tr>
<td align="left">12</td>
<td align="left">PST</td>
<td align="right">49.2</td>
</tr>
<tr>
<td align="left">13</td>
<td align="left">PST</td>
<td align="right">17.1</td>
</tr>
<tr>
<td align="left">14</td>
<td align="left">PST</td>
<td align="right">24.1</td>
</tr>
<tr>
<td align="left">15</td>
<td align="left">PST</td>
<td align="right">36.8</td>
</tr>
<tr>
<td align="left">16</td>
<td align="left">PST</td>
<td align="right">34.0</td>
</tr>
<tr>
<td align="left">17</td>
<td align="left">PST</td>
<td align="right">12.1</td>
</tr>
<tr>
<td align="left">18</td>
<td align="left">PST</td>
<td align="right">20.9</td>
</tr>
<tr>
<td align="left">19</td>
<td align="left">PST</td>
<td align="right">18.2</td>
</tr>
<tr>
<td align="left">20</td>
<td align="left">PST</td>
<td align="right">37.4</td>
</tr>
<tr>
<td align="left">1</td>
<td align="left">EDTA</td>
<td align="right">19.9</td>
</tr>
<tr>
<td align="left">2</td>
<td align="left">EDTA</td>
<td align="right">20.0</td>
</tr>
<tr>
<td align="left">3</td>
<td align="left">EDTA</td>
<td align="right">32.9</td>
</tr>
<tr>
<td align="left">4</td>
<td align="left">EDTA</td>
<td align="right">6.4</td>
</tr>
<tr>
<td align="left">5</td>
<td align="left">EDTA</td>
<td align="right">27.0</td>
</tr>
<tr>
<td align="left">6</td>
<td align="left">EDTA</td>
<td align="right">29.0</td>
</tr>
<tr>
<td align="left">7</td>
<td align="left">EDTA</td>
<td align="right">48.1</td>
</tr>
<tr>
<td align="left">8</td>
<td align="left">EDTA</td>
<td align="right">24.3</td>
</tr>
<tr>
<td align="left">9</td>
<td align="left">EDTA</td>
<td align="right">36.0</td>
</tr>
<tr>
<td align="left">10</td>
<td align="left">EDTA</td>
<td align="right">13.7</td>
</tr>
<tr>
<td align="left">11</td>
<td align="left">EDTA</td>
<td align="right">36.9</td>
</tr>
<tr>
<td align="left">12</td>
<td align="left">EDTA</td>
<td align="right">57.4</td>
</tr>
<tr>
<td align="left">13</td>
<td align="left">EDTA</td>
<td align="right">15.7</td>
</tr>
<tr>
<td align="left">14</td>
<td align="left">EDTA</td>
<td align="right">26.3</td>
</tr>
<tr>
<td align="left">15</td>
<td align="left">EDTA</td>
<td align="right">42.6</td>
</tr>
<tr>
<td align="left">16</td>
<td align="left">EDTA</td>
<td align="right">44.2</td>
</tr>
<tr>
<td align="left">17</td>
<td align="left">EDTA</td>
<td align="right">14.1</td>
</tr>
<tr>
<td align="left">18</td>
<td align="left">EDTA</td>
<td align="right">25.4</td>
</tr>
<tr>
<td align="left">19</td>
<td align="left">EDTA</td>
<td align="right">23.0</td>
</tr>
<tr>
<td align="left">20</td>
<td align="left">EDTA</td>
<td align="right">34.1</td>
</tr>
</tbody>
</table>
<p>Now we see that there is a column for tube-type and a column for the PTH results which we can name accordingly. You can see why this called the “datalong” format.</p>

<pre>names(tube.data.2) &lt;- c(&quot;Subject&quot;, &quot;Tube.Type&quot;, &quot;PTH&quot;)
tube.data.2$Tube.Type &lt;- as.factor(tube.data.2$Tube.Type) #turns tube type into factor</pre>

<h2>Visualize</h2>
<p>Summarize the data:</p>

<pre>summary(tube.data)</pre>

<pre>##     Subject        SST             PST             EDTA      
##  1      : 1   Min.   : 5.70   Min.   : 6.20   Min.   : 6.40  
##  2      : 1   1st Qu.:17.18   1st Qu.:17.85   1st Qu.:19.98  
##  3      : 1   Median :23.35   Median :25.10   Median :26.65  
##  4      : 1   Mean   :24.75   Mean   :25.36   Mean   :28.85  
##  5      : 1   3rd Qu.:32.90   3rd Qu.:32.42   3rd Qu.:36.23  
##  6      : 1   Max.   :44.50   Max.   :49.20   Max.   :57.40  
##  (Other):14</pre>

<p>Let's just have a quick look graphically:</p>

<pre>library(mcr)
plot(mcreg(SST, EDTA,
           method.reg = &quot;PaBa&quot;,
           mref.name = &quot;SST&quot;,
           mtest.name = &quot;EDTA&quot;))</pre>

<p><img alt="plot of chunk unnamed-chunk-6" src="https://i1.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-6-1.png?w=604" data-recalc-dims="1"></p>

<pre>plot(mcreg(SST, PST,
           method.reg = &quot;PaBa&quot;,
           mref.name = &quot;SST&quot;,
           mtest.name = &quot;PST&quot;))</pre>

<p><img alt="plot of chunk unnamed-chunk-6" src="https://i2.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-6-2.png?w=604" data-recalc-dims="1"></p>
<p>And as a boxplot with the points overtop:</p>

<pre>boxplot(PTH ~ Tube.Type,
        data = tube.data.2,
        col = c(&quot;purple&quot;, &quot;lightgreen&quot;, &quot;gold&quot;))
stripchart(PTH ~ Tube.Type,
           vertical = TRUE,
           data = tube.data.2, 
           method = &quot;jitter&quot;,
           add = TRUE,
           pch = 20,
           col = rgb(0,0,0,0.5))</pre>

<p><img alt="plot of chunk unnamed-chunk-7" src="https://i1.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-7-1.png?w=604" data-recalc-dims="1"></p>
<h2>Separate the Wheat from the Chaff</h2>
<p>Now we want to make comparisons to see if these are different. To accomplish this, we will use the <code>aov()</code> function. This requires us to have data formatted “datalong” as it is in the <code>tube.data.2</code> dataframe.</p>

<pre>fit &lt;- aov(PTH ~ Tube.Type + Error(Subject/Tube.Type), data=tube.data.2)</pre>

<p>If you are like me, this syntax is confusing. But it goes like this. <code>PTH</code> is a function of <code>Tube.Type</code> which is straight forward–hence the <code>PTH ~ Tube.Type</code> bit. The error term has the <code>Subject</code> in front of the <code>/</code> and the factor that is different within the subjects (<code>Tube.Type</code>) after the <code>/</code>. That's my grade 2 explanation from reading <a rel="nofollow" href="http://personality-project.org/r/r.guide/r.anova.html" target="_blank">this</a> and <a rel="nofollow" href="http://coltekin.net/cagri/R/r-exercisesse5.html" target="_blank">this</a> and <a rel="nofollow" href="http://www.psych.upenn.edu/~baron/rpsych/rpsych.html#htoc60" target="_blank">this</a>. </p>

<pre>summary(fit)</pre>

<pre>## 
## Error: Subject
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 19   7307   384.6               
## 
## Error: Subject:Tube.Type
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Tube.Type  2  195.9   97.97   22.47 3.63e-07 ***
## Residuals 38  165.7    4.36                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</pre>

<p>This tells us that there is a difference between the groups but it does not specify where the difference is.</p>
<h2>I can't see the difference. Can you see the difference?</h2>
<p>Sorry – I just had to make a pop-culture reference to <a rel="nofollow" href="https://www.youtube.com/watch?v=LHyNCGSYB4s" target="_blank">this</a>. We want to be specific about where the differences are without making a <a rel="nofollow" href="https://en.wikipedia.org/wiki/Type_I_and_type_II_errors" target="_blank">Type I error</a> which might arise if we blindly charge ahead and do multiple paired t-tests. One easy way to accomplish this is to use the <code>pairwise.t.test()</code> function which does corrections for multiple comparisons. You can choose from a number of approaches for adjustment for pairwise comparison. This requires the “response vector” which is PTH and the “grouping factor” which is the tube type.</p>

<pre># choices for p.adjust.method are: c(&quot;holm&quot;, &quot;hochberg&quot;, &quot;hommel&quot;, &quot;bonferroni&quot;, &quot;BH&quot;, &quot;BY&quot;, &quot;fdr&quot;, &quot;none&quot;)
pwt &lt;- pairwise.t.test(tube.data.2$PTH, tube.data.2$Tube.Type, p.adj = &quot;bonferroni&quot;, paired = TRUE)
pwt</pre>

<pre>## 
##  Pairwise comparisons using paired t tests 
## 
## data:  tube.data.2$PTH and tube.data.2$Tube.Type 
## 
##     EDTA    PST    
## PST 0.00083 -      
## SST 7.9e-05 0.35033
## 
## P value adjustment method: bonferroni</pre>

<p>This is pretty easy to understand. There are statistically significant differences found between the EDTA and PST (p = 0.00083) and the EDTA and PST (p = 0.00008) but none between SST and PST (p = 0.35033).</p>
<h2>Conclusion</h2>
<p>Non-statistician's approach to tube-type comparisons, which is also applicable to analyte stability studies. This is a one-way repeated measures ANOVA with one within-subjects factor. There is a great deal more to say on the matter by people who know much more in the citations in the links provided above.</p>

<p><strong>God probably uses datawide format</strong></p>
<p><em>All the nations will be gathered before him, and he will separate the people one from another as a shepherd separates the sheep from the goats. He will put the sheep on his right and the goats on his left.</em></p>
<p>(Matt 25:32-33)</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://labrtorian.com/2017/08/21/compare-tube-types-with-r-repeated-measures-anova/"> The Lab-R-torian</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/P6IMj3ruWms" height="1"></a></div>