---
title: "Lecture 3 - R code"
author: "WILD 502- Jay Rotella"
output: pdf_document
---

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

## Likelihood Values and Sample Size

First, let's take a look at how the likehood values for various values of *p* vary as we change the sample size but keep the proportion of successes the same. 

```{r, message=FALSE}
library(tidyverse) # for dplyr and ggplot2
library(binom) # for binomial confidence intervals
library(knitr) # for 'kable' function for printing nice table
true.p = 0.6
p=seq(from = 0, to = 1, by = 0.01) 

# make a dataframe with combinations of p and N
out <- expand.grid(p = seq(0, 1, 0.01), N = seq(10, 100, 30)) 
# add y and likelihood values for each value of p
out <- out %>% 
  mutate(y = true.p * N, 
         Lp = choose(N, y) * p^(y) * (1-p)^(N-y),
         lnLik = log(Lp))

# plot likelihoods for each value of N
ggplot(out, aes(x = p, y = Lp, color = factor(N))) +
          geom_line() + geom_vline(xintercept = true.p)
```

If we calculate confidence intervals based on the likelihood profile, we can see how the uncertainty in $\hat{p}$ changes as sample size changes.

```{r}
CIs = binom.profile(x = 0.6 * seq(10, 100, 30),
              n = seq(10, 100, 30), 
              conf.level = 0.95)
kable(CIs)
```

As explained in [Chapter 1 of CW](http://www.phidot.org/software/mark/docs/book/pdf/chap1.pdf)(*see page 1-22*), "... for a single parameter, likelihood theory shows that the 2 points 1.92 units down from the maximum of the log likelihood function provide a 95% confidence interval when there is no extra-binomial variation ...". The graphic below shows this for the scenario where N = 40 and y = 24. 

```{r}
out40 <- filter(out, N == 40, between(p, 0.3, 0.8))
max_ln_lik <- max(out40$lnLik)

ggplot(out40, aes(x = p, y = lnLik, color = factor(N))) +
  geom_line() + geom_vline(xintercept = true.p) +
  geom_hline(yintercept = max_ln_lik - 1.92, linetype = "dashed") +
  geom_vline(xintercept = c(CIs[2, 5], CIs[2, 6]), linetype = "dashed")
```


\pagebreak

## Converting between Probability and Log-odds

This simple exercise is intended to give you a better understanding of the connections between the log-odds of an outcome and the  probability of an outcome. Recall that if the probability of an event is 0.25, that the value for the log-odds is calculated as $ln(\frac{0.25}{0.75})$. In R, you can use the `qlogis` function to obtain the log-odds for a given probability, e.g., `qlogis(0.25)` = `r qlogis(0.25)`.

If you know the log-odds, then you can calculate the probability by using $\frac{exp(x)}{1 + exp(x)}$, where *x* represents the log-odds value. The In R, you can use the `plogis` function to obtain the probability for a given log-odds value, e.g., `plogis(-1.098612)` = `r plogis(-1.0986123)`.

Probability values range from 0 to 1. It turns out that for $\frac{exp(x)}{1 + exp(x)}$, values of *x* ranging from -5 to +5 create probabilities that range from just above 0 to very close to 1. Values of *x* ranging from -1 to +1 create probabilities that range from about 0.25 to 0.75. The material below will let you explore the relationships for yourself.

```{r}
log_odds = seq(from = -5, to = 5, by = 0.25)
# use 'plogis' function to calculate exp(x)/(1 + exp(x))
y = plogis(log_odds)
# store log_odds and y in data frame for use with ggplot
d = data.frame(log_odds, y)
head(d, 4)
tail(d, 4)
```

Below, we plot the relationship, so you can see the pattern among the values for log-odds and associated probabilities. You might wonder what happens if you get log-odds values that are very very small (e.g., -24, -147, or -2421) or very big (e.g.,14, 250, or 1250). You should use the `plogis` function on such values (no commas in your numbers, e.g., `plogis(-2421)`) to find out for yourself. 

```{r}
ggplot(d, aes(x = log_odds, y = y)) +
  geom_line() +
  geom_hline(aes(yintercept = 0.5),
             colour = "gray",
             linetype = "dashed") +
             geom_vline(aes(xintercept = 0.0),
             colour = "gray",
             linetype = "dashed") +
  scale_x_continuous(breaks = seq(-5, 5, by = 1))
```