BIO708_Class2_Estimation.Rmd

---
title: "Objective functions and Estimation"
author: "Ian Dworkin"
date: "`r format(Sys.time(),'%d %b %Y')`"
output:
  slidy_presentation: 
    fig_retina: 1
    fig_width: 6
    incremental: no
    keep_md: yes
  beamer_presentation:
    incremental: no
  html_document: 
    keep_md: yes
    number_sections: yes
    toc: yes
  ioslides_presentation: 
    fig_height: 4
    fig_retina: 1
    fig_width: 6
    keep_md: yes
editor_options:
  chunk_output_type: console
---

# BIO708, Thinking about estimation (using Least Squares as an example)

##

<img src="curve_fitting_2x.png" alt="Estimation all depends" style="height: 700px; width:350px;"/>


https://xkcd.com/2048/


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(show.signif.stars=F)
options(digits=3)
```

## set up

```{r librariesToUse}
library(ggplot2)
```

## download data

```{r}
sct_data  <- read.csv("http://beaconcourse.pbworks.com/f/dll.csv",
                       h = T, stringsAsFactors = TRUE)
```
We are going to use this data set alot in this class. Below will download it from the web, but a local copy from here [download](https://www.dropbox.com/s/b9o0jqdhzc8suod/dll_data.csv?dl=0). 

## Data provenance and description

This data is from an [old paper](https://onlinelibrary.wiley.com/doi/10.1111/j.1525-142X.2005.05010.x), but from a kind of experiment we still often use, so it provides a useful template for this.

In this experiment I:

- Was testing a theory about the evolution of canalization, using a system to release cryptic genetic variation. 

- Introduced the *Dll[11]* mutation into about 25 different wild type strains of *Drosophila* via backcrossing. This mutation has some dominant phenotypes affecting leg (among other) traits, including length of segments and the number of sex comb teeth (in males).

- Reared flies from each of these strains and three different temperatures.

- Collected both *Dll[+]* and *Dll[11]* individuals from each of the strains from each of the temperature treatments.

- For each individual, a single first leg was dissected and three leg segments were measured along with the number of sex comb teeth.



## Let's start with a simple visualization using ggplot2

Today we are not going into checking and cleaning the data (we will next week). 

You can just remove any missing data for the moment

```{r}
sct_data <- na.omit(sct_data) 
```


## Simple scatterplor

```{r, echo=F}
ggplot(sct_data, aes(y = SCT, x = tarsus)) +
     geom_jitter(width = 0, height = 0.3, alpha = 0.25)
```

NOTE: I have added jitter (noise) to the SCT values, so we can get a better sense of the observations.

## Compare these two plots

```{r, echo=F}
ggplot(sct_data, aes(y = SCT, x = tarsus, color = genotype)) +
        geom_jitter(width = 0, height = 0.25, alpha = 0.2) +
        geom_smooth( method = loess)
```

## Compare these two plots

```{r, echo=F}
ggplot(sct_data, aes(y = SCT, x = tarsus, color = genotype)) +
        geom_jitter(width = 0, height = 0.25, alpha = 0.2) +
        geom_smooth( method = lm)
```

## Compare these two plots

- Same Data, but telling two different stories?

- Which is "best"?


## Compare these two plots

- Same Data, but telling two different stories?

- Which is "best"?

- Depends on what are goals are.

## What are we doing in each case?

```{r, echo=F}
ggplot(sct_data, aes(y = SCT, x = tarsus, color = genotype)) +
        geom_jitter(width = 0, height = 0.25, alpha = 0.2) +
        geom_smooth( method = lm)
```


## Objective functions and Estimation

>- We will discuss the various goals of estimation (not inference or prediction)
>- More than one way to skin a cat (different ways to estimate things)
>- The need for an objective function links the different ways of estimating things.

## What is the goal of estimation?

```{r, echo = F, include = F}
x <- seq(1, 10, by = 0.4)

y <- rnorm(length(x), 
           mean = 6 + 2*x - 1.5*(x^2) + 0.3*(x^3),  sd = 2*x)
```


```{r}
plot(y ~ x, pch = 20,
    ylab = "response variable", xlab = "predictor variable", 
    main = "What is the underlying process generating this data?")
```


## Fitting a model to the data. What do we mean by "best fit"

```{r}
plot(y ~ x, pch = 20,
    ylab = "response variable", xlab = "predictor variable", 
    main = "What is the underlying process generating this data?")

model_1 <- lm(y ~ poly(x,3))
param <- coef(model_1)
lines(x = x, y = fitted(model_1), col = "purple", lwd=2)
```

## Which do you think is the best fit? Why?

```{r}
plot(y ~ x, pch = 20,
    ylab = "response variable", xlab = "predictor variable", 
    main = "What is the underlying process generating this data?")

model_2 <- lm( y ~ x)
lines(x=x, y = fitted(model_2), col="purple", lwd = 2, lty = 2)
text(x=8,y=0, "Parametric models", col = "black")

model_3 <- lm( y ~ poly(x,2))
lines(x = x, y = fitted(model_3), col="blue", lwd = 2)

model_4 <- lm(y ~ poly(x, 3))
lines(x = x, y = fitted(model_4), col="red", lwd=2, lty=2)

model_22 <- lm(y ~ poly(x, 21))
lines(x = x, y = fitted(model_22), col="grey", lwd = 2)

legend("topleft", bty = "n",
       col = c("purple", "blue", "red", "grey"), lty = 1, 
       legend = c("linear", "quadratic", "cubic", "23rd order polynomial"))
```


## Which do you think is the best fit? Why?

```{r}
plot(y ~ x, pch = 20,
    ylab = "response variable", xlab = "predictor variable", 
    main = "a polynomial of order 3 (cubic)")

lines(x=x, y = fitted(model_2), col="purple", lwd = 2, lty = 2)

lines(x = x, y = fitted(model_3), col="blue", lwd = 2)

lines(x = x, y = fitted(model_4), col="red", lwd=2, lty=2)

lines(x = x, y = fitted(model_22), col="grey", lwd = 2)

text(x=4.25, y=(max(y) -5), cex = 1,
    expression(paste("~N( ", mu, " = ", 6 + 2*x -1.5*x^2 +0.3*x^3, ", ",sigma, " = 2x)")) )
```


## What happens when we use these estimates but for new data (same data generating process)?

```{r}
y2 <- rnorm(length(x), 
            mean = 6 + 2*x - 1.5*(x^2) + 0.3*(x^3),  
            sd= 2*x)

plot(y2 ~ x, pch = 20, 
     ylab = "response", xlab = "predictor",
     main="What happens to the model fit with new data (but same process)")

lines(x=x, y = fitted(model_2), col="purple", lwd=2, lty=2)
lines(x=x, y = fitted(model_3), col="blue", lwd =2, lty=1)
lines(x=x, y = fitted(model_4), col="red", lwd=2)
lines(x=x, y = fitted(model_22), col="grey", lwd=2)
```


## How do we go about determining the method to identify the "best" estimates.

- Let's leave aside how complex a polynomial we want to fit, and just think about a simple relationship that is seemingly linear.


```{r}
y3 <- rnorm(length(x), 
            mean = 6 + 2*x,  
            sd = 1.5)

plot(y3 ~ x, 
     ylab = "response", xlab = "predictor", pch = 20, cex = 2,
     main = "how do we get 'best' estimates of the slope and intercept?")

```


## In other words. Which line is the "line of best fit"?

```{r}
plot(y3 ~ x, 
     ylab = "response", xlab = "predictor", pch = 20, cex = 2,
     main = "how do we get 'best' estimates of the slope and intercept?")

abline(a = 16.7, b = 0, col = "grey", lwd = 1.2, lty = 3)
abline(a = 10, b = 1, col = "blue", lwd = 0.8, lty = 2)
abline(a = 3, b = 3, col = "red", lwd = 0.8, lty = 2)
abline(lm(y3 ~ x), col = "purple", lwd = 1, lty = 1)
```


## We need to decide on criteria that provide a basis to judge "best"
- Take a few minutes and try to think up at least two different criteria that could be used as a basis to assess "best fit"

## There is no one single best - deciding on your objective function.

>- As we will learn, there are many different definitions of "best"
>- We will need additional criterion to consider.
>- For now let us start with one you are most familiar with.

## Least squares estimation as a way of evaluating best fit

- Want to find our "best estimates" for parameters like the intercept ($\hat{\beta}_0$) and slope ($\hat{\beta}_1$)
- What is our objective function?
- We want to minimize a quantity, the sum of squared residuals.
- Usually this quantity is called the residual sum of squares.


$$ rss = \sum_{i=1}^{n}\epsilon^2_i$$
 Where $\epsilon_i$ is the residual (difference) between the *observed* value and the *estimated* value. $n$ is the number of observations you have.
 
## Least squares estimation

- In other words
$$
\epsilon_i = \hat{y}_i - y_i
$$

where

$$
\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i
$$

$\hat{\beta}_0$ is the estimated value for the intercept, $\hat{\beta}_1$ is the estimated value for the slope and $x_i$ is the observed predictor (independent) variable.


## The same thing one more way

$$ rss = \sum_{i=1}^{n}\epsilon^2_i$$
or
$$  = \sum_{i=1}^{n}(\hat{y}_i - y_i)^2$$

or 
$$ = \sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i)^2$$

Are all the same thing. 
- We try different values for $\beta_0$ and $\beta_1$ until we minimize this quantify, our objective function.


## Least squares estimation - why do we want to minimize?

What does it mean to minimize this quantity? 


```{r}
plot(y3 ~ x, 
     ylab = "response", xlab = "predictor", pch = 20, cex = 2,
     main = "how do we get 'best' estimates of the slope and intercept?")

abline(a = 16.7, b = 0, col = "grey", lwd = 1.2, lty = 3)
abline(a = 10, b = 1, col = "blue", lwd = 0.8, lty = 2)
abline(a = 3, b = 3, col = "red", lwd = 0.8, lty = 2)
abline(lm(y3 ~ x), col = "purple", lwd = 1, lty = 1)
```


## Least squares estimation - why do we want to minimize our function for RSS?
>- Minimizing reduces the difference between the "error" and the estimated values.
>- That is a sensible goal (think about the plot) for the RSS.
>- So our objective function is $\sum_{i=1}^{n}\epsilon^2_i$, which we want to minimize.
>- In other words we want to *minimize* the sum of squares residuals

## minimizing our objective function.
- So by this criterion the best estimates for our parameter are found by:

$$ \hat{\beta} = \underset{\beta}{\operatorname{argmin}} \sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i)^2$$

Which is ridiculous notation to say that we try different values of $\beta_0$ and $\beta_1$ to minimize the quantity we are interested in!


## What does this look like in R?

- Let's say we only needed to estimate the slope.

$$ \hat{\beta} = \underset{\beta}{\operatorname{argmin}} \sum_{i=1}^{n}(y_i - \beta_1 x_i)^2$$
- In R we would start by writing the objective function like this:

```{r, echo =T}
ResSumSq <- function(slope) { sum((y -slope*x )^2) } 
```

## Let's do this in R with simulated data

We are going to simulate some data to work with

First our covariate (predictor)
```{r, echo = TRUE}
x <- seq(from = 0, to = 10, by = 0.25)
```

- We are assuming for the time being that these are known, or measured WITHOUT error.

## Now we simulate our response

- We assume the relationship is

$$y \sim N(\mu = 0 + 1.25x, \sigma = 2)$$
  
```{r, echo = TRUE}
y <- rnorm(x, mean = 0 + 1.25*x, sd = 2)
```

- Importantly this means we assume that the intercept is 0
- We don't normally do this (even with simulated data)
- We are just doing this to make our first example a bit simpler

## Plot the relationship of the simulated data

```{r}
plot(y ~ x, pch = 20, col = "blue")
```

## Stochastic simulation

- We will use the `lm` function to fit the regression. We are specificying that the intercept will be zero.

```{r, echo = TRUE}
mod1 <- lm(y ~ 0 + x)
summary(mod1)
```

## estimates

```{r}
plot(y ~ x, pch = 20)
abline(a = 0, b = 1.25, col = "red")
abline(mod1, col = "blue")
```


## residual sum of squares from the model output

```{r, echo = T}
sum(resid(lm(y ~ 0 + x))^2)
```


## Our objective function for least squares

It is surprisingly simple in this case as we saw:

```{r, echo = T}
ResSumSq <- function(slope) { sum((y - slope*x )^2) } 
```


## So how do we find the "best" value

We could of course randomly try different values for the slope $\beta_1$.

- Try this, but for $\beta_1 = 0$


## how do we find the "best" value

- Try this, but for $\beta_1 = 0$

```{r, echo = T}
ResSumSq(slope = 0)
```

## We could try a few values now - Which are better?

- Try this for a few other values (bad or good guesses)

## We could try a few values now - which are better? 

```{r, echo = T}
ResSumSq(slope = 0.25)
```

```{r, echo = T}
ResSumSq(slope = - 1.25)
```

```{r, echo = T}
ResSumSq(slope = 1)
```

## This is really inefficient. So how would you do this more efficiently?
- Can you think about how you might go through lots of candidate values and try?

## how would you do this more efficiently?

- Let's create a variable called slope that we will use to store different numbers

```{r, echo=TRUE}
slope <- seq(-1.5, 3, by = 0.01)

head(slope)
```


## More efficiently trying numbers by "brute force"
- Nice people call this "grid approximation"
- For each candidate value in the vector "slope" calculate Residual sum of squares

```{r, echo = T}
rss_est <- sapply(slope, ResSumSq )

head(rss_est)
```

- Just put the trial values of the slope with rss values
```{r, echo = TRUE}
together <- cbind(slope, rss_est)
```

## candidate slope values and RSS


```{r}
plot(rss_est ~ slope, type = "p", pch = 20, cex = 0.5,
     ylim = c(0, max(rss_est)),
     ylab = "Sum of Squares", xlab = "Candidate slope")
```


## Where does the minimum RSS occur?

```{r, echo = TRUE}
brute_force_estimate <- together[which.min(together[,2]),] 

brute_force_estimate
```

- We can see that the estimate for $\beta_1$, and RSS is close (but not exactly the same) as we found for the linear model. Why?
- How could we improve these estimates using the brute force approach?

## Approximate, but not exactly the same


```{r}
plot(rss_est ~ slope, type = "p", pch = 20, cex = 0.5,
     ylim = c(0, max(rss_est)),
     ylab = "Sum of Squares", xlab = "Candidate slope")

segments(x0 = brute_force_estimate[1], x1 = brute_force_estimate[1],
         y0 = max(rss_est), y1 = brute_force_estimate[2], 
         lwd = 2, col = "purple", lty = 3)

segments(x0 = min(slope), x1 = brute_force_estimate[1],
         y0 = brute_force_estimate[2], y1 = brute_force_estimate[2], 
         lwd = 2, col = "blue", lty = 2)
```


## the computer has way better approaches

- We will use a default R optimizer to take our objective function, and find the minimum


```{r, echo = T}
optim_min <- optim(par = 3, ResSumSq , method="BFGS")

optim_min$par
optim_min$value
```


- From the call to lm the RSS

```{r, echo =T}
sum(resid(lm(y ~ 0 + x))^2)
```

- and the slope

```{r}
coef(mod1)
```


## Let's do a slightly more involved example with real data


```{r, echo =TRUE}
dll_data <- read.csv("http://beaconcourse.pbworks.com/w/file/fetch/35183279/dll.csv", header=TRUE)   #data frame input

dll_data <- na.omit(dll_data) # removing missing values
```


## Fit the model

```{r, echo =TRUE }
dll_lm_1 <- lm(SCT ~ 1 + tarsus, data=dll_data)
```

- The intercept and slope

```{r}
coef(dll_lm_1)
```

- RSS

```{r}
sum(resid(dll_lm_1)^2)
```


## Can you write the objective function for this?

- First, just extract the variables we want to use for the function

```{r, echo = T}
y <- dll_data$SCT
x <- dll_data$tarsus
```

## The objective function we want to minimize the RSS

```{r, echo = TRUE}
MinSumSq2 <- function(b) { 
	intercept <- b[1] # First value in b, corresponding to the intercept
	slope     <- b[2] # second value in b, corresponding to the slope of the model 
	sum((y -intercept -slope*x)^2)  } 
```

- We have just added the intercept (`b[1]`) and slope is now `b[2]`. This strange notation is courtesy of the `optim` function, and is a bit annoying at first. It wants all parameters estimated to be in a single vector (which we called `b`). 
- Once you start thinking in vectors, it is kind of nice.


## Now we can optimize (minimize)

```{r, echo = TRUE}
optim_min_2 <- optim(par=c(0, 0), 
                    MinSumSq2, method="BFGS")
```

- par provides the starting values to try for the optimizer

## What did we get?

```{r, echo=TRUE}
optim_min_2$par
optim_min_2$value
```


## We can double check this by putting these back into our objective function

```{r, echo = TRUE}
MinSumSq2(b = optim_min_2$par)
```

## some things to know.

- This is not actually how R computes these estimates. There are way more (computationally efficient ways of doing it).
- Beyond the scope of this workshop.


## What other objective functions could we consider?


## Why not this?

$$ = \sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i)$$


## How about this?

$$ = \sum_{i=1}^{n} |{ y_i - \beta_0 - \beta_1 x_i}|$$

This is actually called *least absolute deviation* and is used for *robust regression*.

## What is the objective function for the least absolute deviation 


```{r, echo = T}
LeastAbsDev <- function(b) { 
	intercept <- b[1] 
	slope     <- b[2] 
	sum(abs(y -intercept -slope*x))  } 
```

- and we run it


```{r, echo = T}
optim_min_3 <- optim(par=c(0.1, 0.1), 
                     LeastAbsDev, method="BFGS")

optim_min_2$par # LSE
optim_min_3$par # LAD
```

## penalized (regularized) estimation

- Ridge (L2 regularization of estimates).
- We minimize this quantity (the first part is the same RSS as above)
$$ \sum_{i=1}^{n}(y_i - \sum_{j=0}^{k}\beta_j x_{ij})^2 + \lambda \sum_{j=1}^{k}\beta_{j}^2$$
or, equivalently

$$ =  RSS + \lambda \sum_{j=1}^{k}\beta_{j}^2$$
- Lasso regression (L1 regularization of estimates
)
$$ \sum_{i=1}^{n}(y_i - \sum_{j=0}^{k}\beta_j x_{ij})^2 + \lambda \sum_{j=1}^{k}|\beta_{j}|$$
- $\lambda$ is a bit of a special parameter (the weight for the penalty on the coefficients). Usually determined with cross-validation.



## Let's try these out

```{r, echo = TRUE}

fakeRidge <- function(b, lambda = 0.2) { 
	intercept <- b[1] 
	slope     <- b[2] 
	(sum((y -intercept -slope*x)^2) + lambda * sum(b^2))}

fakeLasso <- function(b, lambda = 0.2) { 
	intercept <- b[1] 
	slope     <- b[2] 
	(sum((y -intercept -slope*x)^2) + lambda * sum(abs(b)))}
```

## How do these compare?

```{r, echo = TRUE}
optim_min_4 <- optim(par=c(0.1, 0.1), 
                     fakeRidge, method="BFGS")

optim_min_5 <- optim(par=c(0.1, 0.1), 
                     fakeLasso, method="BFGS")
optim_min_2$par #LSE
optim_min_4$par # Ridge (L2)
optim_min_5$par # lasso (L1)
```


## note

- For real lasso or ridge we do need to use cross-validation to determine optimal lambda.
- Think what would happen if we just allowed lambda to vary with the constraint of minimizing RSS?