---
title: Prediction Analysis for WSL
author: Lukas Graz
date: 2025-02-13
---

```{r}
#| code-fold: false
source("R/data_prep.R")
```

## Setup
```{r}
#| code-fold: true
#| code-summary: "Benchmark helper function"
library(mlr3verse, quietly = TRUE)

mse <- msrs(c("regr.mse"))

if (!interactive())
  lgr::get_logger("mlr3")$set_threshold("warn")

get_benchi_table <- function(tasks, nfolds = 5) {
  set.seed(123)
  learners <- lrns(c("regr.featureless", "regr.lm", "regr.xgboost", "regr.ranger"))
  learners$regr.xgboost$param_set$set_values(
    eta = 0.03, 
    nrounds = 300, 
    max_depth = 2
  )

  benchi <- xfun::cache_rds({
    benchmark(benchmark_grid(
      tasks, 
      learners, 
      rsmp("cv", folds = nfolds)
    ))
  }, 
  file = "benchmark.rds", 
  dir = "cache/",
  hash = list(tasks, nfolds)
  )
  
  res <- tidyr::pivot_wider(benchi$aggregate(mse), 
    id_cols = task_id,
    names_from = learner_id,
    values_from = regr.mse
  ) |> as.data.frame()
  
  rownames(res) <- res$task_id
  res <- res[, -1]
  colnames(res) <- gsub("regr.", "", colnames(res))
  stopifnot(any(colnames(res) == "featureless"))
  res <- 1 - res / res$featureless
  res[, -1, drop = FALSE] |> round(3)
}
```

Testing prediction quality of GIS_vars -> Mediators -> PRS_vars using 

- Linear models 
- Random forests (default parameters)
- XGBoost (with parameter tuning)
- FASSO (not shown since inferior)

**GIS Variables:**

```{r}
#| echo: false
GIS_vars
```

**Mediators:**

```{r}
#| echo: false
Mediator_vars
```

### PRS ~ GIS
```{r}
#| code-fold: true
tasks_GIS <- lapply(PRS_vars, \(y) 
  as_task_regr(
    subset(Dmlr, select = c(y, GIS_vars)),
    target = y,
    id = y
  ))
get_benchi_table(tasks_GIS) 
```
GIS variables alone show poor predictive performance.

### PRS ~ GIS + Mediators
```{r}
#| code-fold: true
tasks_GIS_MED <- lapply(PRS_vars, \(y) 
  as_task_regr(
    subset(Dmlr, select = c(y, Mediator_vars, GIS_vars)),
    target = y,
    id = y
  ))
get_benchi_table(tasks_GIS_MED) 
```

### PRS ~ Mediators
```{r}
#| code-fold: true
tasks_MED <- lapply(PRS_vars, \(y) 
  as_task_regr(
    subset(Dmlr, select = c(y, Mediator_vars)),
    target = y,
    id = y
  ))
get_benchi_table(tasks_MED) 
```

### Mediators ~ GIS
```{r}
#| code-fold: true
tasks_MED_by_GIS <- lapply(Mediator_vars, \(y) 
  as_task_regr(
    subset(Dmlr, select = c(y, GIS_vars)),
    target = y,
    id = y
  ))
get_benchi_table(tasks_MED_by_GIS)
```

### Legacy Code
```{r}
#| code-summary: "XGBoost - Parameter Tuning"
#| code-fold: true
#| eval: false

# Get parameter estimates for XGBoost
t <- as_task_regr(
  subset(Dmlr, select = c("FEELNAT", GIS_vars)),
  target = "FEELNAT"
)

l <- lrn("regr.xgboost",
  nrounds = 500  # More iterations due to lower learning rate
)

# Create search space
ps <- ps(
  max_depth = p_int(2, 3),
  eta = p_dbl(0.001, 0.3, tags = "logscale")
)

# Setup tuning
instance <- ti(
  task = t,
  learner = l,
  resampling = rsmp("cv", folds = 3),
  measure = msr("regr.mse"),
  terminator = trm("none"),
  search_space = ps
)

# Grid search
tuner <- mlr3tuning::tnr("grid_search")
tuner$optimize(instance)
instance$result
```

```{r}
#| code-fold: true
#| code-summary: "Superseded setup with mlr3"
#| echo: false
#| eval: false
#| message: false
#| warning: false

library(randomForest)

fit <- lm(as.formula(paste0(
    "cbind(", paste(PRS_vars, collapse = ", "), ")",
    " ~ ",
    paste(Mediator_vars, collapse = " + ")
  )), 
  D)
coef(fit) |> round(2)

rsq.lm <- sapply(summary(fit), \(x) x$r.sq)
rsq.rf <- sapply(PRS_vars, \(x) {
  rf <- randomForest(as.formula(paste0(
    x, " ~ ", paste(Mediator_vars, collapse = " + ")
  )), 
  D, na.action = na.omit
  ) 
  rf$rsq[500]
})

cbind(lm = rsq.lm, rf = rsq.rf) |> round(2)
```

```{r}
#| code-fold: true
#| code-summary: "First basic setup with mlr3"
#| echo: false
#| message: false
#| warning: false
#| eval: false

autoplot(mytsk1, type = "pairs")
mytsk1 <- as_task_regr(
  subset(Dmlr, select = c("FA", Mediator_vars, GIS_vars)),
  feature = c(Mediator_vars, GIS_vars),
  target = "FA",
  id = "bla"
)

lrn_xgb <- lrn("regr.xgboost")
lrn_avg <- lrn("regr.featureless")
splits <- partition(mytsk1)
lrn_xgb$train(mytsk1, splits$train)$predict(mytsk1, splits$test)$score(mse)
lrn_avg$train(mytsk1, splits$train)$predict(mytsk1, splits$test)$score(mse)
rr <- resample(mytsk1, lrn_xgb, cv3)
rr$aggregate(mse)

learners <- lrns(c("regr.featureless", "regr.lm", "regr.xgboost", "regr.ranger"))
learners$regr.xgboost$param_set$set_values(eta = 0.03, nrounds = 300, max_depth = 2)
learners <- c("regr.featureless", "regr.lm")
```
