Statistics review 10: Further nonparametric methods
===================================================

R code accompanying
`paper <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC420045/pdf/cc2836.pdf>`__

Key learning points
-------------------

-  Nonparametric methods for testing differences between more than two
   groups or treatments

.. code:: r

    suppressPackageStartupMessages(library(tidyverse))

.. code:: r

    options(repr.plot.width=4, repr.plot.height=3)

Kruskal–Wallis test
-------------------

Nonparametric alternative to one-way analysis of variance.

Manual calculation
~~~~~~~~~~~~~~~~~~

.. code:: r

    ct <- c(7,1,2,6,11,8)
    m <- c(4,7,16,11,21)
    ns <- c(20,25,13,9,14,11)

.. code:: r

    icu <- c(rep("ct", length(ct)), rep("m", length(m)), rep("ns", length(ns)))
    days <- c(ct, m, ns)
    
    df <- data.frame(icu=icu, days=days)

.. code:: r

    head(df)



.. raw:: html

    <table>
    <thead><tr><th scope=col>icu</th><th scope=col>days</th></tr></thead>
    <tbody>
    	<tr><td>ct</td><td> 7</td></tr>
    	<tr><td>ct</td><td> 1</td></tr>
    	<tr><td>ct</td><td> 2</td></tr>
    	<tr><td>ct</td><td> 6</td></tr>
    	<tr><td>ct</td><td>11</td></tr>
    	<tr><td>ct</td><td> 8</td></tr>
    </tbody>
    </table>



.. code:: r

    ggplot(df, aes(x=icu, y=days)) + geom_boxplot()





.. image:: SR10_Further_nonparametric_methods_files/SR10_Further_nonparametric_methods_8_1.png


.. code:: r

    df1 <- df %>% mutate(rank=rank(days))
    head(df1)



.. raw:: html

    <table>
    <thead><tr><th scope=col>icu</th><th scope=col>days</th><th scope=col>rank</th></tr></thead>
    <tbody>
    	<tr><td>ct  </td><td> 7  </td><td> 5.5</td></tr>
    	<tr><td>ct  </td><td> 1  </td><td> 1.0</td></tr>
    	<tr><td>ct  </td><td> 2  </td><td> 2.0</td></tr>
    	<tr><td>ct  </td><td> 6  </td><td> 4.0</td></tr>
    	<tr><td>ct  </td><td>11  </td><td>10.0</td></tr>
    	<tr><td>ct  </td><td> 8  </td><td> 7.0</td></tr>
    </tbody>
    </table>



.. code:: r

    df2 <- df1 %>% group_by(icu) %>% summarize(r = sum(rank), n=n())
    df2



.. raw:: html

    <table>
    <thead><tr><th scope=col>icu</th><th scope=col>r</th><th scope=col>n</th></tr></thead>
    <tbody>
    	<tr><td>ct  </td><td>29.5</td><td>6   </td></tr>
    	<tr><td>m   </td><td>48.5</td><td>5   </td></tr>
    	<tr><td>ns  </td><td>75.0</td><td>6   </td></tr>
    </tbody>
    </table>



.. code:: r

    N <- dim(df1)[1]
    
    T <- (12/(N*(N+1)))*sum((df2$r^2/df2$n)) - 3*(N+1)

.. code:: r

    round(T, 2)



.. raw:: html

    6.9


.. code:: r

    df <- dim(df2)[1]- 1
    df



.. raw:: html

    2


.. code:: r

    round(1 - pchisq(T, df), 4)



.. raw:: html

    0.0317


Using built in function
~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    kruskal.test(days ~ icu, data=df1)



.. parsed-literal::

    
    	Kruskal-Wallis rank sum test
    
    data:  days by icu
    Kruskal-Wallis chi-squared = 6.9442, df = 2, p-value = 0.03105



Dunn test for multiple comparisons
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    suppressPackageStartupMessages(library(FSA))

.. code:: r

    dunnTest(df1$days, df1$icu)


.. parsed-literal::

    Dunn (1964) Kruskal-Wallis multiple comparison
      p-values adjusted with the Holm method.
    



.. parsed-literal::

      Comparison          Z    P.unadj      P.adj
    1     ct - m -1.5691321 0.11661717 0.23323435
    2    ct - ns -2.6090676 0.00907893 0.02723679
    3     m - ns -0.9185163 0.35834862 0.35834862


Ad-hoc comparisons with Wilcoxon
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    pairwise.wilcox.test(df1$days, df1$icu, p.adjust.method = "holm")


.. parsed-literal::

    Warning message in wilcox.test.default(xi, xj, paired = paired, ...):
    “cannot compute exact p-value with ties”Warning message in wilcox.test.default(xi, xj, paired = paired, ...):
    “cannot compute exact p-value with ties”Warning message in wilcox.test.default(xi, xj, paired = paired, ...):
    “cannot compute exact p-value with ties”


.. parsed-literal::

    
    	Pairwise comparisons using Wilcoxon rank sum test 
    
    data:  df1$days and df1$icu 
    
       ct    m    
    m  0.338 -    
    ns 0.031 0.464
    
    P value adjustment method: holm 


The Jonckheere–Terpstra test
----------------------------

For comparisons where the treatment group is ordinal.

We re-use the same data set, assuming that the ICU ordering is ct, m,
ns.

Using a package
~~~~~~~~~~~~~~~

.. code:: r

    head(df1)



.. raw:: html

    <table>
    <thead><tr><th scope=col>icu</th><th scope=col>days</th><th scope=col>rank</th></tr></thead>
    <tbody>
    	<tr><td>ct  </td><td> 7  </td><td> 5.5</td></tr>
    	<tr><td>ct  </td><td> 1  </td><td> 1.0</td></tr>
    	<tr><td>ct  </td><td> 2  </td><td> 2.0</td></tr>
    	<tr><td>ct  </td><td> 6  </td><td> 4.0</td></tr>
    	<tr><td>ct  </td><td>11  </td><td>10.0</td></tr>
    	<tr><td>ct  </td><td> 8  </td><td> 7.0</td></tr>
    </tbody>
    </table>



.. code:: r

    df1$icu <- factor(df1$icu, levels=c("ct", "m","ns"))

.. code:: r

    str(df1)


.. parsed-literal::

    'data.frame':	17 obs. of  3 variables:
     $ icu : Factor w/ 3 levels "ct","m","ns": 1 1 1 1 1 1 2 2 2 2 ...
     $ days: num  7 1 2 6 11 8 4 7 16 11 ...
     $ rank: num  5.5 1 2 4 10 7 3 5.5 14 10 ...


.. code:: r

    library(clinfun)

.. code:: r

    jonckheere.test(df1$days, as.numeric(df1$icu), nperm=10000)



.. parsed-literal::

    
    	Jonckheere-Terpstra test
    
    data:  
    JT = 77, p-value = 0.0102
    alternative hypothesis: two.sided



.. code:: r

    df1$days



.. raw:: html

    <ol class=list-inline>
    	<li>7</li>
    	<li>1</li>
    	<li>2</li>
    	<li>6</li>
    	<li>11</li>
    	<li>8</li>
    	<li>4</li>
    	<li>7</li>
    	<li>16</li>
    	<li>11</li>
    	<li>21</li>
    	<li>20</li>
    	<li>25</li>
    	<li>13</li>
    	<li>9</li>
    	<li>14</li>
    	<li>11</li>
    </ol>



.. code:: r

    as.numeric(df1$icu)



.. raw:: html

    <ol class=list-inline>
    	<li>1</li>
    	<li>1</li>
    	<li>1</li>
    	<li>1</li>
    	<li>1</li>
    	<li>1</li>
    	<li>2</li>
    	<li>2</li>
    	<li>2</li>
    	<li>2</li>
    	<li>2</li>
    	<li>3</li>
    	<li>3</li>
    	<li>3</li>
    	<li>3</li>
    	<li>3</li>
    	<li>3</li>
    </ol>



The Friedman Test
-----------------

Extension of the sign test for matched pairs and is used when the data
arise from more than two related samples. The Friedman test is the
non-parametric alternative to the one-way ANOVA with repeated measures.
Used as a two-way ANOVA with a completely balanced design.

.. code:: r

    A <- c(6,9,10,14,11)
    B <- c(9,16,14,14,16)
    C <- c(10,16,22,40,17)
    D <- c(16,32,67,19,60)
    
    df3 <- data.frame(Patient = 1:5, A=A, B=B, C=C, D=D)

.. code:: r

    df3.rank <- t(apply(df3[,2:5], 1, rank))
    df3.rank <- data.frame(df3.rank)
    df3.rank



.. raw:: html

    <table>
    <thead><tr><th scope=col>A</th><th scope=col>B</th><th scope=col>C</th><th scope=col>D</th></tr></thead>
    <tbody>
    	<tr><td>1.0</td><td>2.0</td><td>3.0</td><td>4  </td></tr>
    	<tr><td>1.0</td><td>2.5</td><td>2.5</td><td>4  </td></tr>
    	<tr><td>1.0</td><td>2.0</td><td>3.0</td><td>4  </td></tr>
    	<tr><td>1.5</td><td>1.5</td><td>4.0</td><td>3  </td></tr>
    	<tr><td>1.0</td><td>2.0</td><td>3.0</td><td>4  </td></tr>
    </tbody>
    </table>



.. code:: r

    df4 <- df3.rank %>% summarise_each("sum")
    df4



.. raw:: html

    <table>
    <thead><tr><th scope=col>A</th><th scope=col>B</th><th scope=col>C</th><th scope=col>D</th></tr></thead>
    <tbody>
    	<tr><td>5.5 </td><td>10  </td><td>15.5</td><td>19  </td></tr>
    </tbody>
    </table>



.. code:: r

    k <- dim(df3.rank)[2] 
    b <- dim(df3.rank)[1]
    
    T <- (12/(b*k*(k+1)))*sum(df4^2) - 3*b*(k+1)

.. code:: r

    T



.. raw:: html

    12.78


.. code:: r

    round(1 - pchisq(T, k-1), 4)



.. raw:: html

    0.0051


.. code:: r

    df3



.. raw:: html

    <table>
    <thead><tr><th scope=col>Patient</th><th scope=col>A</th><th scope=col>B</th><th scope=col>C</th><th scope=col>D</th></tr></thead>
    <tbody>
    	<tr><td>1 </td><td> 6</td><td> 9</td><td>10</td><td>16</td></tr>
    	<tr><td>2 </td><td> 9</td><td>16</td><td>16</td><td>32</td></tr>
    	<tr><td>3 </td><td>10</td><td>14</td><td>22</td><td>67</td></tr>
    	<tr><td>4 </td><td>14</td><td>14</td><td>40</td><td>19</td></tr>
    	<tr><td>5 </td><td>11</td><td>16</td><td>17</td><td>60</td></tr>
    </tbody>
    </table>



.. code:: r

    df5 <- df3 %>% gather(treatment, score, A:D)
    head(df5)



.. raw:: html

    <table>
    <thead><tr><th scope=col>Patient</th><th scope=col>treatment</th><th scope=col>score</th></tr></thead>
    <tbody>
    	<tr><td>1 </td><td>A </td><td> 6</td></tr>
    	<tr><td>2 </td><td>A </td><td> 9</td></tr>
    	<tr><td>3 </td><td>A </td><td>10</td></tr>
    	<tr><td>4 </td><td>A </td><td>14</td></tr>
    	<tr><td>5 </td><td>A </td><td>11</td></tr>
    	<tr><td>1 </td><td>B </td><td> 9</td></tr>
    </tbody>
    </table>



.. code:: r

    str(df5)


.. parsed-literal::

    'data.frame':	20 obs. of  3 variables:
     $ Patient  : int  1 2 3 4 5 1 2 3 4 5 ...
     $ treatment: chr  "A" "A" "A" "A" ...
     $ score    : num  6 9 10 14 11 9 16 14 14 16 ...


First way to call Friedman test
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: r

    friedman.test(score ~ treatment | Patient, data=df5)



.. parsed-literal::

    
    	Friedman rank sum test
    
    data:  score and treatment and Patient
    Friedman chi-squared = 13.312, df = 3, p-value = 0.004007



Second way to call Friedman test
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: r

    friedman.test(df5$score, df5$treatment, df5$Patient)



.. parsed-literal::

    
    	Friedman rank sum test
    
    data:  df5$score, df5$treatment and df5$Patient
    Friedman chi-squared = 13.312, df = 3, p-value = 0.004007



.. code:: r

    head(chickwts)



.. raw:: html

    <table>
    <thead><tr><th scope=col>weight</th><th scope=col>feed</th></tr></thead>
    <tbody>
    	<tr><td>179      </td><td>horsebean</td></tr>
    	<tr><td>160      </td><td>horsebean</td></tr>
    	<tr><td>136      </td><td>horsebean</td></tr>
    	<tr><td>227      </td><td>horsebean</td></tr>
    	<tr><td>217      </td><td>horsebean</td></tr>
    	<tr><td>168      </td><td>horsebean</td></tr>
    </tbody>
    </table>



Exercise
--------

**1**. Practice using the non-parametric tests on the ``chickwts`` data
set.