---
title: "Lab 2 - Plot Models 4, 9, 10 and 11"
author: "WILD 502 - Jay Rotella"
output: pdf_document
---

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

To help promote understanding of the patterns in models 4, 9, 10, and 11, consider the following plots. First, we'll set up a data frame for making predictions, which has the full range of values of length for each sex. We then use functions in the `dplyr` package to add the product of sex * length. 

```{r, message=FALSE}
library(dplyr)
library(ggplot2)

dat <- expand.grid(int = 1, sex = c(0, 1), 
                   length = seq(from = 108, to = 136, by = 2))
dat <- dat %>% 
  mutate(sex.x.length = sex * length)
head(dat)
```

Next, we'll create the predicted values of log_odds of survival and of survival for each of the 4 models. To do so, we first create the design matrix for the model. The code below displays the first few rows of each design matrix so you can see the model structure clearly. Next, the code produces $X\beta$ (estimated log-odds of *S*) and $\frac{exp{(X\beta)}}{1 + exp{(X\beta)}}$ (estimated *S*). Then, we'll store all of the information on the estimates and which model they are from with the original data in a long data frame for use in plotting.

### S(length) or Model 4

```{r}
# S(length) or Model 4 or m4
m4.DM <- model.matrix(~ length, data = dat)
head(m4.DM, 4)

m4.beta <- matrix(c(-10.2340400, 0.0799874), nrow = 2, ncol=1)
m4.LO <- m4.DM %*% m4.beta
m4 <- data.frame(LO = m4.LO, S = plogis(m4.LO), mod = "S.len")
m4 <- cbind(dat[, 2:3], m4) # store sex, length, log_odds, S, and model name
```
\pagebreak

### S(sex + length) or Model 9

```{r}
# S(sex + length) or Model 9 or m9
m9.DM <- model.matrix(~ sex + length, data = dat)
head(m9.DM, 4)

m9.beta <- matrix(c(-12.6958420, -1.0554989, 0.1040179), nrow = 3, ncol=1)
m9.LO <- m9.DM %*% m9.beta
m9 <- data.frame(LO = m9.LO, S = plogis(m9.LO), mod = "S.sex.plus.len")
m9 <- cbind(dat[, 2:3], m9) # store sex, length, log_odds, S, and model name
```

### S(length + sex*length) or Model 10

```{r}
# S(length + sex*length) or Model 10 or m10
m10.DM <- model.matrix(~ length + length:sex, data = dat)
head(m10.DM, 4)

m10.beta <- matrix(c(-13.295076, 0.1090090, -0.0087417), nrow = 3, ncol=1)
m10.LO <- m10.DM %*% m10.beta
m10 <- data.frame(LO = m10.LO, S = plogis(m10.LO), mod = "S.len:sex")
m10 <- cbind(dat[, 2:3], m10) # store sex, length, log_odds, S, and model name
```

### S(Sex + Length + Sex*Length) or Model 11

```{r}
# S(sex + length + sex*length) or Model 11 or m11
m11.DM <- model.matrix(~ sex + length + length:sex, data = dat)
head(m11.DM, 4)

m11.beta <- matrix(c(-22.724123, 17.579626, 0.1858879, -0.1507537), nrow = 4, ncol=1)
m11.LO <- m11.DM %*% m11.beta
m11 <- data.frame(LO = m11.LO, S = plogis(m11.LO), mod = "S.sex.x.len")
m11 <- cbind(dat[, 2:3], m11) # store sex, length, log_odds, S, and model name
```
\pagebreak

### Plots

Sex-specific fitted lines for Models 9 (sex-specific intercepts, 1 slope) and Model 10 (1 intercept, sex-specific slopes) look alike. Sex-specific lines for Model 11 (sex-specific intercepts and slopes) do look different from those for Models 9 and 10, and $\Delta{AIC_c}$ provides a bit more support for that model. However, the uncertainty associated with each model's estimates makes it hard to select a single model among these 3.   

```{r, echo=FALSE, fig.width = 6, fig.height = 4}
dat.all <- rbind(m4, m9, m10, m11)
ggplot(dat.all, aes(x = length, y = LO, group = factor(sex))) + 
  geom_line(aes(color = factor(sex))) + facet_wrap(~ mod, nrow = 2)
```


```{r, echo=FALSE, fig.width = 6, fig.height = 4}
dat.all <- rbind(m4, m9, m10, m11)
ggplot(dat.all, aes(x = length, y = S, group = factor(sex))) + 
  geom_line(aes(color = factor(sex))) + facet_wrap(~ mod, nrow = 2)
```
\pagebreak

The plots below use a length range of 0 to 136, which covers values well below the minimum observed value of 108 cm (observed range highlighted in blue). This was done for educational purposes only to show the differences in intercepts and slopes more clearly on the plots of log-odds for Models 9 and 10; you would not want to predict survival for length values for which you have no observations.

```{r, echo = FALSE, , fig.width= 8, fig.height= 4}
# much bigger range of length values
dat <- expand.grid(int = 1, sex = c(0, 1), 
                   length = seq(from = 0, to = 136, by = 2))
dat <- dat %>% 
  mutate(sex.x.length = sex * length)

# S(length) or Model 4 or m4
m4.DM <- model.matrix(~ length, data = dat)
m4.beta <- matrix(c(-10.2340400, 0.0799874), nrow = 2, ncol=1)
m4.LO <- m4.DM %*% m4.beta
m4 <- data.frame(LO = m4.LO, S = plogis(m4.LO), mod = "S.len")
m4 <- cbind(dat[, 2:3], m4) # store sex, length, log_odds, S, and model name

# S(sex + length) or Model 9 or m9
m9.DM <- model.matrix(~ sex + length, data = dat)
m9.beta <- matrix(c(-12.6958420, -1.0554989, 0.1040179), nrow = 3, ncol=1)
m9.LO <- m9.DM %*% m9.beta
m9 <- data.frame(LO = m9.LO, S = plogis(m9.LO), mod = "S.sex.plus.len")
m9 <- cbind(dat[, 2:3], m9) # store sex, length, log_odds, S, and model name

# S(length + sex*length) or Model 10 or m10
m10.DM <- model.matrix(~ length + length:sex, data = dat)
m10.beta <- matrix(c(-13.295076, 0.1090090, -0.0087417), nrow = 3, ncol=1)
m10.LO <- m10.DM %*% m10.beta
m10 <- data.frame(LO = m10.LO, S = plogis(m10.LO), mod = "S.len:sex")
m10 <- cbind(dat[, 2:3], m10) # store sex, length, log_odds, S, and model name

# S(sex + length + sex*length) or Model 11 or m11
m11.DM <- model.matrix(~ sex + length + length:sex, data = dat)
m11.beta <- matrix(c(-22.724123, 17.579626, 0.1858879, -0.1507537), nrow = 4, ncol=1)
m11.LO <- m11.DM %*% m11.beta
m11 <- data.frame(LO = m11.LO, S = plogis(m11.LO), mod = "S.sex.x.len")
m11 <- cbind(dat[, 2:3], m11) 

# store sex, length, log_odds, S, and model name
dat.all <- rbind(m4, m9, m10, m11)

# make plots
ggplot(dat.all, aes(x = length, y = LO, group = factor(sex))) + 
  geom_line(aes(color = factor(sex))) + facet_wrap(~ mod, nrow = 2) +
  annotate("rect", xmin = 108, xmax = 136, ymin = -Inf, ymax = Inf, 
           alpha = 0.2, fill = "cornflowerblue") 
```


```{r, echo = FALSE, fig.width= 8, fig.height= 4}
ggplot(dat.all, aes(x = length, y = S, group = factor(sex))) + 
  geom_line(aes(color = factor(sex))) + facet_wrap(~ mod, nrow = 2) +
  annotate("rect", xmin = 108, xmax = 136, ymin = -Inf, ymax = Inf, 
           alpha = 0.2, fill = "cornflowerblue") 
```

\