Statistics review 1: Presenting and summarising data
====================================================

R code to accompany
`paper <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC137399/pdf/cc1455.pdf>`__

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

-  The first step in any analysis is to describe and summarize the data
-  Understand the following

   -  qualitative data (unordered and ordered) and quantitative data
      (discrete and continuous)
   -  how these types of data can be represented figuratively
   -  the two important features of a quantitative dataset (location and
      variability)
   -  the measures of location (mean, median and mode)
   -  the measures of variability (range, interquartile range, standard
      deviation and variance)
   -  common distributions of clinical data
   -  simple transformations of positively skewed data.

.. code:: r

    suppressPackageStartupMessages(library(tidyverse))

Set default plot size to 4" by 3"
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: r

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

Summary Statistics
------------------

Read in data set from file
~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    df1 <- read.csv('data/haemoglobin.csv', header=FALSE, col.names=c("hb"))

Working with data frames
^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: r

    class(df1)



.. raw:: html

    'data.frame'


.. code:: r

    head(df1)



.. raw:: html

    <table>
    <thead><tr><th scope=col>hb</th></tr></thead>
    <tbody>
    	<tr><td>5.4</td></tr>
    	<tr><td>8.2</td></tr>
    	<tr><td>6.4</td></tr>
    	<tr><td>8.3</td></tr>
    	<tr><td>6.4</td></tr>
    	<tr><td>8.3</td></tr>
    </tbody>
    </table>



.. code:: r

    str(df1)


.. parsed-literal::

    'data.frame':	48 obs. of  1 variable:
     $ hb: num  5.4 8.2 6.4 8.3 6.4 8.3 7 8.6 7.1 8.8 ...


.. code:: r

    dim(df1)



.. raw:: html

    <ol class=list-inline>
    	<li>48</li>
    	<li>1</li>
    </ol>



.. code:: r

    summary(df1)



.. parsed-literal::

           hb        
     Min.   : 5.400  
     1st Qu.: 8.750  
     Median : 9.800  
     Mean   : 9.869  
     3rd Qu.:10.800  
     Max.   :14.100  


Dataframe indexing
^^^^^^^^^^^^^^^^^^

.. code:: r

    df1$hb



.. raw:: html

    <ol class=list-inline>
    	<li>5.4</li>
    	<li>8.2</li>
    	<li>6.4</li>
    	<li>8.3</li>
    	<li>6.4</li>
    	<li>8.3</li>
    	<li>7</li>
    	<li>8.6</li>
    	<li>7.1</li>
    	<li>8.8</li>
    	<li>7.3</li>
    	<li>8.9</li>
    	<li>7.7</li>
    	<li>9.1</li>
    	<li>8.1</li>
    	<li>9.3</li>
    	<li>9.3</li>
    	<li>9.9</li>
    	<li>9.4</li>
    	<li>9.9</li>
    	<li>9.4</li>
    	<li>9.9</li>
    	<li>9.4</li>
    	<li>10.1</li>
    	<li>9.4</li>
    	<li>10.3</li>
    	<li>9.5</li>
    	<li>10.3</li>
    	<li>9.7</li>
    	<li>10.4</li>
    	<li>9.7</li>
    	<li>10.4</li>
    	<li>10.5</li>
    	<li>11.9</li>
    	<li>10.5</li>
    	<li>12.3</li>
    	<li>10.6</li>
    	<li>12.6</li>
    	<li>10.8</li>
    	<li>12.7</li>
    	<li>10.8</li>
    	<li>13</li>
    	<li>11.3</li>
    	<li>13.3</li>
    	<li>11.7</li>
    	<li>14</li>
    	<li>11.7</li>
    	<li>14.1</li>
    </ol>



.. code:: r

    df1[,1]



.. raw:: html

    <ol class=list-inline>
    	<li>5.4</li>
    	<li>8.2</li>
    	<li>6.4</li>
    	<li>8.3</li>
    	<li>6.4</li>
    	<li>8.3</li>
    	<li>7</li>
    	<li>8.6</li>
    	<li>7.1</li>
    	<li>8.8</li>
    	<li>7.3</li>
    	<li>8.9</li>
    	<li>7.7</li>
    	<li>9.1</li>
    	<li>8.1</li>
    	<li>9.3</li>
    	<li>9.3</li>
    	<li>9.9</li>
    	<li>9.4</li>
    	<li>9.9</li>
    	<li>9.4</li>
    	<li>9.9</li>
    	<li>9.4</li>
    	<li>10.1</li>
    	<li>9.4</li>
    	<li>10.3</li>
    	<li>9.5</li>
    	<li>10.3</li>
    	<li>9.7</li>
    	<li>10.4</li>
    	<li>9.7</li>
    	<li>10.4</li>
    	<li>10.5</li>
    	<li>11.9</li>
    	<li>10.5</li>
    	<li>12.3</li>
    	<li>10.6</li>
    	<li>12.6</li>
    	<li>10.8</li>
    	<li>12.7</li>
    	<li>10.8</li>
    	<li>13</li>
    	<li>11.3</li>
    	<li>13.3</li>
    	<li>11.7</li>
    	<li>14</li>
    	<li>11.7</li>
    	<li>14.1</li>
    </ol>



.. code:: r

    df1[1,]



.. raw:: html

    5.4


.. code:: r

    df1[5:10,]



.. raw:: html

    <ol class=list-inline>
    	<li>6.4</li>
    	<li>8.3</li>
    	<li>7</li>
    	<li>8.6</li>
    	<li>7.1</li>
    	<li>8.8</li>
    </ol>



Measuring location
~~~~~~~~~~~~~~~~~~

Using a custom function
^^^^^^^^^^^^^^^^^^^^^^^

.. code:: r

    sr.mean <- function(x) {
        sum(x)/length(x)
    }

.. code:: r

    sr.mean(df1$hb)



.. raw:: html

    9.86875


Using built-in functions
^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: r

    df1 %>% summarize(mean=mean(hb), median=median(hb), mode=mode(hb))



.. raw:: html

    <table>
    <thead><tr><th scope=col>mean</th><th scope=col>median</th><th scope=col>mode</th></tr></thead>
    <tbody>
    	<tr><td>9.86875</td><td>9.8    </td><td>numeric</td></tr>
    </tbody>
    </table>



Visualizing the data distribution
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: r

    ggplot(df1, aes(x=hb)) +
    geom_histogram(binwidth=0.5, fill="grey", color="darkgrey")





.. image:: SR01_Presenting_and_summarising_data_files/SR01_Presenting_and_summarising_data_25_1.png


Measuring variability
~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    range(df1$hb)



.. raw:: html

    <ol class=list-inline>
    	<li>5.4</li>
    	<li>14.1</li>
    </ol>



.. code:: r

    quantile(df1$hb, c(0.25, 0.75))



.. raw:: html

    <dl class=dl-horizontal>
    	<dt>25%</dt>
    		<dd>8.75</dd>
    	<dt>75%</dt>
    		<dd>10.8</dd>
    </dl>



.. code:: r

    sd(df1$hb)



.. raw:: html

    1.97291837512189


.. code:: r

    var(df1$hb)



.. raw:: html

    3.89240691489362


.. code:: r

    df1 %>% summarize(min=min(hb), 
                     max=max(hb), 
                     iqr=quantile(hb, 0.75)- quantile(hb, 0.25), 
                     sd=sd(hb),
                     var=var(hb))



.. raw:: html

    <table>
    <thead><tr><th scope=col>min</th><th scope=col>max</th><th scope=col>iqr</th><th scope=col>sd</th><th scope=col>var</th></tr></thead>
    <tbody>
    	<tr><td>5.4     </td><td>14.1    </td><td>2.05    </td><td>1.972918</td><td>3.892407</td></tr>
    </tbody>
    </table>



Using a convenience function
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    summary(df1)



.. parsed-literal::

           hb        
     Min.   : 5.400  
     1st Qu.: 8.750  
     Median : 9.800  
     Mean   : 9.869  
     3rd Qu.:10.800  
     Max.   :14.100  


Common distributions and simple transformations
-----------------------------------------------

**Exercise**. Read in the ``urea.csv`` data file into a data frame
``df2`` and name the column ``urea``.

.. code:: r

    df2 <- read.csv("data/urea.csv", header=FALSE, col.names=c("urea"))

.. code:: r

    head(df2, n=3)



.. raw:: html

    <table>
    <thead><tr><th scope=col>urea</th></tr></thead>
    <tbody>
    	<tr><td>16.007049</td></tr>
    	<tr><td>13.647212</td></tr>
    	<tr><td> 6.653046</td></tr>
    </tbody>
    </table>



.. code:: r

    g <- ggplot(df2, aes(x=urea)) 
    g <- g + geom_histogram(binwidth=1, fill="grey", color="darkgrey")
    g





.. image:: SR01_Presenting_and_summarising_data_files/SR01_Presenting_and_summarising_data_38_1.png


Skewness
^^^^^^^^

We say the data has a positive or right skew. This name comes from the
fact that there is a statistical measure called skewness that is
positive for long right tails, and negative for long left tails.

.. code:: r

    install.packages("e1071", repos = "http://cran.r-project.org")
    library(e1071)


.. parsed-literal::

    
    The downloaded binary packages are in
    	/var/folders/3l/tbmzdkss71152d8t9n1f8nx40000gn/T//RtmpdeIh9P/downloaded_packages


.. code:: r

    skewness(df2$urea)



.. raw:: html

    1.80179083039024


Log transform of the data
^^^^^^^^^^^^^^^^^^^^^^^^^

.. code:: r

    df2['trans'] = log(df2['urea'])

.. code:: r

    g <- ggplot(df2, aes(x=trans)) 
    g <- g + geom_histogram(binwidth=0.3, fill="grey", color="darkgrey")
    g





.. image:: SR01_Presenting_and_summarising_data_files/SR01_Presenting_and_summarising_data_44_1.png


.. code:: r

    head(df2)



.. raw:: html

    <table>
    <thead><tr><th scope=col>urea</th><th scope=col>trans</th></tr></thead>
    <tbody>
    	<tr><td>16.007049</td><td>2.773029 </td></tr>
    	<tr><td>13.647212</td><td>2.613535 </td></tr>
    	<tr><td> 6.653046</td><td>1.895075 </td></tr>
    	<tr><td> 5.107674</td><td>1.630744 </td></tr>
    	<tr><td>19.325193</td><td>2.961410 </td></tr>
    	<tr><td>10.141074</td><td>2.316594 </td></tr>
    </tbody>
    </table>



Finding geometric mean
^^^^^^^^^^^^^^^^^^^^^^

.. code:: r

    gm  <- function(x) {
        return(exp(mean(log(x))))
    }

.. code:: r

    df2 %>% summarise(mean=mean(urea), geom.mean=gm(urea))



.. raw:: html

    <table>
    <thead><tr><th scope=col>mean</th><th scope=col>geom.mean</th></tr></thead>
    <tbody>
    	<tr><td>10.97242</td><td>8.504802</td></tr>
    </tbody>
    </table>



Exercise
~~~~~~~~

**1** Load the file "ph.csv" into a data frame called ``df``. Label the
column ``ph``.


**2**. Plot a histogram of the data. Calculate the skewness of the
``ph`` column.


**3** Left skewed data can sometimes be made more "normal" by an
exponential transformation. That is, if the original data is :math:`x`,
the transformed data is :math:`e^x`. Create another column named
``trasn`` with the transformed data.


**4**. Plot a histogram of the transformed data.


**5** Write your own function called ``sr.sd`` to calculate the standard
deviation using the formula in Table 3.


**6**. Create a new table using the ``summarise`` function with the mean
and standard deviation of both original and transformed data values.
This data frame should have 1 row and 4 columns named ``orig.mean``,
``orig.sd``, ``trans.mean``, ``trans.sd``.