Contour plot of likelihood under multicollinearity

multico.R

library(ggplot2)
library(dplyr)
library(mvtnorm)
set.seed(1)
x1 = rnorm(10, 0, 1)
x2 = x1+rnorm(10, 0, 0.1)
y = x1+x2 + rnorm(10, 0, 1)

print(cor(x1,x2))

fit1 = lm(y~x1+x2-1)
summary(fit1)

rmse = function(b1,b2,y,x1,x2){
    mean(dnorm(y, b1*x1+b2*x2, 1, log=TRUE))
}

dfB = expand.grid(b1 = seq(-15, 15, length.out = 200),
                  b2 = seq(-15, 15, length.out = 200)) %>%
  rowwise() %>%
  mutate(lp = rmse(b1,b2,y,x1,x2))

ggplot(dfB, aes(x=b1, y=b2, fill=exp(lp)))+
  geom_tile() + 
  theme_bw()

summa1 = summary(fit1)

rand = data.frame(rmvnorm(10000, fit1$coefficients,summa1$cov.unscaled))

ggplot(dfB)+
  geom_point(data=rand, aes(x=x1, y=x2), alpha=0.1, size=0.1) +
  geom_contour(aes(x=b1, y=b2, z=exp(lp), colour=after_stat(level))) + 
  scale_colour_viridis_c()+
  theme_bw(16)+labs(x="coef1", y="coef2")

ggsave("contor.png")