---
title: "ST2304 - Statistical Modelling for Biologists/Biotechnologists"
author: "Bob O'Hara"
date: "bob.ohara@ntnu.no"
output:
  beamer_presentation: default
  ioslides_presentation: default
  slidy_presentation: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

## Administration Matters

- Reference Group
- Blackboard



## Text Books

_Introductory statistics with R_ by Peter Dalgaard

![](Dalgaard.jpeg)

## Text Books

Grafen & Hails: Modern Statistics for the Life Sciences

![](GrafenHails.jpg)


## Previous Knowledge

What do you know?

(so I don't assume too much!)

## Mathematics

- logarithms/exponentials
- differeniation
- integration

## Probability

- Conditional Probability
- Expectation E(X)
- Indepedence
- Probability distribution


## Statistics

- Averages
- Variance, standard deviation
- Covariance, correlation
- t-tests
- regression (fitting a straight line)

## Programming

- Languages
    - C, BASIC, Python, Java, Python
    - R
- Object Oriented Programming


## Why do we need statistical modelling?

Because the world is complex!

## Face width and Fitness

Men with wider faces are more aggressive

But does this affect their fitness?

- survival in the Winter War
- number of (legitimate) offspring

## Face width and Fitness

But other things affect fitness
- regiment
- rank
- survival affects number of offspring

We need a model to disentangle all these 

## Face width and Fitness

Men with wider faces had more offspring, and might have slightly lower survival

![](FaceParams.png)

## Face width and Fitness

But thinner face -> higher rank

![](FaceRankFigure.png)


## Statistical Modelling: Brain size in birds

- What sort of question can we ask?

- What is likely to influence brain size?
    - actual brain size
    - what we measure


## Data Analysis

Different stages

Outline the full work flow

Modelling is only a part of this!

## First, collect the data

![](Johan.jpg)


## Format the data

* Get it into the right format so you can use it
* For big problems, use a data base
* Store the data in a convenient format, e.g. .csv

```{R ReadData, echo = TRUE}
BirdBrains <- read.csv("ece32961-sup-0001-TableS1.csv")
names(BirdBrains) <- gsub("Brain.*", "Brain", names(BirdBrains))
names(BirdBrains)[grep("resou", names(BirdBrains))] <- "Body"
```


## Look at the data

- Plot, calculate some summaries 
- Spot any possible problems
    - Typos etc.

```{R PlotData, fig.height = 4, echo = FALSE}
par(mfrow=c(1,1), mar=c(4.1,4.1,1,1))
plot(BirdBrains$Body, BirdBrains$Brain, type="p", xlab="Body Mass (g)", ylab="Brain Mass (g)")
```


## Look at the data


Nicer as a log-log plot?

```{R logPlotData, echo = FALSE}
par(mfrow=c(1,1), mar=c(4.1,4.1,1,1))
plot(BirdBrains$Body, BirdBrains$Brain, type="p", xlab="Body Mass (g)", ylab="Brain Mass (g)", log ="xy")
```


## Start the modelling

- What are the questions you are asking?
- What are problems in the data that you need to worry about?
- Do we need additional data?


## Start the modelling

How does body size affect brain size in birds?
    - is is similar to mammals?
    
Does the relationship differ between different groups of bird?


## Fit the model

Get the computer to do it
Look at the output

```{r FitMod, echo = TRUE}
bird.mod <- lm(log(Brain) ~ log(Body), data=BirdBrains)
summary(bird.mod)


# bird.mod2 <- lm(log(Brain) ~ log(Body) + Order, data=BirdBrains)
# order.mod <- lm(log(Brain) ~ log(Body), data=OrderBrains)
# 
# OrderBrains <- aggregate(BirdBrains[, c("Brain", "Body")], list(BirdBrains$Order), mean)
# 
# par(mfrow=c(1,1), mar=c(4.1,4.1,1,1))
# plot(BirdBrains$Body, BirdBrains$Brain, type="p", xlab="Body Mass (g)", ylab="Brain Mass (g)", log ="xy", 
#      col=1+(BirdBrains$Order=="Passeriformes"))
# points(OrderBrains$Body, OrderBrains$Brain, col=4, pch=16)
# 
# 
# par(mfrow=c(4,6), mar=c(1.1,1.1,3,1))
# sapply(levels(BirdBrains$Order), function(order, df) {
#   Use.df <- df[df$Order==order,]
#   plot(df$Body, df$Brain, type="n", log ="xy", main=order, cex=0.7)
#   points(Use.df$Body, Use.df$Brain)
# }, df=BirdBrains)
# 
# 
# Order.mods <- t(sapply(levels(BirdBrains$Order), function(order, df) {
#   Use.df <- df[df$Order==order,]
#   mod <- lm(log(Brain) ~ log(Body), data=Use.df)
#   c(coef(mod)["log(Body)"], confint(mod)["log(Body)",], N = nrow(Use.df))
# }, df=BirdBrains))
# Order.mods <- Order.mods[!is.na(Order.mods[,"log(Body)"]),]
# 
# par(mar=c(4.1,7.5,1,1))
# plot(Order.mods[,"log(Body)"], 1:nrow(Order.mods), xlim=c(0,1.5), lwd=2, yaxt="n", 
#      xlab="Coefficient", ylab="")
# segments(Order.mods[,"2.5 %"], 1:nrow(Order.mods), Order.mods[,"97.5 %"], 1:nrow(Order.mods), lwd=sqrt(Order.mods[,"N"])/2, lend=3)
# abline(v=coef(bird.mod)["log(Body)"], lty=3)
# axis(2, rownames(Order.mods), at=1:nrow(Order.mods), las=1)

```

## The Fitted Model

```{r PlotBrainsMod, echo = TRUE, fig.height=4, fig.width=7}
par(mfrow=c(1,1), mar=c(4.1,4.1,1,1))
plot(log(BirdBrains$Body), log(BirdBrains$Brain), type="p", xlab="log(Body Mass)", ylab="log(Brain Mass)")
abline(bird.mod, col=2)
text(2, 3, paste0("log(Body Mass) = ", round(coef(bird.mod)[1], 2), " + ", round(coef(bird.mod)[2], 2), " log(Brain Mass)"), adj=0)
```


## Check the Model

- Does it fit well?
- Outliers?
- Assumptions OK?
e.g. linearity
- Can you understand the results?

## Check the Model

```{r PlotMod, echo = TRUE, fig.height=4.5}
par(mar=c(4.1,4.1,0.1,1))
plot(fitted(bird.mod), resid(bird.mod), type="p", xlab="Fitted Value", ylab="Residual")

```

## Interpret the model

```{r MammalMod, echo = FALSE}
MammalBrains <- read.table("mammals.dat", header = TRUE)
mammal.mod <- lm(log(brain) ~ log(body), data=MammalBrains)

bird.mod2 <- lm(log(Brain) ~ log(Body) + Order, data=BirdBrains)
```
Back to biology!

The slope is `r round(coef(bird.mod)[2], 2)`, with a 95% confidence interval of `r round(confint(bird.mod, "log(Body)")[1], 2)` to `r round(confint(bird.mod, "log(Body)")[2], 2)`.  

This is less than it is for mammal, where the slope is `r round(coef(mammal.mod)[2], 2)` (95% confidence interval: `r round(confint(mammal.mod, "log(body)")[1], 2)` to `r round(confint(mammal.mod, "log(body)")[2], 2)`

Part of the reason for this is variation between bird orders. If we correct for order, the slope is `r round(coef(bird.mod2)[2], 2)`, 95% confidence interval: `r round(confint(bird.mod2, "log(Body)")[1], 2)` to `r round(confint(bird.mod2, "log(Body)")[2], 2)`. 

## A Quick Overview of R

A statistics programme

The language (S) designed for statistics

Object oriented

Interpreted language

## Basic Syntax: Assignment

```{r Assign, echo = TRUE}
x <- 2
x
(y <- 1:3)
(z <- c(4.5,6.7,-3))
```

## A Calculator

```{r Sums, echo = TRUE}

x + y

exp(z)

z/y # = z[1]/y[1] z[2]/y[2] z[3]/y[3] 
```


## Simple Plotting

```{r PlotYZ, echo = TRUE}
plot(y, z, col=2)

```

## Functions & Objects

exp(z) & plot(y, z) are functions

y & z is an _arguments_: they are given to the function, which does something with them

The function returns an object: 
- exp() returns a vector 
- plot() returns NULL (i.e. nothing: the plot is a side-effect)


## Statistical Models

lm() is a more complicated function

We give it 'log(Brain) ~ log(Body)' and BirdBrains as arguments

```{r MammalMod3, echo = TRUE}
bird.mod <- lm(log(Brain) ~ log(Body), data=BirdBrains)
```

It returns something horrid

## A big object

```{r MammalMod4, echo = TRUE}
str(bird.mod)
```

## The power of objects

We usually don't need to know what is in an object: there are functions to get us what we want

```{r BirdModObj, echo = TRUE}
coef(bird.mod)
summary(bird.mod)$coefficients
```

## Data Types

Our basic data can be of diffeent types, e.g. integer, real number, text, factor, logical 


```{r MammalMod5, echo = TRUE}
(txt <- c("thing", "stuff", "thing")) # text
(f <- factor(c("thing", "stuff", "thing"))) # factor
(lg <- c(TRUE, FALSE)) # logical
  
```
## Using data types: indices

```{r MammalMod6a, echo = TRUE}
(a1 <- -5:5)
a1[1]
a1[c(3,5,6)]

```

## Using data types: subsetting

```{r MammalMod6, echo = TRUE}
(logical.a1 <- a1 > 0)
(logical.a2 <- a1%%2==0)
a1[logical.a1]
a1[logical.a1 & logical.a2]

```


## Running an analysis

R has functions to read in data

```{r MammalMod7, echo = TRUE}
MammalBrains <- read.table("mammals.dat", header = TRUE)
```

and write it out

```{r MammalMod7a, echo = TRUE}
write.csv(MammalBrains, "mammals.csv")
save(MammalBrains, x, file="mammals.RData")
```


## Saving plots

We can save plots in files: pdf, png, jpg, etc etc

```{r MammalMod8, echo = TRUE, fig.height=2.3, fig.width=4}
par(mar=c(4.1,4.1,0,0))
plot(MammalBrains$body, MammalBrains$brain)

png("MammalBrains.png")
plot(MammalBrains$body, MammalBrains$brain)
dev.off()
```

## Functions

We can write our own functions....

```{r Function, echo = TRUE, fig.height=2.3, fig.width=4}
HelloWorld <- function(str) {
  if(!is.character((str))) 
    stop("str should be a character string you idiot")
  cat(paste0("Hello World, ", str, "!\n"))
}
HelloWorld("Orpheus")
```

## Packages

Packages put together a lot of functions (and objects etc.)

A lot of these are on CRAN

```{r Packages, echo = TRUE, fig.height=2.3, fig.width=4}
library(sp)
data(meuse)
coordinates(meuse) <- c("x", "y")
par(mar=c(2,2,1,1))
plot(meuse)

```


## Getting and Using R

You can download R from CRAN: http://www.r-project.org/

A lot of people use RStudio: http://rstudio.org/

## Finally

If you want to start using R now:

https://www.math.ntnu.no/emner/TMA4268/2018v/1Intro/Rbeginner.html