---
title: "Logistic regression"
author: "Julian Faraway"
date: "`r format(Sys.time(), '%d %B %Y')`"
output:
  html_document
---

```{r echo=FALSE}
library(knitr)
opts_chunk$set(comment=NA, warning=FALSE, fig.path="/tmp/Figs/")
```

Logistic regression is useful for modeling a binary response.

We start by loading the `tidyverse` and the data.

```{r}
library(tidyverse)
load("hospdata.rda")
```

In our example dataset, gender is the only binary variable. It does not make sense to predict this
so we create another binary:

```{r}
hospdata %>% mutate(pneumonia = factor(ifelse(str_sub(Illness,1,3) == "J18","yes","no"))) -> hospdata
hospdata
```

The `pneumonia` variable indicates whether the patients was admitted for this condition as `yes`
and `no` if not.

A *factor* variable is a special kind of character variable that R likes to use in modeling. It is more convenient
for some purpose - for example, we ask for a summary:

```{r}
summary(hospdata)
```

which comes out nicely. But sometimes factor variables cause problems elsewhere so it is better to stick with character variables by default.

The summary confirms that `actividad` now has only two levels. There are methods that can model non-binary responses but these are more complicated.

# Fitting the model

Logistic regression is an example of a *Generalized Linear Model* (GLM) - it generalizes the normally distributed response found in the linear model to other kinds of response. This is done by specifying the *family* - in this case, binomial:

```{r}
gmod = glm(pneumonia ~ Gender + Age + PM10, data=hospdata, family = binomial())
```

The logistic model takes the form:
```
log(p/(1-p)) = beta0 + beta1 * x1 + beta2 * x2 + ... + error
```

It models the probability *p* of the second level of the response (which is `yes`). The term 
*log(p/(1-p))* is called the *log-odds*. Odds are an alternative representation of probability often used
in gambling. It is more convenient to work with odds so we exponeniate:

```{r}
exp(coef(gmod))
```

We can see that the odds of having pneumonia are multiplied by 1.21 if the patient is male. In other words, men are more likely to get pneumonia than women.  We see that older people are more likely to have pneumonia and relatively more cases will be observed when the PM10 is higher.

## Exercises

Make `Gender` the response and `pneumonia + Age + PM10` the predictors. (This is not practical but do
it for the sake of the exercise). You will need to make `Gender` a factor with

```
hospdata %>%  mutate(Gender = factor(Gender)) -> hospdata
```

1. Does the model predict the probability of being a man or a woman?

2. What is the effect of `pneumonia` on the odds of the gender of the patient?

# Hypothesis testing

We can examine the model output

```{r}
summary(gmod)
```

and see that both Gender and PM10 are not statistically significant. We might propose a
model only using the age:

```{r}
gmodp = glm(pneumonia ~ Age, data=hospdata, family = binomial())
```

We would like to compare this with the original model. Unfortunately, there is a complication because of the missing values. We construct a data set that only uses the variables in the models and only uses complete cases:

```{r}
hospdata %>% select(pneumonia, Age, PM10, Gender) %>% na.omit() -> hospdata2
```

We refit both models on the complete data.

```{r}
gmod = glm(pneumonia ~ Gender + Age + PM10, data=hospdata2, family = binomial())
gmodp = glm(pneumonia ~ Age, data=hospdata2, family = binomial())
```


We can compare this to the original model with:

```{r}
anova(gmodp,gmod,test="Chi")
```

Since the p-value is large, we pick the smaller model (which only uses Age as a predictor)



## Exercises

Use hypothesis testing to reduce the model of the previous exercise. Which variables are selected? If you want to fit a model with
no predictors, use the formula `Gender ~ 1`.

# Prediction

We can make predictions. For example, we predict the outcome for a patient, aged 50:

```{r}
predict(gmodp, data.frame(Age=50))
```

The answer is on the log-odds scale which may not be useful. We ask for a probability instead:

```{r}
predict(gmodp, data.frame(Age=50), type="response")
```

We see a 15% chance of pneumonia.

We can make several predictions at the same time. Let's predict the probability of being occupied
for all ages between 20 and 80 both for women and men:

```{r}
pv = tibble(Age=0:80,pred=predict(gmodp,tibble(Age=0:80),type="response"))
pv
```

We can make a plot of these probabilities:

```{r}
ggplot(pv, aes(x=Age,y=pred)) + geom_line() + ylab("P(pneumonia)")
```

We can see how the probability of having pneumonia increases with Age.


## Exercises

Using the model, `Gender ~ pneumonia + Age + PM10`, predict the probability of being a woman given Age=50, PM10 = 100, pneumonia="yes".



