---
title: 'Lab #10: Updating a prediction model'
output:
  html_document:
    df_print: paged
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
set.seed(24)
#install.packages("rms")
library(rms)
library(ranger)
library(ROCit)
```

# Logistic recalibration of a prediction model

1. Create a function `make_data` that simulates data according to a logistic model. The function accepts the number of observations `n`, the true logistic model coefficients `beta`, and the true intercept `intercept`. Generate the X matrix where the variables are independently distributed, and the outcomes according to the logistic model. Return a data frame.
```{r}
make_data <- function(n, beta, intercept=0){
  p <- length(beta)
  x <- matrix(rnorm(n * p), ncol = p)
  beta_vec <- matrix(beta, ncol=1)
  logit <- x %*% beta_vec + intercept
  mu <- 1/(1 + exp(-logit))
  y <- rbinom(n=n, size=1, prob=mu)
  data.frame(x=x, y=y)
}

make_data(n=10, beta=c(1,1))
```

2. Create 100 training observations and create a 10-dimensional vector `beta` where the first 3 coefficients are -1, the next 3 are 1, and the last 4 are zero. Set the intercept to zero.
```{r}
beta <- c(-1,-1,-1,1,1,1,0,0,0,0)
data_original <- make_data(n=100, beta=beta, intercept=0)
data_original
```


3. Let's fit a random forest using `ranger` that predicts probabilities and call it `rf`. Use the default options.
```{r}
rf <- ranger(as.factor(y) ~ ., data=data_original, probability = T)
rf
```

4. Create a test dataset with 1000 observations. What is the AUC of `rf`?
```{r}
test_original <- make_data(n=1000, beta=beta, intercept=0)
test_pred <- predict(rf, test_original)$predictions[,2]
roc_original <- rocit(test_pred, test_original$y)
plot(roc_original)
ciAUC(roc_original)
```

5. Plot the calibration curve for the random forest using the `val.prob` function from the `rms` package. Do not do any recalibration just yet, so set the argument `logistic.cal=F`.
```{r}
val_res <- val.prob(c(test_pred), test_original$y, logistic.cal = F, g = 10)
```

6. Let's fit a logistic model to recalibrate the scores predicted from the random forest. In fact, `val_res` returns the slope and intercept from running recalibration. Use this to recalibrate the probabilities from `rf`. What does the calibration curve look like now?
```{r}
recalib_log_odds <- val_res["Intercept"] + val_res["Slope"] * log(test_pred/(1 - test_pred))
recalib_prob <- 1/(1 + exp(-recalib_log_odds))
val_res <- val.prob(recalib_prob, test_original$y, logistic.cal = F, g = 10)
```

# Adjusting the model for a new site

7. Let's generate data from a new site. Suppose the associations are similar but different in this population. Define a new set of coefficients for this new site, `beta2`, where the first 3 coefficients are now -0.9 instead. Also change the intercept to 0.1. Create 40 new observations from this new site and call it `new_data`.

```{r}
beta2 <- c(rep(-0.9, 3), beta[-seq(3)])
beta2
new_data <- make_data(n=40, beta=beta2, intercept = 0.1)
new_data
```

8. What is the AUC of `rf` in the new population? Has it changed?
```{r}
new_pred <- predict(rf, new_data)$predictions[,2]
roc_new <- rocit(new_pred, new_data$y)
ciAUC(roc_new)
```

9. Plot the calibration curve of the original model in this new population, again with `logistic.cal=F`. How does it look?
```{r}
val.prob(new_pred, new_data$y, logistic.cal = F, g=5)
```

10. Let's try recalibrating the model using logistic regression. Create a recalibration dataset and fit a logistic model. Make sure the logistic regression model is applied to the predicted logit, not the probability
```{r}
recalib_df <- data.frame(
  logit=log(new_pred/(1 - new_pred)),
  y=new_data$y
)
recalib_glm <- glm(y ~ ., data=recalib_df, family=binomial, offset=logit)
summary(recalib_glm)
```

11. Based on the above results, would you advocate using logistic recalibration?
```{r}
```

12. For comparison, let's try refitting a random forest only using only the data available at the new site. Name the new model `new_rf`.
```{r}
new_rf <- ranger(as.factor(y) ~ ., data=new_data, probability = T)
```

13. Generate a new batch of data at the new site with 200 observations. Use `val.prob` to evaluate the AUC and calibration of `new_rf` on `new_data_big`. Compare it to the recalibrated version of `rf`.
```{r}
new_data_big <- make_data(n=500, beta=beta2, intercept=0.1)

new_rf_pred <- predict(new_rf, new_data_big)$predictions[,2]
val.prob(new_rf_pred, new_data_big$y, logistic.cal = F, g=10)

rf_pred <- predict(rf, new_data_big)$predictions[,2]
rf_pred_df <- data.frame(logit=log(rf_pred/(1 - rf_pred)))
recalib_rf_prob <- predict(recalib_glm, rf_pred_df, type="response")
val.prob(recalib_rf_prob, new_data_big$y, logistic.cal = F, g=10)
```