---
title: Hypothesis Testing
author: Lukas Graz
date: 2025-02-13
include-in-header:
  - text: |
      <style>
      .quarto-notebook .cell-container .cell-decorator {
        display: none;
      }
      </style>
---

```{r}
#| code-fold: true
source("R/data_prep.R")
options(digits = 3)

# interactive <- function() FALSE
# D_trn$HM_NOISE_nrm <- scale(D_trn$HM_NOISE)
# D_tst$HM_NOISE_nrm <- scale(D_tst$HM_NOISE)
## No suspicious patterns in missing data
# mice::md.pattern(D[c(Mediator_vars, GIS_vars)[nNAs>0]], plot = FALSE)
# mice::md.pattern(D[PRS_orig_vars])
```

## Linear Modeling
### Imputation with MissForest on Training Data
```{r}
#| code-fold: true
#| code-summary: "Number of NAs in Mediators and GIS variables"
sapply(D[Mediator_vars], \(x) sum(is.na(x)))
sapply(D[GIS_vars], \(x) sum(is.na(x)))
```

```{r}
#| code-fold: true
#| code-summary: "Impute missing values using MissForest"
# Mediator imputation
D_trn[Mediator_vars] <- xfun::cache_rds({
  missForest(as.matrix(D_trn[Mediator_vars]))
  }, 
  file = "Mediator_imputation.rds", 
  dir = "cache/",
  hash = list(as.matrix(D_trn[Mediator_vars]))
)$ximp |> as.data.frame()

# GIS imputation (missForest)
D_trn[GIS_vars] <- xfun::cache_rds({
  missForest(as.matrix(D_trn[GIS_vars]))
  }, 
  file = "GIS_imputation.rds", 
  dir = "cache/",
  hash = list(as.matrix(D_trn[GIS_vars]))
)$ximp |> as.data.frame()
```


### Scaling Test Data
```{r}
#| code-fold: true
#| code-summary: "Scaling variables and show old scale"
all_vars <- c(Mediator_vars, GIS_vars, PRS_vars)
old_scale <- t(sapply(D_tst[c(all_vars)], \(x) 
  c(mean = mean(x, na.rm = TRUE), sd = sd(x, na.rm = TRUE))))

D_tst[c(all_vars)] <- lapply(D_tst[c(all_vars)], scale)

old_scale
```


### Testing VIF
```{r}
#| code-fold: true
#| code-summary: "VIF: PRS ~ Mediators + GIS_vars  (NO interaction)"
car::vif(fit_PRS_MED <- lm(as.formula(paste0(
  PRS_vars[1], 
  " ~ ", 
  paste(Mediator_vars, collapse = " + "), " + ",
  paste(GIS_vars,      collapse = " + ")
)), D_trn)) |> summary()
```

```{r}
#| code-fold: true
#| code-summary: "VIF: PRS ~ (Mediators + GIS_vars)^2 (WITH interaction)"
suppressMessages(
car::vif(fit_PRS_MED <- lm(as.formula(paste0(
  PRS_vars[1], 
  " ~ ", 
  "(", paste(Mediator_vars, collapse = " + "), 
     " + ", paste(GIS_vars, collapse = " + "), 
  ")^2"
)), D_trn))) |> summary()
```

Since we model _with_ interactions later, the latter VIF are relevant for us. Given that they are very high (c.f. median and max), we would have no hope of finding any significant results in the full interaction model. Therefore, we will first perfom a variable selection, to reduce the VIF and enable us to find significant effects.


### All Interactions: Mediators ~ (GIS)^2
```{r}
Res3 <- list()
for (mediator in Mediator_vars) {
  intercept_model <- lm(as.formula(paste0(
    mediator, " ~ 1")), D_trn)
  step_model <- step(intercept_model, 
    scope = as.formula(paste0(
      mediator, " ~ ", 
      "(", paste(GIS_vars, collapse = " + "), ")^2"
    )),
    trace = FALSE, k = log(nrow(D_trn))
  )
  Res3[[mediator]] <- lm(formula(step_model), D_tst)
}
lapply(Res3, summary)
```

### All Interactions: PRS ~ (Mediators + GIS)^2
```{r}
Res4 <- list()
for (prs in PRS_vars) {
  intercept_model <- lm(as.formula(paste0(
    prs, " ~ 1")), D_trn)
  step_model <- step(intercept_model, 
    scope = as.formula(paste0(
      prs, " ~ ", 
      "(", paste(GIS_vars, collapse = " + "), " + ", 
      paste(Mediator_vars, collapse = " + "), ")^2"
    )),
    trace = FALSE, k = log(nrow(D_trn))
  )
  Res4[[prs]] <- lm(formula(step_model), D_tst)
}
lapply(Res4, summary)
```

### Legacy Code
```{r}
#| code-fold: true
#| code-summary: "Variable selection over multiple y (not desired)"
#| echo: false
#| eval: false

# library(grpreg)
# fit <- cv.grpreg(X = model.matrix(fit_MED_GIS)[,-1], y = D_trn[Mediator_vars[1:2]])
# coef(fit) |> t()
# fit$beta
# plot(fit)
```

```{r}
#| code-fold: true
#| code-summary: "Linear models with mice"
#| echo: false
#| eval: false

Y <- D[]

library(mice, quietly = TRUE)
library(car, quietly = TRUE)
library(miceadds, quietly = TRUE)
data(nhanes2, package = "mice")
set.seed(9090)

mi.res <- miceadds::mice.1chain(nhanes2, burnin = 4, iter = 20, Nimp = 8)
an2a <- miceadds::mi.anova(mi.res = mi.res, formula = "bmi ~ age * chl")

mod1 <- with(mi.res, stats::lm(bmi ~ age * chl))
mod0 <- with(mi.res, stats::lm(bmi ~ age + chl))

mitml::testModels(model = mod1$analyses, null.model = mod0$analyses, method = "D1")
mitml::testModels(model = mod1$analyses, null.model = mod0$analyses, method = "D2")

an2b <- miceadds::mi.anova(mi.res = mi.res, formula = "bmi ~ age * chl", type = 3)
```


```{r}
#| code-fold: true
#| code-summary: "Mediators ~ GIS"
#| echo: false
#| eval: false
stop("this should not run")
Res1 <- list()
for (mediator in Mediator_vars) {
  full_model <- lm(as.formula(paste0(
    mediator, " ~ ", 
    "HM_NOISE_nrm * (", paste(GIS_vars, collapse = " + "), ")"
  )), D_trn)
  small_model <- step(full_model, trace = FALSE, k = log(nrow(D_trn)))
  Res1[[mediator]] <- lm(formula(small_model), D_tst)
}
lapply(Res1, summary)
```

```{r}
#| code-fold: true
#| code-summary: "PRS ~ Mediators"
#| echo: false
#| eval: false
Res2 <- list()
for (mediator in Mediator_vars) {
  full_model <- lm(as.formula(paste0(
    mediator, " ~ ", 
    "HM_NOISE_nrm * (", paste(GIS_vars, collapse = " + "), ")"
  )), D_trn)
  small_model <- step(full_model, trace = FALSE, k = log(nrow(D_trn)))
  Res2[[mediator]] <- lm(formula(small_model), D_tst)
}
lapply(Res2, summary)
```

