```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,tidy=TRUE,message=FALSE,warning=FALSE)
showsol<-TRUE
```


# Interactive lecture - first week

## Theoretical questions - with and without use of R

### Problem 1: Model assumptions

a) What are the model assumptions for a binary regression? 
b) Which link function and response function is used for the logit model? 
c) What is the difference between the logit model and a logistic regression?

---

### Problem 2: Log-likelihood.  

a) What is the definition of the log-likelihood?

b) For the logit model the log-likelihood is
$$l(\boldsymbol{\beta})=\sum_{j=1}^G[\ln \binom{n_j}{y_j}+ y_j \ln \pi_j-y_j\ln(1-\pi_j)+n_j\ln(1-\pi_j)]$$
for grouped data. Explain how we have arrived at this formula?

c) Write the version of the loglikelihood for individual data (i.e. $n_j=1$ and $G=n$).

d) Where is $\boldsymbol{\beta}$ in the loglikelihood in c? Rewrite this to be a function of $\boldsymbol{\beta}$.  

e) Why can we ignore the normalizing constant (what is the constant?) in the case of $n_j = 1 \ \forall j$? Considering what the log-likelihood is used for, why can we ignore the normalizing constant in all cases (i.e., also when $n_j \neq 1$)?

f) What does this graph of $l$ look like as a function of $\boldsymbol{\beta}$ for the beetle data? First discuss shortly and then to aid you in answering this we look at the loglikelihood for the beetle data. Read the R code, discuss what is done and work on interpreting the final graph - in particular comment on the yellow ridge in the plot.

The beetle data has only one covariate (in addition to the intercept) - so this means that we have $\boldsymbol{\beta}=(\boldsymbol{\beta}_0,\boldsymbol{\beta}_1)$. Look at the following code and explain what is done - remark: we have used the $n_i=1$ version of the loglikelihood here.

```{r}
library(investr)
library(ggplot2)
library(viridis)

# from aggregated to individual data (because these data were aggregated)
ldose <- rep(investr::beetle$ldose, investr::beetle$n)
y <- NULL
for (i in 1:8) y = c(y, rep(0, investr::beetle$n[i] - investr::beetle$y[i]), rep(1, investr::beetle$y[i]))
beetleds = data.frame(killed = y, ldose = ldose)

loglik <- function(par, args){
  y <- args$y; x <- args$x; n <- args$n
  res <- sum(y*x%*%par - n*log(1 + exp(x%*%par)))
  return(res)
}

loglik(c(1,1), args = list(y = beetleds$killed, 
                           x = cbind(rep(1, nrow(beetleds)), beetleds$ldose), 
                           n = rep(1, nrow(beetleds))))

loglikmat <- matrix(NA, nrow = 100, ncol = 100)
loglikframe <- data.frame()
beta_0 <- seq(-90,-30, length.out = 100)
beta_1 <- seq(20, 50, length.out = 100)

for (i in 1:length(beta_0)){
  for (j in 1:length(beta_1)){
    
    loglikmat[i,j] <- loglik(c(beta_0[i], beta_1[j]), args = list(y = beetleds$killed, 
                                                                  x = cbind(rep(1, nrow(beetleds)), beetleds$ldose), 
                                                                  n = rep(1, nrow(beetleds))))
    
    loglikframe <- rbind(loglikframe, c(beta_0[i], beta_1[j], loglikmat[i,j]))
  
  }
}
names(loglikframe) <- c("beta_0", "beta_1", "loglik")
head(loglikframe)

ggplot(data = loglikframe, mapping = aes(x = beta_0, y = beta_1, z = loglik)) + geom_raster(aes(fill = exp(0.0001*loglik))) +
  geom_point(data = loglikframe[which.max(loglikframe$loglik),], mapping = aes(x = beta_0, y = beta_1), 
             size = 5, col = "red", shape = 21, stroke = 2) + scale_shape(solid = FALSE) +
  scale_fill_viridis() + geom_contour(col = "black")
```

Comments to the code: for the `loglik` function we have two arguments: par= the parameters to be estimated, and args=a list with data. The reason for only having these two arguments is that it is easier to use when we later perform optimization (with `optim`) of the loglikelihood to find the ML estimates.

---

### Problem 3: Score function

a) What is the definition of the score function? What is the dimension of the score function?
b) Derive the score function for the logit model (individual data). The result should be
$$s(\boldsymbol{\beta})=\sum_{i=1}^n {\bf x}_i (y_i-\pi_i)=\sum_{i=1}^n {\bf x}_i (y_i-\frac{\exp({\bf x}_i^T\boldsymbol{\beta})}{1+\exp({\bf x}_i^T\boldsymbol{\beta})})$$
c) What do we need the score function for? 

---

### Problem 4: Fisher information. 

(We did not cover this in the lecture week 1, but we know one of the definitions from Module 2. Either you skip Problem 4 and move to Problem 5, or you look at the section ["Properties of the score function"](#propscore), and ["The expected Fisher information matrix $F(\boldsymbol{\beta})$"](#covscore) together.)

a) What is the definition of the expected (and the observed) Fisher information matrix? What is the dimension of thise matrix (matrices)?
b) What is the role of these matrices in ML estimation?

c) For the logit model with grouped data the expected and the observed Fisher information matrix are equal and given as

$$F(\boldsymbol{\beta})=\sum_{j=1}^G {\bf x}_j {\bf x}_j^T n_j \pi_j (1-\pi_j)$$
Where is $\boldsymbol{\beta}$ in this expression?

d) Write the version of the expected Fisher information for individual data (i.e. $n_j=1$ and $G=n$). 

---

### Problem 5: Maximum likelihood

To find the ML estimate for $\boldsymbol{\beta}$ we may either use the function `glm` or optimize the log-likelihood manually. We will do both.

a) First we use the `glm` function in R, and we also check that the individual and the grouped data give the same parameter estimates for the $\boldsymbol{\beta}$. Read the R-code, notice the different input structures and check the results.

```{r}
# the beetle.ds was made above
fitind=glm(killed ~ ldose, family = "binomial", data = beetleds) # individual data
summary(fitind)
```

```{r}
fitgrouped=glm(cbind(y, n-y) ~ ldose, family = "binomial", data = investr::beetle) # grouped data. response is #success AND #fails (here we have defined a dead beetle as a success)
summary(fitgrouped)
```

b) What is the default convergence criterion for the glm? (Note: IRWLS used in `glm` - more in Module 5.)

c) We implemented the log-likelihood as a function in item 2 above. Now we will use this together with the `optim` function on the beetle data set to optimize the loglikelihood. Read the R-code, take notice of how we put data into the `args` slot and how the optimization is called with `optim`. (In Compulsory exercise 2 you will use this in a Poisson regression.)

```{r}
loglik_gr <- function(par, args) {
    
    y <- args$y
    x <- args$x
    n <- args$n
    
    res <- y %*% x - t(t(n * x) %*% ((1 + exp(-x %*% par))^(-1)))
    return(res)
}

opt <- optim(c(-60, 30), fn = loglik, gr = loglik_gr, args = list(y = beetleds$killed, 
    x = cbind(rep(1, nrow(beetleds)), beetleds$ldose), n = rep(1, nrow(beetleds))), 
    control = list(fnscale = -1), hessian = TRUE,method="BFGS")
opt 

sqrt(diag(-solve(opt$hessian))) # calculate the standard deviations of the parameters
```

---

### Problem 6: Interpreting results

a) Interpret the estimated $\boldsymbol{\beta}$´s. Odds ratio is useful for this. 
b) Plot the predicted probability of a beetle dying against the dosage and discuss what you see. (Yes, since this is the last question you may try to program by yourself!)

<!-- # Summing up and moving on? -->

<!-- ## Aim of binary regression:  -->
<!-- 1. Construct a model to help understand the relationship between a "success probability" and one or several explanatory variables. The response measurements are binary (present/absent, true/false, healthy/diseased). -->
<!-- 2. Use the model for estimation and prediction of success probabilites. -->

<!-- Our two examples: -->

<!-- * probability of killing a beetle, for given dose chemicals -->

<!-- * probability for getting infant respiratory disease, covariates are sex of infant (boy, girl) and method of feeding (bottle, supplement, breast) -->

---


# References for further reading

* A. Agresti (2015): "Foundations of Linear and Generalized Linear Models." Wiley.
* A. J. Dobson and A. G. Barnett (2008): "An Introduction to Generalized Linear Models", Third edition. 
* J. Faraway (2015): "Extending the Linear Model with R", Second Edition. <http://www.maths.bath.ac.uk/~jjf23/ELM/>
* P. McCullagh and J. A. Nelder (1989): "Generalized Linear Models". Second edition.


