mvn_lm.Rmd

---
title: "MVN + Regression Interactions"
author: "Eoin Travers"
date: "`r Sys.Date()`"
output: html_document
---

https://twitter.com/realHollanders/status/1605462357697667072?s=20&t=sfGV7iiPavt6DwCrMdks-w

```{r}
library(tidyverse)
library(MASS)
```

```{r}
n = 200
X = mvrnorm(n, mu = c(0, 0), Sigma = S, empirical = T)
x1 = X[,1]; x2 = X[,2]

# Fixed parameters
b1 = 1
b2 = 1
b12 = -.8
y_err = .1

simulate = function(b12, y_err){
  y = x1 * b1 + x2 * b2 + x1 * x2 * b12
  if(y_err > 0) y = y + rnorm(n, 0, y_err)
  df = data.frame(x1, x2, y)
  cor(df)
}

cor_obs = simulate(b12, y_err)
cor_obs
```

```{r}
evaluate = function(b12, y_err = globalenv()$y_err){
  cor_pred = simulate(b12, y_err)
  sum((cor_obs - cor_pred)^2 )
}

b12_vals = seq(-1, 1.01, .01)
loss_vals = map_dbl(b12_vals, evaluate, y_err = y_err) # NB - Is robust to providing wrong y_err

b12_estimate = optimize(evaluate, c(-2, 2), y_err = y_err)$minimum
b12_estimate
```


```{r}
plot(b12_vals, loss_vals, 'l',
     xlab = 'Interaction term', ylab = 'Loss')
abline(v = b12)
abline(v = b12_estimate, lty = 'dashed')
```