---
title: "Statistical Inference"
author: "Bob O'Hara"
date: "bob.ohara@ntnu.no"
output: beamer_presentation
---

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

## Administration Matters

- Reference Group
- Blackboard



## Likelihood (again)



## Single normal distribution


## 2 means


## Regression


## What it is


# Fitting by ML


## What the curve looks like

## Assumptions of the model: linearity




## The general idea

Our data is $Y$. We sample it randomly from a population, according to some probability $P(Y | X, \theta)$

- $X$ are (fixed) covariates
- $\theta$ are parameters that describe the population

$P(Y | X, \theta)$ is a model that describes how we sample

We want an estimator of $\theta$

## An Example

What proportion of the earth is land?

If we have a globe, how can we estimate what proportion in land and what proportion sea?

(plant cover is a real example of this problem)

## Sampling The Earth

Toss the globe around
When you catch it. put your finger on a point, and say whether it lands on the land or sea
Then toss it to someone else

We will record the number of times we get Land or Sea, and use this as an estimate of the proportion of the glode that is land

## The Model

Each observation is a sample from the real world

- "Bernoulli trial"

We observe $N$ trials, of which $n$ are land, and $(N-n)$ are water

$n$ follows a binomial distribution, with an unknown $p$ (the population-level mean)

## The Model

$$
Pr(n = r | N, p) = \frac{N !}{r! (N-r)!} p^r (1-p)^{N-r}
$$

We want an estimator of $p$, which is the "true" population-level mean

## Terminology


We call $p$ the _estimand_: this is what we want an estimator of

We will call the _estimator_ $\hat{p}$

Because we will get $\hat{p}$ by maximising the likelihood, we call it the _maximum likelihood estimator_ (MLE).


## Likelihood

$Pr(n = r | N, p)$ as a function of $p$: call it the _likelihood_: $L(p|N, r)$

```{R BinomLhood, fig.height= 5}
N <- 10
p <- 0.4
r <- 4

pp <- seq(0,1,length=500)
lhood <- dbinom(x = r, size = N, prob = pp)

par(mar=c(4.1,4.1,1,1), lwd=2, cex=1.5, las=1)
plot(pp, lhood, type="l", xlab = "Probability", ylab="Likelihood")
```

Find the value of $p$ that maximises $Pr(n = r | N, p)$

## The Philosophy

The likelihood is a data generating mechanism: it is a statistical model

We assume that the data are random, and the parameters (and model) are fixed

We want to find the parameters which are most likely to give rise to the data

 - we maximise the likelihood

## Maximising the likelihood: two tricks

The maximum of the likelihood is at the same value of $p$ as the maximum of the log-likelihood

i.e. we work with $\log(L(p|N, r)) = l(p|N, r)$

$$
l(p | n) = \log(N !) - \log(r!) - \log((N-r)!) + r \log(p) + (N-r) \log(1-p)
$$

Becuse this is a function of $p$, not $N$ or $r$, several terms are constants:

$$
l(p | n) = r \log(p) + (N-r) \log(1-p) + C
$$

## The log-likelihood

```{R BinomlogLhood, fig.height= 5}

par(mfrow=c(1,2), mar=c(3,1.1,3,0.1), oma=c(2,2,0,0), lwd=2, cex=1.5, las=1)
plot(pp, lhood, type="l", main="Likelihood", xlab="", ylab="", yaxt="n")
plot(pp, log(lhood), type="l", main="log Likelihood", xlab="", ylab="", yaxt="n")
mtext("Probability", 1, outer = TRUE, las=0, cex=2, line=0)
mtext("Density", 2, outer = TRUE, las=0, line=0, cex=2)

```

## Maximising the log likelihood

Differentiate with respect to p

$$
l(p | n) = r \log(p) + (N-r) \log(1-p) + C
$$
using the chain rule for the second term ($u=1-p$):

$$
\frac{d \log(1-p)}{dp} = \frac{d \log u}{du}\frac{d u}{dp} = \frac{1}{u} (-1) = -\frac{1}{1-p}
$$
we get

$$
\frac{d l(p | n)}{dp} = \frac{r}{p} - \frac{N-r}{1-p}
$$

## In Figures


```{R BinomLhoodFigs, fig.height= 6}
Grad.dbinom=function(x, size, prob) {
  x/prob - (size - x)/(1-prob)
}
par(mfrow=c(1,3), mar=c(3,1.1,3,0.1), oma=c(2,2,0,0), lwd=2, cex=1.5, las=1)
plot(pp, lhood, type="l", main="Likelihood", xlab="", ylab="", yaxt="n")
plot(pp, log(lhood), type="l", main="log Likelihood", xlab="", ylab="", yaxt="n")
plot(pp, Grad.dbinom(x = r, size = N, prob = pp), type="l", main="Gradient of\nLikelihood", 
     xlab="", ylab="", yaxt="n")
abline(h=0, lty=3, col=2)
mtext("Probability", 1, outer = TRUE, las=0, cex=2, line=0)
mtext("Density", 2, outer = TRUE, las=0, line=0, cex=2)

```

## Maximising

Set the gradient to 0:

$$
0 = \frac{r}{p} - \frac{N-r}{1-p}
$$
So 

$$
\frac{p}{1-p} = \frac{r}{N-r}
$$

i.e. the _odds_ of success are equal to the ratio of successes to failure. We can re-write this as $\frac{1-p}{p} = \frac{N-r}{r}$, and re-arrange to get
xtable
$$
\hat{p} = \frac{r}{N}
$$

## So...

We have (analytically) maximised the likelihood to get an estimator of $p$

$$
\hat{p} = \frac{r}{N}
$$

In more complicated problems we do the same thing, but sometimes the maximisation is done numerically (or even through simulation)

But we always use the log-likelihood & ignore the normalising constants

## What happens if we take another sample?

??


## What happens if we take another sample?

Each sample gives us a different $\hat{p}$

$p$ is fixed, and the data are random, so $\hat{p}$ is a property of the data

We can sample repeatedly many times, and each time get a different $\hat{p}$

The likelihood is the distribution of $\hat{p}$

## Samples

```{R RepeatedPsNumbers}
N.rpt <- 5e2
p.rpt <- 0.4
r.rpt <- rbinom(1e3, N.rpt, p.rpt)/N.rpt
r.tab <- table(r.rpt)
nn <- 0:N.rpt
lhood.rpt <- dbinom(x = nn, size = N.rpt, prob = 0.4)
```

_N_ = `r N.rpt`, _p_ = `r round(p.rpt, 1)`

```{R RepeatedPsPlot}

plot(r.tab, lwd=15, lend=4, col="lightblue")
# hist(r.reps)
abline(v=p, lwd=3, col=2)
lines(nn/max(nn), lhood.rpt*max(r.tab)/max(lhood.rpt), lwd=2, lty=2)

```

## Uncertainty

Because different samples give different estimates, we want to quantify this
- suggest plausible values

What summaries could we use?

## Standard Errors 

Standard Deviations of Statistics: $s$

Binomial variance of $n$: $Var(n|N, p) =Np(1-p)$

Our statistic: $n/N$

$Var(n/N) = 1/N^2 Var(n) = p(1-p)/N$

Standard error:

$s = \sqrt{p(1-p)/N}$

## Confidence Intervals

Standard errors are useful, but what values of $p$ are likely?

Confidence intervals give the range of values in which the statistic is likely

- usually we use 95%

A 95% confidence interval is one that has a 95% probability of containing the estimate of $p$ if the estimate $\hat{p}$ is the true value.

## Constructing a Confidence interval

```{R CI}
x.norm <- seq(-3,3,length=50)
y.norm <- dnorm(x.norm)
par(mar=c(1,1,1,1), lwd=2)
plot(x.norm, y.norm, type="l")
LowerCI <- x.norm < -1.96
UpperCI <- x.norm > 1.96
polygon(c(x.norm[LowerCI], rev(x.norm[LowerCI])),
        c(y.norm[LowerCI], rep(0,sum(LowerCI))),
        col = "red")
polygon(c(x.norm[UpperCI], rev(x.norm[UpperCI])),
        c(y.norm[UpperCI], rep(0,sum(UpperCI))),
        col = "red")
text(0, max(y.norm)/2, "95%", cex=5)
text(mean(x.norm[UpperCI]), max(y.norm)/4, "2.5%", cex=3, col=2)
text(mean(x.norm[LowerCI]), max(y.norm)/4, "2.5%", cex=3, col=2)

```


## Confidence Interval for a Binomial

More difficult, because data are discrete
```{R RepeatedPsNumbers2}
N.small <- 5
p.small <- 0.4
r.small <- rbinom(1e3, N.small, p.small)/N.small
r.tab <- table(r.small)
nn <- 0:N.small
lhood.small <- dbinom(x = nn, size = N.small, prob = 0.4)
```

e.g. N = 5


_N_ = `r N.small`, _p_ = `r round(p.small, 1)`

```{R RepeatedPsPlot2}

Lines.off <- max(r.tab)/max(lhood.small)
plot(r.tab, lwd=15, lend=4, col="lightblue")
# hist(r.reps)
abline(v=p, lwd=3, col=2)
lines(nn/max(nn), lhood.small*Lines.off, lwd=2, lty=2)
text(0, lhood.small[1]*Lines.off+50, 
     round(lhood.small[1]/sum(lhood.small), 2))
text(1, rev(lhood.small)[1]*Lines.off+50, 
     round(rev(lhood.small)[1]/sum(lhood.small), 2))

```
With larger $N$, the problem is smaller

## Approximations

In general, when we have a lot of data the likelihood looks like a normal distribution

So we can use use that to make an approximation

## Approximating the Likelihood

The likelihood is approximately Normal with mean equal to the MLE, and standard deviation equal to the standard error
$$
l(p|N, n) \sim N(\hat{p}, \sqrt{\hat{p}(1-\hat{p})/N})
$$

```{R BinomApproxlogLhood, fig.height= 5}
N100 <- 100; p100 <- 0.4; r100 <- 40

phat <- r/N
stderr <- sqrt(phat*(1-phat)/N)
lhood.approx <- dnorm(pp, phat, stderr)
lhood.approx <- lhood.approx/max(lhood.approx)

lhood100 <- dbinom(x = r100, size = 100, prob = pp)
stderr100 <- sqrt(phat*(1-phat)/100)
lhood.approx100 <- dnorm(pp, phat, stderr100)
lhood.approx100 <- lhood.approx100/max(lhood.approx100)

par(mfrow=c(1,2), mar=c(3,1.1,3,0.1), oma=c(2,2,0,0), lwd=2, cex=1.5, las=1)
plot(pp, lhood, type="l", main="N = 10", xlab="", ylab="", yaxt="n")
lines(pp, lhood.approx*max(lhood), col=2)

plot(pp, lhood100, type="l", main="N = 100", xlab="", ylab="", yaxt="n")
lines(pp, lhood.approx100*max(lhood100), col=2)
legend(0.53, max(lhood100)*0.8, c("Exact", "Approximation"), col=1:2, lty=1, cex=0.6)
mtext("Probability", 1, outer = TRUE, las=0, cex=2, line=0)
mtext("Density", 2, outer = TRUE, las=0, line=0, cex=2)

```


## Approximating the Confidence Interval

We can use athis approximation, so the interval is

$$
\left( \hat{p} - 1.96 \sqrt{\hat{p}(1-\hat{p})/N},   \hat{p} + 1.96 \sqrt{\hat{p}(1-\hat{p})/N} \right)
$$

So, for $N=10$, $n=4$ the interval is (0.4 - `r round(sqrt(p*(1-p)/N),2)`, 0.4 + `r round(sqrt(p*(1-p)/N),2)`), i.e. (`r round(0.4-sqrt(p*(1-p)/N),2)`, `r round(0.4+sqrt(p*(1-p)/N),2)`).

With $N=100$, $n=40$ the interval is smaller: (`r round(0.4-sqrt(p100*(1-p100)/N100),2)`, `r round(0.4+sqrt(p100*(1-p100)/N100),2)`).

## Approximating the Confidence Interval

```{R PlotBinomApproxlogLhood, fig.height= 5}
CI10 <- 0.4 + 1.96*c(-stderr, stderr)
CI100 <- 0.4 + 1.96*c(-stderr100, stderr100)
Dens.x10 <- seq(CI10[1], CI10[2], length=500)
Dens.l10 <- dnorm(4/10, Dens.x10, stderr)
Dens.l10 <- Dens.l10/max(Dens.l10)

Dens.x100 <- seq(CI100[1], CI100[2], length=500)
Dens.l100 <- dnorm(4/10, Dens.x100, stderr100)
Dens.l100 <- Dens.l100/max(Dens.l100)

par(mfrow=c(1,2), mar=c(3,1.1,3,0.1), oma=c(2,2,0,0), lwd=2, cex=1.5, las=1)
plot(pp, lhood, type="l", main="N = 10", xlab="", ylab="", yaxt="n")
polygon(c(Dens.x10[1], Dens.x10, Dens.x10[length(Dens.x10)]), c(0, Dens.l10, 0)*max(lhood), col = "lightblue", border=NA)
lines(pp, lhood.approx*max(lhood), col=2)
lines(pp, lhood, col=1)

plot(pp, lhood100, type="l", main="N = 100", xlab="", ylab="", yaxt="n")
polygon(c(Dens.x100[1], Dens.x100, Dens.x100[length(Dens.x100)]), c(0, Dens.l100, 0)*max(lhood100), col = "lightblue", border=NA)
lines(pp, lhood.approx100*max(lhood100), col=2)
lines(pp, lhood100, col=1)

legend(0.53, max(lhood100)*0.8, c("Exact", "Approximation"), col=1:2, lty=1, cex=0.6)
mtext("Probability", 1, outer = TRUE, las=0, cex=2, line=0)
mtext("Density", 2, outer = TRUE, las=0, line=0, cex=2)

```

## How well are we doing?

The confidence interval is an interval where, if we repeat the same experiment many thimes, we have a 95% probability of getting an estimate inside the interval.

We can check this!

## How well are we doing?

```{R CheckCoverage, fig.height=5}
SimulatedData10 <- rbinom(n=1e4, size = 10, prob = 0.4)
SimulatedData100 <- rbinom(n=1e4, size = 100, prob = 0.4)

par(mfrow=c(1,2), mar=c(3,1.1,3,0.1), oma=c(2,2,0,0), lwd=2, cex=1.5, las=1, yaxt="n")
plot(table(SimulatedData10/10), lwd=10, lend=2, main="N = 10", 
     panel.first = {rect(Dens.x10[1], -100, Dens.x10[length(Dens.x10)], 10e4, col="lightblue", border=NA)})

plot(table(SimulatedData100/100), lwd=10, lend=2, main="N = 100", 
     panel.first = {rect(Dens.x100[1], -100, Dens.x100[length(Dens.x100)], 10e4, col="lightblue", border=NA)})

```

For N = 10 we have a proportion `r mean(SimulatedData10/10>Dens.x10[1] & SimulatedData10/10<Dens.x10[length(Dens.x100)])` within the CI, 

For N = 100 the proportion is `r mean(SimulatedData100/100>Dens.x100[1] & SimulatedData100/100<Dens.x100[length(Dens.x100)])` within the CI, 


## Summary

We want to estimate parameters 
- summaries statistics

Estiamtors are functions of the data: the estimand is fixed

Likelihood gives the probability of the data (and thus estimators)

We can summarise the uncertainty with standard errors and confidence intervals

Our example was the binomial distribution

## Next Week

Regression: fitting models with straight lines

- multiple parameters
- normal likelihood
- model fit




