---
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)
```

## Normal Likelihood and Regression

```{r GetData}
Data <- read.csv("../Data/EggBodySize.csv")

Data$IsKiwi <- grepl("Apteryx", Data$species)
Data$IsSwift <- grepl("Apus", Data$species)
Data$IsOwl <- grepl("Tyto", Data$species)

Data$logEggMass <- log(Data$Egg.mass..g.)
Data$logFemaleMass <- log(Data$female.mass..g.)
lnEggMass <- Data$logFemaleMass[!is.na(Data$logFemaleMass)]

```

We want to get to multiple regression, ANOVA etc. etc, but we will build up to it

Start with estimating the parmeters of a simple normal distribution


## Normal Likelihoods 


Some data: egg size in `r nrow(Data)` species of bird

```{r HistData}
par(mar=c(4.1,4.1,1,1))
hist(lnEggMass, xlab="Log egg mass", main="", lwd=3, col="lightgreen")

```

## The Normal Distribution

$$
f(x| \mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(x - \mu)^2}{2\sigma^2} }
$$


```{r NormalCurve, fig.height=5}
par(mar=c(4.1,4.1,1,1), lwd=3)
# curve(dnorm, from = -4, to = 4, xlab="X", ylab="Density")
NormPts <- seq(min(lnEggMass), max(lnEggMass), length=50)
FittedNorm <- dnorm(NormPts, mean(lnEggMass), sd(lnEggMass))

par(mar=c(4.1,4.1,1,1))
hist(lnEggMass, xlab="Log egg mass", main="", 
              freq=FALSE, lwd=3, col="lightgreen")
lines(NormPts, FittedNorm, col=2, lwd=2)
```

## Our Task

Estimate the parameters of this distribution
- estimate $\hat{\mu}$ and $\hat{\sigma}^2$

These parameters tell us about how big eggs are, and how much variation there is
- more complicated models have parameters, these parameters should say something useful scientifically


## Normal Probabilities

The normal is continuous, so $Pr(X=x) = 0$, and we talk about probabilities of being between 2 values:

$$
Pr(a < x < b | \mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \int_{a}^{b}{e^{\frac{1}{2\sigma^2} (x - \mu)^2}} dx 
$$
```{r NormalProb, fig.height=5}
a <- 0.5; b <- 1.7
xx <- seq(a, b, length=100)
yy <- dnorm(xx)
par(mar=c(4.1,4.1,1,1), lwd=3)
curve(dnorm, from = -4, to = 4, xlab="X", ylab="Density", main="", ylim=c(0, 1.04*dnorm(0)), yaxs="i",
      panel.first={polygon(c(xx[1], xx, xx[length(xx)]), c(0, yy, 0), border=NA, col="lightgreen")})
axis(1, at=c(a, b), labels=c("a", "b"), cex.axis=2, padj=1)
```

## Why the Normal Distribution is Nice

The shape of the normal does not change with the parameters. 

- Changing the mean moves the whole distribution
- Changing the variance stretches the distribution

Lots of things look normal

- Central Limit Theorem


## The Normal Likelihood

For one data point the likelihood is the probability density function ("pdf"):

$$
p(x_1|\mu, \sigma^2) =  \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(x_1 - \mu)^2}{2\sigma^2} }
$$


```{r NormalLhood, fig.height=2, fig.width=4}

par(mar=c(4.1,4.1,1,1))
plot(NormPts, FittedNorm, type="l", col=2, lwd=2, ann=FALSE, 
     panel.first={segments(lnEggMass[1], 0, 
                           lnEggMass[1], 
                           dnorm(lnEggMass[1], mean(lnEggMass), sd(lnEggMass)), col=4)})
```



## The Normal Likelihood: lots of data

With many data poits, we assume they are independent (given the parameters). The likelihood is simply the product of each of the likelihoods:
$$
p(x_1, x_2, x_3, ..., x_n|\mu, \sigma^2) = p(x_1) p(x_2) p(x_3)... p(x_n) = \prod_{i=1}^n {p(x_i|\mu, \sigma^2)}
$$

```{r NormalLhoods, fig.height=2, fig.width=4}
par(mar=c(4.1,4.1,1,1))
plot(NormPts, FittedNorm, type="l", col=2, lwd=2, ann=FALSE, 
     panel.first={segments(lnEggMass[1:10], rep(0, 10), 
                           lnEggMass[1:10], 
                           dnorm(lnEggMass[1:10], mean(lnEggMass), sd(lnEggMass)), col=4)})

```


## The Full Normal Likelihood

Some notation: just like we use $\sum$ to mean a summation, we use $\prod$ to mean a product (i.e. multiplying terms together)

$$
p(x_1, x_2, x_3, ..., x_n|\mu, \sigma^2) = \prod_{i=1}^n {p(x_i|\mu, \sigma^2)}
$$

Plugging in the normal pdf we get

$$
p(x_1, x_2, x_3, ..., x_n|\mu, \sigma^2) =  \frac{1}{(2 \pi \sigma^2)^\frac{n}{2}} \prod_{i=1}^n e^{-\frac{(x_i - \mu)^2}{2\sigma^2} }
$$

i.e. scary

## The Normal Likelihood

We can make this simpler by looking at the log-likelihood

$$
l(\mathbf{x}|\mu, \sigma^2) = \sum_{i=1}^n{l(x_i|\mu, \sigma^2)} = -\frac{n}{2} \log{2 \pi} - \frac{n}{2} \log{\sigma^2} - \sum_{i=1}^n {\frac{(x_i - \mu)^2}{2\sigma^2} }
$$

and removing constants

$$
l(\mathbf{x}|\mu, \sigma^2) = - \frac{n}{2} \log{\sigma^2} - \sum_{i=1}^n {\frac{(x_i - \mu)^2}{2\sigma^2} }
$$

So this is just quadratic in $x_i$

## Maximising the likelihood

We want estimates of $\mu$ and $\sigma^2$, so we need to find values of $\mu$ and $\sigma^2$ which give the maximum of the likelihood

(Quick Break??)

## Maximising the likelihood

There are two parameters: $\mu$ and $\sigma^2$

We can find the MLE of $\lbrace \mu, \sigma^2 \rbrace$ by looking at these separately, i.e. solve $\frac{\partial l(\mathbf{x}|\mu, \sigma^2)}{\partial \mu}$ and $\frac{\partial l(\mathbf{x}|\mu, \sigma^2)}{\partial \sigma^2}$

This is mathematically more involved. We do it by finding $\hat{\mu}$, the the LME for $\mu$:

$$
\begin{aligned}
l(\mathbf{x}|\mu, \sigma^2) & = - \frac{n}{2} \log{\sigma^2} - \frac{1}{2\sigma^2} \sum_{i=1}^n {(x_i^2 - 2x_i \mu + \mu^2)} \\
 & = - \frac{n}{2} \log{\sigma^2} - \frac{1}{2\sigma^2} \left( \sum_{i=1}^n x_i^2 -  2\mu \sum_{i=1}^n x_i + n \mu^2 \right)
\end{aligned}
$$

Differnetiate w.r.t (w.r.t = "with respect to") $\mu$, set to zero, re-arrange, and pray

## The Solutions

The MLE for the mean is 

$$
\hat{\mu} = \frac{\sum_{i=1}^n x_i}{n}
$$

For the variance, differentiate w.r.t $\sigma^2$, set to zero, re-arrange, pray some more, and get

$$
\hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \hat{\mu})^2
$$

For details, do it yourself or see https://www.statlect.com/fundamentals-of-statistics/normal-distribution-maximum-likelihood 


## Comments

The estimate $\hat{\mu}$ is just the sample mean, and $\hat{\sigma}^2$ is just the sample variance

- the whole distribution can be summarised by these two statistics

$\hat{\sigma}^2$ has $n$ as a denominator, not $n-1$

- because we assume the MLE for $\hat{\mu}$: using (n-1) is better because it takes into account the uncertainty


## What it looks like

```{R LHood, fig.height=4}
# mu <- 0; s2 <- 1; n <- 20
# dat <- rnorm(n, mu, sqrt(s2))

Pars <- expand.grid(
  mu = seq(min(lnEggMass), max(lnEggMass), length=1201), 
  s2 = seq(0.01,5, length=500)
)

# sum(!is.finite(Pars$lhood))
Pars$lhood <- apply(Pars, 1, function(ppars, data) {
  sum(dnorm(data, ppars["mu"], sqrt(ppars["s2"]), log=TRUE))
}, data=lnEggMass)

mu.hat <- mean(lnEggMass)
s2.hat <- (length(lnEggMass)-1)*var(lnEggMass)/length(lnEggMass)


MLE <- which(Pars$lhood==max(Pars$lhood))
Pars$ML.mu <- zapsmall(Pars$mu-Pars$mu[MLE])==0
Pars$ML.s2 <- zapsmall(Pars$s2-Pars$s2[MLE])==0


Pars$Explhood <- exp(Pars$lhood-max(Pars$lhood))
par(mfrow=c(1,2), oma=c(0.1,2,1,0), mar=c(4.1,2.1,1,1), ann=TRUE, lwd=3)
plot(Pars$mu[Pars$ML.s2], Pars$Explhood[Pars$ML.s2], type="l", 
     xlab=expression(mu), ylab="", main="Mean")
points(lnEggMass, runif(length(lnEggMass), 0, max(Pars$Explhood)/2), col="grey70", cex=0.5, pch=18)
abline(v=mu.hat, lty=3, col=2); 
plot(Pars$s2[Pars$ML.mu], Pars$Explhood[Pars$ML.mu], type="l", 
     xlab=expression(sigma^2), ylab="", main="Variance")
abline(v=s2.hat, lty=3, col=2); 
mtext("Likelihood", 2, outer=TRUE)

```

## Joint Likelihood

The likelihoods for $\hat{\mu}$ and $\hat{\sigma}^2$ are independent

```{R JointLHood, results='hide'}
library(lattice)
levelplot(Explhood ~ mu + s2, data=Pars, 
          xlim=c(5,7), ylim=c(2,5), 
          xlab=list(label=expression(mu), cex=2), 
          ylab=list(label=expression(sigma^2), cex=2), rot=0)
trellis.focus("panel", 1, 1, highlight=FALSE)
 lpoints(Pars$mu[Pars$lhood==max(Pars$lhood)], Pars$s2[Pars$lhood==max(Pars$lhood)], pch=4, cex=2, col=1)
trellis.unfocus()

```

## Uncertainty

Different realisations of the data give different results

```{R RepeatedSamples, results='hide'}
rep.y <- rep(1:10, each=10)
rep.x <- rnorm(length(rep.y))

MLEs <- data.frame(
  y = aggregate(rep.y, list(as.factor(rep.y)), mean)$x,
  hat.mu = aggregate(rep.x, list(as.factor(rep.y)), mean)$x,
  hat.s2 = aggregate(rep.x, list(as.factor(rep.y)), function(x) (length(x)-1)*var(x)/length(x))$x
)

par(cex=1.5)
plot(rep.x, rep.y, type="n", xlab="Simulated Data", ylab="Replicate")
abline(v=0)
apply(MLEs, 1, function(mle) {
  # segments(mle["hat.mu"]-sqrt(mle["hat.s2"]), mle["y"], mle["hat.mu"]+sqrt(mle["hat.s2"]), mle["y"], 
  #          lwd=5, col="lightgreen")
  segments(0, mle["y"], sqrt(mle["hat.s2"]), mle["y"], lend=3, 
           lwd=10, col="lightgreen")
  abline(h=mle["y"], lty=3)
  points(mle["hat.mu"], mle["y"], col="darkgreen", cex=2, pch=16)
})
points(rep.x, rep.y, pch=4, col="blue", cex=1)

```

We can summarise this with standard errors and confidence intervals

## Back to maths

The first differential of the log-likelihood w.r.t. $\mu$ is
$$
\frac{\delta l}{\delta \mu} = -\frac{1}{\sigma^2} \left(\sum_{i=1}^nx_i  - n \mu \right)
$$
If we differentiate again w.r.t $\hat{\mu}$ we get 

$$
\frac{\delta^2 l}{\delta \mu^2} = -\frac{n}{\sigma^2}
$$

The second differential is the curvature of the likelihood, and is the negative inverse of the standard error

This turns out to be a general assypmtotic approximation

- i.e. usually a good approximation, with enough data

## Curves in Plots

```{R Curvature}

dlhood <- function(x, mean, sd) {
  -(sum(x) - length(x)*mean)/sd
}

Datall <- scale(Data$logEggMass)
DataSome <- sample(Datall, size=100)
stderr.all <- sd(Datall)/sqrt(length(Datall))
stderr.some <- sd(DataSome)/sqrt(length(DataSome))
# ymax <- dnorm(0, sd=min(stderr.all, stderr.some, log=TRUE))
Range.y <- range(dnorm(Datall[abs(Datall)<1], mean=mean(Datall), sd=stderr.all, log=TRUE))
At.legend.y <- Range.y[1] + 0.6*(Range.y[2]-Range.y[1])
par(mfrow=c(2,1))
curve(dnorm(x=mean(Datall), mean=x, sd=stderr.all, log=TRUE), 
      from = -1, to = 1, xlab="X", ylab="Density", main="Likelihood") #, ylim=c(0, 1.04*ymax), yaxs="i")
curve(dnorm(x=mean(DataSome), mean=x, sd=stderr.some, log=TRUE), 
      from = -1, to = 1, col=2, add=TRUE)
legend(-0.4, At.legend.y, c("Full data", "100 points"), lty=1, col=1:2)

curve(dlhood(x=Datall, mean=x, sd=stderr.all), 
      from = -1, to = 1, xlab="X", ylab="Density", main="Slope of log-likelihood")
curve(dlhood(x=DataSome, mean=x, sd=stderr.some), 
      from = -1, to = 1, col=2, add=TRUE)
abline(h=0, lty=3)
```

## Confidence Intervals for $\hat{\mu}$

$\hat{\mu}$ has a Normal likelihood with standard error $\hat{\sigma}^2/\sqrt{n}$

- proof not shown!

So we can use a Normal distribution to create the (asymptotic) ML confidence intervals

95% Confidence interval = $(\hat{\mu} - 1.96 \sqrt{\hat{\sigma}^2/n}, \hat{\mu} + 1.96 \sqrt{\hat{\sigma}^2/n})$

## Confidence Intervals for $\hat{\sigma}^2$

The (asymptotic) standard error for $\hat{\sigma}^2$ is $\sqrt{2/n} \hat{\sigma}^2$, so the confidence interval is _approximately_

$$
(\hat{\sigma}^2 - 1.96 \sqrt{\frac{2}{n}} \hat{\sigma}^2, \hat{\sigma}^2 + 1.96\sqrt{\frac{2}{n}} \hat{\sigma}^2)
$$

## Why asymptotic?

The standard error is calculated at the MLEs of the other parameters

So any uncertainty in $\hat{\sigma}^2$ is ignored

This is also why the confidence interval is normally distributed, not a t-distribution

When $N$ is large, it is OK

For $\sigma^2$, the likelihood is not symmetric, but at reasonable sample sizes this does not matter

## Coverage 

```{R SimCoverage, cache=TRUE}
SimIsInCI <- function(N, mu=NULL, sd=NULL) {
  IsIn <- function(x, mn, stderr) {
    (mn - 1.96*stderr < x) & (mn + 1.96*stderr > x)
  }
  if(is.null(mu)) mu <- 0
  if(is.null(sd)) sd <- 1
  x <- rnorm(N, mean=mu, sd=sd)
  mu.hat <- mean(x)
  var.hat <- (N-1)*var(x)/N
  mu.cover <- IsIn(x = mu, mn = mu.hat, stderr = sqrt(var.hat/N))
  var.cover <- IsIn(x = sd, mn = var.hat, stderr = sqrt(2*var.hat/N))
  c(mu = mu.cover, sigma2 = var.cover)
}
CalcCoverage <- function(n, NR) {
  reps <- replicate(NR, SimIsInCI(n))
  c(n=n, apply(reps, 1, mean))
}

N <- seq(10, 1000, by = 50)
NRep <- 1e6
library(parallel)
Coverage <- t(mcmapply(CalcCoverage, n=N, NR=NRep, mc.cores = 4))
```

The coverage is the proportion of times the confidence interval contains the true value. For a 95% confidence interval, this should be 95%.

```{R Coverage, fig.height=5, cache=TRUE}
par(mfrow=c(1,2), mar=c(2,2,3,1), oma=c(2,2,0,0))
plot(Coverage[,"n"], Coverage[,"mu"], xlab="", ylab="", ylim=range(Coverage[,c("mu", "sigma2")]),
     main=expression(paste("Coverage for ", mu)))
abline(h=0.95, col=2)
plot(Coverage[,"n"], Coverage[,"sigma2"], xlab="", ylab="", ylim=range(Coverage[,c("mu", "sigma2")]), 
     main=expression(paste("Coverage for ", sigma^2)))
abline(h=0.95, col=2)
mtext("Sample Size", 1, outer=TRUE)
mtext("Coverage", 2, outer=TRUE)
```


## Does the model fit?

How could the model be wrong?

The assumptions:

- data follow a normal distribution
- data have a constant mean and variance
- data are independent


## Does the model fit?

This is a simple model, but it could still fit badly. For example, this data does not look like a normal distribution

```{R UniformData}
FittedNorm <- function(x, len=50) {
  x.pts <- seq(min(x), max(x), length=len)
  norm.pts <- dnorm(x.pts, mean(x), sd(x))
  cbind(x=x.pts, y=norm.pts)
}
PlotNorms <- function(x, len, ...) {
  Pts <- FittedNorm(x=x, len=50)
  hist(x, freq=FALSE, lwd=3, col="lightgreen", ...)
  lines(Pts[,"x"], Pts[,"y"], col=2, lwd=2)
}
N <- 500
x.unif <- runif(N)

par(mar=c(4.1,4.1,1,1))
PlotNorms(x.unif, xlab="Uniform Data", main="")

```

## Checking Normality

We should check if our model is good, i.e. are the assumptions OK?

If not, we might have to change the model


## Checking Normality

For our bird egg data, the main question is whether the data are normally distributed

Later we will ask if we can improve the model by trying to explain the variation

## Some Examples

Some examples of good/bad fit

```{R NormalExamples, fig.height=4}
x.norm <- rnorm(N)
x.t <- rt(N, df=2)
x.gamma <- rgamma(N, 2)
x.outlier <- c(x.norm[-1], 10)

NormPts <- FittedNorm(x=x.norm, len=50)
tPts <- FittedNorm(x=x.t, len=50)
GammaPts <- FittedNorm(x=x.gamma, len=50)
OutlierPts <- FittedNorm(x=x.outlier, len=50)

par(mfrow=c(1,5), mar=c(4.1,4.1,3,1))
PlotNorms(x.norm, xlab="Normal Data", main="Normal: Good!")
PlotNorms(x.unif, xlab="Uniform Data", main="Uniform")
PlotNorms(x.gamma, xlab="Gamma Data", main="Gamma (Chi-squared)")
PlotNorms(x.t, xlab="t Data", main="\"Student's\" t")
PlotNorms(x.outlier, xlab="t Data", main="Normal With Outlier")

```


## Quantile Plots

If we sample $n$ points from some distribution, we can work out what the expected values are. One way to do this is to split the distribution up into equal slices, and put the expected data where the slices are

- one point: the expected point is the mean (or median!)
- two points: 33% & 67% quantiles
- three points: 25%, 50% (median) & 75% quantiles

```{R Quantiles, fig.height=3}
par(mfrow=c(1,4))
curve(dnorm, from = -4, to = 4, xlab="X", ylab="Density", main="", ylim=c(0, 1.04*dnorm(0)), yaxs="i")
lines(c(0,0), c(0, dnorm(0)), lwd=4, col=4)

twoPts <- qnorm(c(1/3, 2/3))
curve(dnorm, from = -4, to = 4, xlab="X", ylab="Density", main="", ylim=c(0, 1.04*dnorm(0)), yaxs="i")
segments(twoPts, c(0, 0), twoPts, dnorm(twoPts, c(0, 0)), lwd=4, col=4)

threePts <- qnorm((1:3)/4)
curve(dnorm, from = -4, to = 4, xlab="X", ylab="Density", main="", ylim=c(0, 1.04*dnorm(0)), yaxs="i")
segments(threePts, c(0, 0), threePts, dnorm(threePts, c(0, 0)), lwd=4, col=4)

fourPts <- qnorm((1:4)/5)
curve(dnorm, from = -4, to = 4, xlab="X", ylab="Density", main="", ylim=c(0, 1.04*dnorm(0)), yaxs="i")
segments(fourPts, c(0, 0), fourPts, dnorm(fourPts, c(0, 0)), lwd=4, col=4)

```

## Normal Probability Plots

With lots of points, the quantiles should be about the same as the (ordered) points

So, we can line them up

If we plot them against each other they should lie on the 1:1 line


## Normal probability Plots

The normal distribution is nice because it has the same shape: we can move the x-axis or stretch it, but the shape is the same.

So we can plot the ordered data against the expected quantiles of a N(0,1) distribution with the same number of points

## Normal probability Plots

```{R QuantilePlots, fig.height=5}
par(mfrow=c(2,5), mar=c(4.1,4.1,3,1))
PlotNorms(x.norm, xlab="Normal Data", main="Good!")
PlotNorms(x.unif, xlab="Uniform Data", main="Thin Tails")
PlotNorms(x.gamma, xlab="Gamma Data", main="Skewed")
PlotNorms(x.gamma, xlab="t Data", main="Thick Tails")
PlotNorms(x.outlier, xlab="t Data", main="An Outlier")

qqnorm(x.norm); qqline(x.norm)
qqnorm(x.unif); qqline(x.unif)
qqnorm(x.gamma); qqline(x.gamma)
qqnorm(x.t); qqline(x.t)
qqnorm(x.outlier); qqline(x.outlier)

```

## Normal probability Plots

```{R EggQuantilePlots, fig.height=6}
par(mar=c(4.1,4.1,1,1), cex=2)
qqnorm(lnEggMass, col=4); qqline(lnEggMass)

```

Suggests tails are too thin for normal, but not massively so


## What we have Covered

The normal likelihood

- used a lot

Maximum likelihood with > 1 parameter

- maximumse w.r.t. each parameter

Standard errors and asymptotic confidence intervals

Model Checking

- does the model fit the data?


## Next Week

Regression

Allowing $\mu$ to change

More modelling, less inference