---
title: Predict RL via HM
author: Lukas Graz
date: 2025-07-10
---

## Setup + Preprocessing

```{r}
suppressPackageStartupMessages({
  library(dplyr)
  library(ggplot2)
  library(kernlab)
  library(mgcv)  # for GAM models - smoothing
  library(gridExtra)  # for combining plots
})
```

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

D$RL_NDVI     <- pmax(0, D$RL_NDVI)
D_trn$RL_NDVI <- pmax(0, D_trn$RL_NDVI)
D_tst$RL_NDVI <- pmax(0, D_tst$RL_NDVI)

D$LANG <- as.factor(D$LANG)
D_trn$LANG <- as.factor(D_trn$LANG)
D_tst$LANG <- as.factor(D_tst$LANG)

D$SEX <- as.factor(D$SEX)
D_trn$SEX <- as.factor(D_trn$SEX)
D_tst$SEX <- as.factor(D_tst$SEX)

D$JNYTIME <- D$JNYTIME + quantile(D$JNYTIME[D$JNYTIME > 0], 0.05, na.rm=TRUE)/2
D_trn$JNYTIME <- D_trn$JNYTIME + quantile(D_trn$JNYTIME[D_trn$JNYTIME > 0], 0.05, na.rm=TRUE)/2
D_tst$JNYTIME <- D_tst$JNYTIME + quantile(D_tst$JNYTIME[D_tst$JNYTIME > 0], 0.05, na.rm=TRUE)/2

D$DISTKM <- D$DISTKM + quantile(D$DISTKM[D$DISTKM > 0], 0.05, na.rm=TRUE)/2
D_trn$DISTKM <- D_trn$DISTKM + quantile(D_trn$DISTKM[D_trn$DISTKM > 0], 0.05, na.rm=TRUE)/2
D_tst$DISTKM <- D_tst$DISTKM + quantile(D_tst$DISTKM[D_tst$DISTKM > 0], 0.05, na.rm=TRUE)/2


D$SPEED_log <- log(D$DISTKM / D$JNYTIME)
D_trn$SPEED_log <- log(D_trn$DISTKM / D_trn$JNYTIME)
D_tst$SPEED_log <- log(D_tst$DISTKM / D_tst$JNYTIME)


# logical_vars <- c("ALONE","WITH_DOG","WITH_KID","WITH_PAR","WITH_PNT","WITH_FND")
# D[logical_vars] <- lapply(D[logical_vars], as.numeric)
# D_trn[logical_vars] <- lapply(D_trn[logical_vars], as.numeric)
# D_tst[logical_vars] <- lapply(D_tst[logical_vars], as.numeric)


```


```{r}
#| code-fold: true
nams <- names(D)
nams[grep("HM|RL", nams)]

vars <- c(
  # "HM_NOISELVL",
  # "HM_COORDX","HM_COORDY","RL_COORDX","RL_COORDY","RL_GCOORD",
  # "RL_GCOORDN","RL_GCOORDW",
  "HM_NDVI","HM_NOISE","RL_NDVI","RL_NOISE",
  "ALONE","WITH_DOG","WITH_KID",
  "WITH_PAR","WITH_PNT","WITH_FND",
  "LANG","AGE","SEX",
  "DISTKM_sqrt", "JNYTIME_sqrt", "SPEED_log"
) 


summary(D[vars])



D_trn <- D_trn[vars]
D_tst <- D_tst[vars]

D_trn <- D_trn[complete.cases(D_trn), ]
D_tst <- D_tst[complete.cases(D_tst), ]

D_trn[] <- lapply(D_trn[], \(x) if (is.numeric(x)) scale(x) else x)
D_tst[] <- lapply(D_tst[], \(x) if (is.numeric(x)) scale(x) else x)
```

<!-- 
```{r}
#| include: false
#| eval: false
#| echo: false
imp_vars <- grep("^RL", vars, value=TRUE, invert=TRUE)
D_trn[imp_vars] <- xfun::cache_rds({
  missForest((D_trn[imp_vars]))
  }, 
  file = "RL-HM.rds", 
  dir = "cache/",
  hash = list(as.matrix(D_trn[imp_vars]))
)$ximp |> as.data.frame()


Res3 <- list()
for (formula_str in L_formulas) {
  # intercept_model <- lm(as.formula(paste0(
  #   mediator, " ~ 1")), D_trn)
  # step_model <- step(intercept_model, 
  #   scope = formula_str,
  #   trace = FALSE, k = log(nrow(D_trn))
  # )

}
(ResSum3 <- lapply(Res3, summary))
``` 

-->


```{r}
cor(cbind(D_trn[c("DISTKM_sqrt", "JNYTIME_sqrt", "SPEED_log")]))
```

Remove DIST_sqrt, as it is highly correlated with the other two.

## RL_NDVI
```{r}
#| label: lm-summary-rl-ndvi
lm_ndvi <- lm(RL_NDVI ~ (HM_NDVI + HM_NOISE + 
  #  ALONE + WITH_DOG + WITH_KID + WITH_PAR + WITH_PNT + WITH_FND +
   LANG + AGE + SEX +
   SPEED_log + JNYTIME_sqrt)^2, D_trn)
step_ndvi <- step(lm_ndvi, trace = FALSE, k = log(nrow(D_trn)))
summary(fit <- lm(formula(step_ndvi), D_tst))
```

- R² = 0.08
- Higher HM_NDVI corresponds to slightly higher RL-NDVI
- higher JNYTIME_sqrt corresponds to slightly higher RL-NDVI
- The faster (or further) you travel to RL, the more the RL_NDVI differs from HM_NDVI (negative interaction effect)

```{r}
ggplot(D_tst, aes(x = JNYTIME_sqrt, y=HM_NOISE, col = RL_NOISE)) +
  geom_jitter(width=0.07, height = 0.1)
```

## RL_NOISE
```{r}
#| label: lm-summary-rl-noise
lm_noise <- lm(RL_NOISE ~ (HM_NDVI + HM_NOISE + 
  #  ALONE + WITH_DOG + WITH_KID + WITH_PAR + WITH_PNT + WITH_FND +
   LANG + AGE + SEX +
   SPEED_log + JNYTIME_sqrt)^2, D_trn)
step_noise <- step(lm_noise, trace = FALSE, k = log(nrow(D_trn)))
summary(lm(formula(step_noise), D_tst))
```

- R² = 0.184
- Participants can't completely escape HM_NOISE (HM_NOISE positive predictor)
- LANGItalians have it louder (than LANG de/fr)
- Longer JNYTIME_sqrt leads to lower NOISE


```{r}
#| code-fold: true
#| message: false


# Fit Gaussian Process
X <- as.matrix(D_tst[, c("JNYTIME_sqrt", "HM_NOISE")])
y <- D_tst$RL_NOISE

# Fit GP with RBF kernel
gp_model <- gausspr(X, y, kernel = "rbfdot", kpar = "automatic")


# Create prediction grid
x_range <- range(D_tst$JNYTIME_sqrt)
y_range <- range(D_tst$HM_NOISE)
grid <- expand.grid(
  JNYTIME_sqrt = seq(x_range[1], x_range[2], length.out = 50),
  HM_NOISE = seq(y_range[1], y_range[2], length.out = 50)
)

# Predict on grid
grid$predicted_RL_NOISE <- predict(gp_model, as.matrix(grid[, 1:2]))

# Get combined range for both scales
combined_range <- range(c(D_tst$RL_NOISE, grid$predicted_RL_NOISE))
```


```{r}
#| label: fig-predicted-rl-noise-w-gp
#| code-fold: true
#| message: false

# Plot with matching color scales
ggplot() +
  geom_raster(data = grid, aes(x = JNYTIME_sqrt, y = HM_NOISE, fill = predicted_RL_NOISE)) +
  geom_jitter(data = D_tst, aes(x = JNYTIME_sqrt, y = HM_NOISE, col = RL_NOISE, shape = LANG), 
              width = 0.07, height = 0.1, alpha = 0.7) +
  scale_fill_viridis_c(name = "Predicted\nRL_NOISE", limits = combined_range) +
  scale_color_viridis_c(name = "Actual\nRL_NOISE", limits = combined_range) +
  theme_minimal()
```





```{r}
D$PRS <- rowMeans(cbind(D$FA,D$BA,D$EC,D$ES))

summary_fun <- function(x) c(mean = mean(x, na.rm = TRUE), 
                              sd = sd(x, na.rm = TRUE), 
                              median = median(x, na.rm = TRUE))

aggregate(PRS~ HM_NOISELVL, D, summary_fun)
aggregate(FA ~ HM_NOISELVL, D, summary_fun)
aggregate(BA ~ HM_NOISELVL, D, summary_fun)
aggregate(EC ~ HM_NOISELVL, D, summary_fun)
aggregate(ES ~ HM_NOISELVL, D, summary_fun)
```




## Plots for soundscape description


```{r}
#| label: fig-soundscape-description
#| code-fold: true
#| warning: false
#| fig-height: 8
#| fig-width: 7

# Define colors
colors <- c("#E31A1C", "#1F78B4", "#33A02C", "#FF7F00", "#6A3D9A", 
           "#B15928", "#A6CEE3", "#B2DF8A", "#FDBF6F")

# Create Plot 1: Dominating Sounds
suppressWarnings({
p1 <- ggplot(D, aes(x = HM_NOISE)) +
  geom_smooth(aes(y = LDOMAUD1, color = "Nature sounds"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LDOMAUD2, color = "Human sounds"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LDOMAUD3, color = "Traffic sounds"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LDOMAUD4, color = "Other technical noises"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  scale_color_manual(values = colors[1:4]) +
  labs(
    x = "Home Noise Level (HM_NOISE)",
    y = "Response Value",
    color = "Sound Type",
    title = "Dominating Sounds"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    legend.title = element_text(size = 11, face = "bold"),
    legend.text = element_text(size = 9),
    axis.title = element_text(size = 11, face = "bold"),
    panel.grid.minor = element_blank()
  )
})

# Create Plot 2: Soundscape Attributes
p2 <- ggplot(D, aes(x = HM_NOISE)) +
  geom_smooth(aes(y = LSOUNDS1, color = "Pleasant"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSOUNDS2, color = "Chaotic"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSOUNDS3, color = "Vibrant"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSOUNDS4, color = "Uneventful"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSOUNDS5, color = "Tranquil"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSOUNDS6, color = "Bothering"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSOUNDS7, color = "Eventful"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSOUNDS8, color = "Monotone"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSOUNDS9, color = "Loud"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  scale_color_manual(values = colors[1:9]) +
  labs(
    x = "Home Noise Level (HM_NOISE)",
    y = "Response Value",
    color = "Attribute",
    title = "Soundscape Attributes"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    legend.title = element_text(size = 11, face = "bold"),
    legend.text = element_text(size = 9),
    axis.title = element_text(size = 11, face = "bold"),
    panel.grid.minor = element_blank()
  )

# For the noise annoyance variables, we need to handle the character values
# Convert character annoyance values to numeric
convert_annoyance <- function(x) {
  case_when(
    x == "00" ~ 0,
    x == "01" ~ 1,
    x == "02" ~ 2,
    x == "03" ~ 3,
    x == "04" ~ 4,
    x == "05" ~ 5,
    x == "06" ~ 6,
    x == "07" ~ 7,
    x == "08" ~ 8,
    x == "09" ~ 9,
    x == "10" ~ 10,
    x == "No annoyance" ~ 0,
    x == "Response scale mid-point" ~ 5,
    x == "Very annoying" ~ 10,
    TRUE ~ as.numeric(x)
  )
}


# Create a temporary data frame with converted values for noise annoyance
D_temp <- D |>
  mutate(
    LSANNOY1_num = convert_annoyance(LSANNOY1),
    LSANNOY2_num = convert_annoyance(LSANNOY2),
    LSANNOY3_num = convert_annoyance(LSANNOY3),
    LSANNOY4_num = convert_annoyance(LSANNOY4),
    LSANNOY5_num = convert_annoyance(LSANNOY5),
    LSANNOY6_num = convert_annoyance(LSANNOY6),
    LSANNOY7_num = convert_annoyance(LSANNOY7)
  )

# Create Plot 3: Noise Annoyance
p3 <- ggplot(D_temp, aes(x = HM_NOISE)) +
  geom_smooth(aes(y = LNOISE, color = "In General"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSANNOY1_num, color = "Road Traffic"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSANNOY2_num, color = "Public Transport"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSANNOY3_num, color = "Train"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSANNOY4_num, color = "Plane"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSANNOY5_num, color = "Freetime"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSANNOY6_num, color = "Music of others"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  geom_smooth(aes(y = LSANNOY7_num, color = "Works"), method = "gam", formula = y ~ s(x), se = TRUE, alpha = 0.3, size = 1.2) +
  scale_color_manual(values = colors[1:8]) +
  labs(
    x = "Home Noise Level (HM_NOISE)",
    y = "Response Value",
    color = "Annoyance Source",
    title = "Noise Annoyance"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    legend.title = element_text(size = 11, face = "bold"),
    legend.text = element_text(size = 9),
    axis.title = element_text(size = 11, face = "bold"),
    panel.grid.minor = element_blank()
  )

xl <- xlim(30, 70)

# Combine all plots vertically
suppressWarnings({
  combined_plot <- grid.arrange(p1+xl, p2+xl, p3+xl, ncol = 1, 
                             top = "GAM Smoothing Lines: Sound Variables vs Home Noise Level")

})
# Display the combined plot
# print(combined_plot)
ggsave("cache/combinded_plot.pdf", combined_plot, height=8, width = 7)
```

- [ ] one could add rugs to the plots above
