What is a statistical model?

What is a statistical model?

Mathematical model

Statistical model

What is a statistical model?

• stochastic mathematical expression
• low-dimensional summary
• relates one or more dependent random variables to one or more independent variables

What is a statistical model?

• embodies a data generation process along with the distributional assumptions underlying this generation

• incorporates uncertainty

• \(response = model + error\)

What is the purpose of statistical modelling?

• describe relationships / effects
• estimate effects
• predict outcomes

How do we estimate model parameters? - \(y_i \sim{} \beta_0 + \beta_1 x_i\)

What criterion do we use to assess best fit?

• Depends on how we assume Y is distributed

If we assume \(y_i\) is drawn from a normal (gaussian) distribution…

• Ordinary Least Squares OLS

• parameters
  • location (mean)
  • spread (variance) - uncertainty

Least squares

Least squares

Least squares

Least squares

Least squares

Least squares

Least squares

Least squares

Least squares

Least squares

Least squares

Least squares estimates

• Minimize sum of the squared residuals

Least squares estimates

\[ y_i = \beta_0\times 1 + \beta_1\times x_i \]

\[ \underbrace{\begin{bmatrix} y_1\\ y_2\\ y_3\\ ...\\ y_i \end{bmatrix}}_{Y} = \underbrace{\begin{bmatrix} 1&x_1\\ 1&x_2\\ 1&x_3\\ ...\\ 1&x_i \end{bmatrix}}_{\boldsymbol{X}} \underbrace{\vphantom{\begin{bmatrix} y_1\\ y_2\\ y_3\\ ...\\ y_i \end{bmatrix}}\begin{bmatrix} \beta_0&\beta \end{bmatrix}}_{\boldsymbol{\beta}} \]

Least squares estimates

\[ \begin{align} Y =& \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\\ \boldsymbol{\varepsilon} =& Y - \boldsymbol{X}\boldsymbol{\beta} \end{align} \]

Minimize resid sum of squares (\(\boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon}\))

\[ \tiny \begin{bmatrix} \varepsilon_1&\varepsilon_2&\varepsilon_3&...&\varepsilon_n \end{bmatrix} \begin{bmatrix} \varepsilon_1\\ \varepsilon_2\\ \varepsilon_3\\ ...\\ \varepsilon_n\\ \end{bmatrix} = \begin{bmatrix} \varepsilon_1\times\varepsilon_1 + \varepsilon_2\times\varepsilon_2 + \varepsilon_3\times\varepsilon_3 + ... + \varepsilon_n\times\varepsilon_n \end{bmatrix} \]

Least squares estimates

\[ \begin{align} \boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon} =& (Y - \boldsymbol{X}\boldsymbol{\beta})^T(Y - \boldsymbol{X}\boldsymbol{\beta})\\ \boldsymbol{\beta} =& (\boldsymbol{X}^T\boldsymbol{X})^{-1} \boldsymbol{X}^T Y \end{align} \]

## Generate the model matrix
X = model.matrix(~1 + x, data = dat)
## Solve for beta
beta = solve(t(X) %*% X) %*% t(X) %*% dat$y
beta

##                 [,1]
## (Intercept) 2.650667
## x           4.574606

Least squares estimates

• Minimize sum of the squared residuals
• Simple matrix algebra

Provided data (and residuals):

• Gaussian
• Equally varied
• Independent

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

• Normality
• Homogeneity of variance
• Linearity
• Independence

What do we do, if the data do not satisfy the assumptions?

\[ \begin{align} y_i =& \beta_0 + \beta_1 x_i + \varepsilon_i\\ OR\\ y_i \sim{}& N(\mu_i, \sigma^2)\\ \mu_i =& \beta_0 + \beta_1 x_i \end{align} \]

• model embodies data generation processes
• pertains to:
  • effects (linear predictor)
  • distribution

• measurements typically away from zero

What about density (count / area)?

• no mass at 0

zero-bound variables with large var.

\(f(x\mid s, a) = \frac{1}{(s^a\Gamma(a))} x^(a-1) e^{-(x/s)}\)

\(a=shape, s=scale\)

\(\mu = as, \sigma^2 = as^2\)

What about abundance?

Count data

\(f(x\mid \lambda) = \frac{e^{-\lambda}\lambda^x}{x!}\)

\(\mu = \sigma^2 = \lambda = df\)

\(\theta~(dispersion) = \frac{\sigma^2}{\mu} = 1\)

Count data

\(f(x\mid \mu, \omega) = \frac{\Gamma(x + \omega)}{\Gamma(\omega)x!}\times\frac{\mu^x\omega^\omega}{(\mu + \omega)^{\mu+\omega}}\)

\(\theta~(dispersion) = 1/\omega\)

\(\omega = -\frac{\mu^2}{\mu - \sigma^2} \hspace{1cm} \theta = 0~(when~\omega=\infty)\)

What about binary and proportions?

Proportions or Presence/absence

\(f(x\mid n, p) = \binom{n}{p}p^x(1-p)^{n-x}\)

\(\mu = np, \sigma^2 = np(1-p)\)

for presence/absence \(n=1\)

What about percentages?

Continuous between 0 and 1

\(f(x\mid a, b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} x^{a-1} (1-x)^{b-1}\)

\(\mu = \frac{a}{a+b}, \sigma^2 = \frac{ab}{(a + b)^2.(a + b + 1)}\)

• must consider zero-one inflated

\[Y = \beta_0 + \beta_1x_1~+~...~+~\beta_px_p~+~e\] \[\underbrace{g(\mu)}_{Link~~function} = \underbrace{\beta_0 + \beta_1x_1~+~...~+~\beta_px_p}_{Systematic}\]

• Random component. \[Y \sim{} Dist(\mu, ...)\]
• Systematic component \([-\infty,\infty]\)
• Link function (\(g()\)) \[g(\mu) = \beta_0 + \beta_1x_1~+~...~+~\beta_px_p\]

Linear model is just a special case

\[Y = \beta_0 + \beta_1x_1~+~...~+~\beta_px_p~+~e\] \[\underbrace{g(\mu)}_{Link~~function} = \underbrace{\beta_0 + \beta_1x_1~+~...~+~\beta_px_p}_{Systematic}\]

• Random component. \[Y \sim{} N(\mu,\sigma^2)\]
• Systematic component \([-\infty,\infty]\)
• Link function (\(g()\)) \[I(\mu) = \beta_0 + \beta_1x_1~+~...~+~\beta_px_p\]

Gaussian: \[ P(x_i|\mu, \sigma) = \frac{1}{\sqrt{\sigma^2 2\pi}} exp^{-\frac{1}{2}((x_i-\mu)/\sigma)^2} \]

Likelihood:

\[ L(X|\theta) = \prod^n_{i=1} f'(x_i|\theta) \]

Product of the likelihood of each observation given parameters.

Log-likelihood:

\[ LL(X|\theta) = \sum^n_{i=1} log(f'(x_i|\theta)) \]

• Brute force
• Simplex based
• Derivative based
• Stochastic global
• Iterative re-weighted least squares

\[ \boldsymbol{\beta} = (\boldsymbol{X}^T\boldsymbol{W}\boldsymbol{X}^{-1})\boldsymbol{X}^T\boldsymbol{W}\boldsymbol{z} \]

glm(y~1, data=dat)

## 
## Call:  glm(formula = y ~ 1, data = dat)
## 
## Coefficients:
## (Intercept)  
##       27.81  
## 
## Degrees of Freedom: 9 Total (i.e. Null);  9 Residual
## Null Deviance:       2033 
## Residual Deviance: 2033  AIC: 85.53