Statistics review 5: Comparison of means
========================================

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

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

-  The T test

.. code:: r

    suppressPackageStartupMessages(library(tidyverse))
    library(pwr)

.. code:: r

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

Comparison of a single mean with a hypothesized value
-----------------------------------------------------

.. code:: r

    hb <- c(8.1, 10.1, 12.3, 9.7, 11.7, 11.3,
            11.9, 9.3, 13.0, 10.5, 8.3, 8.8,
            9.4, 6.4, 5.4)

.. code:: r

    n <- length(hb)
    n



.. raw:: html

    15


.. code:: r

    mu <- round(mean(hb), 1)
    mu



.. raw:: html

    9.7


.. code:: r

    sigma <- round(sd(hb), 1)
    sigma



.. raw:: html

    2.2


.. code:: r

    se <- sigma/sqrt(n)
    se



.. raw:: html

    0.568037557443755


.. code:: r

    k <- qt(0.975, n-1)
    k



.. raw:: html

    2.1447866879178


.. code:: r

    ci <- c(mu - k*se, mu + k*se)
    round(ci, 1)



.. raw:: html

    <ol class=list-inline>
    	<li>8.5</li>
    	<li>10.9</li>
    </ol>



.. code:: r

    population_mean = 15.0

.. code:: r

    t = (mu - population_mean)/se
    round(t, 1)



.. raw:: html

    -9.3


Interpretation of T score
~~~~~~~~~~~~~~~~~~~~~~~~~

The observed mean is 9.3 standard errors below the hypothesized
population mean.

Getting a p value
~~~~~~~~~~~~~~~~~

.. code:: r

    2 * (1 - pt(abs(t), df = n-1))



.. raw:: html

    2.18878727586969e-07


In Practice: Using the t test
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    t.test(x = hb - population_mean)



.. parsed-literal::

    
    	One Sample t-test
    
    data:  hb - population_mean
    t = -9.4587, df = 14, p-value = 1.853e-07
    alternative hypothesis: true mean is not equal to 0
    95 percent confidence interval:
     -6.444537 -4.062130
    sample estimates:
    mean of x 
    -5.253333 



Comparison of two means arising from paired data
------------------------------------------------

.. code:: r

    before <- c(39.7, 59.1, 56.1, 57.7, 60.6, 37.8, 58.2, 33.6, 56.0, 65.3)
    after <- c(52.9, 56.7, 61.9, 71.4, 67.7, 50.0, 60.7, 51.3, 59.5, 59.8)
    
    df <- data.frame("Subject" = 1:length(before), "Before" = before, "After" = after)

.. code:: r

    df <- df %>% mutate(Difference = After - Before)
    df



.. raw:: html

    <table>
    <thead><tr><th scope=col>Subject</th><th scope=col>Before</th><th scope=col>After</th><th scope=col>Difference</th></tr></thead>
    <tbody>
    	<tr><td> 1  </td><td>39.7</td><td>52.9</td><td>13.2</td></tr>
    	<tr><td> 2  </td><td>59.1</td><td>56.7</td><td>-2.4</td></tr>
    	<tr><td> 3  </td><td>56.1</td><td>61.9</td><td> 5.8</td></tr>
    	<tr><td> 4  </td><td>57.7</td><td>71.4</td><td>13.7</td></tr>
    	<tr><td> 5  </td><td>60.6</td><td>67.7</td><td> 7.1</td></tr>
    	<tr><td> 6  </td><td>37.8</td><td>50.0</td><td>12.2</td></tr>
    	<tr><td> 7  </td><td>58.2</td><td>60.7</td><td> 2.5</td></tr>
    	<tr><td> 8  </td><td>33.6</td><td>51.3</td><td>17.7</td></tr>
    	<tr><td> 9  </td><td>56.0</td><td>59.5</td><td> 3.5</td></tr>
    	<tr><td>10  </td><td>65.3</td><td>59.8</td><td>-5.5</td></tr>
    </tbody>
    </table>



.. code:: r

    apply(df[,2:4], 2, mean)



.. raw:: html

    <dl class=dl-horizontal>
    	<dt>Before</dt>
    		<dd>52.41</dd>
    	<dt>After</dt>
    		<dd>59.19</dd>
    	<dt>Difference</dt>
    		<dd>6.78</dd>
    </dl>



.. code:: r

    mu <- mean(df$Difference)
    se <- (sd(df$Difference)/sqrt(10))
    t <- mu/se
    round(t, 2)



.. raw:: html

    2.87


.. code:: r

    p <- 2*(1 - pt(abs(t), 9))
    round(p, 5)



.. raw:: html

    0.01853


.. code:: r

    k <- qt(0.975, 9)
    round(k, 2)



.. raw:: html

    2.26


.. code:: r

    ci <- c(mu - k*se, mu + k*se)
    round(ci, 1)



.. raw:: html

    <ol class=list-inline>
    	<li>1.4</li>
    	<li>12.1</li>
    </ol>



The paired and one-sample T test are the same
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    t.test(df$After, df$Before, paired=TRUE)



.. parsed-literal::

    
    	Paired t-test
    
    data:  df$After and df$Before
    t = 2.8681, df = 9, p-value = 0.01853
    alternative hypothesis: true difference in means is not equal to 0
    95 percent confidence interval:
      1.432413 12.127587
    sample estimates:
    mean of the differences 
                       6.78 



.. code:: r

    t.test(df$Difference)



.. parsed-literal::

    
    	One Sample t-test
    
    data:  df$Difference
    t = 2.8681, df = 9, p-value = 0.01853
    alternative hypothesis: true mean is not equal to 0
    95 percent confidence interval:
      1.432413 12.127587
    sample estimates:
    mean of x 
         6.78 



Comparison of two means arising from unpaired data
--------------------------------------------------

.. code:: r

    n.x <- 119
    n.y <- 117
    mu.x <- 81
    mu.y <- 95
    sd.x <- 18
    sd.y <- 19

.. code:: r

    k <- qnorm(0.975)
    round(k, 2)



.. raw:: html

    1.96


.. code:: r

    ci.x <- c(mu.x - k*sd.x/sqrt(n.x), mu.x + k*sd.x/sqrt(n.x))
    round(ci.x, 1)



.. raw:: html

    <ol class=list-inline>
    	<li>77.8</li>
    	<li>84.2</li>
    </ol>



.. code:: r

    ci.y <- c(mu.y - k*sd.y/sqrt(n.y), mu.x + k*sd.y/sqrt(n.y))
    round(ci.y, 1)



.. raw:: html

    <ol class=list-inline>
    	<li>91.6</li>
    	<li>84.4</li>
    </ol>



.. code:: r

    pooled.se <- function(s1, s2, n1, n2) {
        sqrt(((n1-1)*s1^2 + (n2-1)*s2^2)/(n1+n2-2)) * sqrt(1/n1 + 1/n2)
    }

.. code:: r

    ps <- pooled.se(sd.x, sd.y, n.x, n.y)
    round(ps, 2)



.. raw:: html

    2.41


.. code:: r

    mu <- mu.x - mu.y
    mu



.. raw:: html

    -14


.. code:: r

    ci <- c(mu - k*ps, mu + k*se)
    round(ci, 1)



.. raw:: html

    <ol class=list-inline>
    	<li>-18.7</li>
    	<li>-9.4</li>
    </ol>



.. code:: r

    t <- mu/ps
    round(t, 2)



.. raw:: html

    -5.81


.. code:: r

    2 * (1 - pt(abs(t), df = n.x+n.y-2))



.. raw:: html

    2.01074956684977e-08


The two-sample T test
~~~~~~~~~~~~~~~~~~~~~

Normally you would just two vectors of data to the t.test function. We
will simulate these since the original data is not provided just to
illustrate.

.. code:: r

    x <- rnorm(n.x, mu.x, sd.x)
    y <- rnorm(n.y, mu.y, sd.y)
    t.test(x, y)



.. parsed-literal::

    
    	Welch Two Sample t-test
    
    data:  x and y
    t = -6.0677, df = 230.89, p-value = 5.268e-09
    alternative hypothesis: true difference in means is not equal to 0
    95 percent confidence interval:
     -21.26554 -10.84025
    sample estimates:
    mean of x mean of y 
     79.24579  95.29868 



Assumptions and limitations
---------------------------

-  The one sample t-test requires that the data have an approximately
   Normal distribution, whereas the paired t-test requires that the
   distribution of the differences are approximately Normal.
-  The unpaired t-test relies on the assumption that the data from the
   two samples are both Normally distributed, and has the additional
   requirement that the SDs from the two samples are approximately
   equal.

Testing for normality
~~~~~~~~~~~~~~~~~~~~~

This is most commonly done visually with a QQ-plot of percentiles of
sample versus percentiles of a standard normal distribution. If the
sample distribution is normal, you should see a straight diagonal line.

An alternative is the Shapiro-Wilk test for normality.

For testing equal variance, you can use the ``var.test`` function in R.

.. code:: r

    df <- data.frame(x=x)
    ggplot(df, aes(sample=x)) + stat_qq()





.. image:: SR05_Comparison_of_means_files/SR05_Comparison_of_means_44_1.png


.. code:: r

    shapiro.test(df$x)



.. parsed-literal::

    
    	Shapiro-Wilk normality test
    
    data:  df$x
    W = 0.99162, p-value = 0.6905



.. code:: r

    var.test(x, y)



.. parsed-literal::

    
    	F test to compare two variances
    
    data:  x and y
    F = 0.81959, num df = 118, denom df = 116, p-value = 0.283
    alternative hypothesis: true ratio of variances is not equal to 1
    95 percent confidence interval:
     0.5692599 1.1793311
    sample estimates:
    ratio of variances 
             0.8195921 



Exercise
--------

.. code:: r

    suppressPackageStartupMessages(library(MASS))

.. code:: r

    data(anorexia)

.. code:: r

    head(anorexia)



.. raw:: html

    <table>
    <thead><tr><th scope=col>Treat</th><th scope=col>Prewt</th><th scope=col>Postwt</th></tr></thead>
    <tbody>
    	<tr><td>Cont</td><td>80.7</td><td>80.2</td></tr>
    	<tr><td>Cont</td><td>89.4</td><td>80.1</td></tr>
    	<tr><td>Cont</td><td>91.8</td><td>86.4</td></tr>
    	<tr><td>Cont</td><td>74.0</td><td>86.3</td></tr>
    	<tr><td>Cont</td><td>78.1</td><td>76.1</td></tr>
    	<tr><td>Cont</td><td>88.3</td><td>78.1</td></tr>
    </tbody>
    </table>



.. code:: r

    str(anorexia)


.. parsed-literal::

    'data.frame':	72 obs. of  3 variables:
     $ Treat : Factor w/ 3 levels "CBT","Cont","FT": 2 2 2 2 2 2 2 2 2 2 ...
     $ Prewt : num  80.7 89.4 91.8 74 78.1 88.3 87.3 75.1 80.6 78.4 ...
     $ Postwt: num  80.2 80.1 86.4 86.3 76.1 78.1 75.1 86.7 73.5 84.6 ...


.. code:: r

    ggplot(anorexia %>% gather(Time, Wt, -Treat), aes(y=Wt, x=Treat, fill=Time)) + 
    geom_boxplot()





.. image:: SR05_Comparison_of_means_files/SR05_Comparison_of_means_52_1.png


**1**. There are 3 treatment groups for the anorexia data set. Calculate
the confidence intervals and p-values for each treatment using the
appropriate T test.


**2** Check if the normality assumptions are met for each treatment
group.