# Lihelihood vs density ---- 

# Likelihood is not a probability function
set.seed(7)
y <- rnorm(5, mean = 1, sd = 1)
lik <- function(theta)
  sapply(theta,
         function(t) prod(dnorm(y, mean = t, sd = 1)))

curve(dnorm(x, 1,1), xlim = c(-100,100))
# Area under the density in y (theta fixed)
integrate(dnorm, -100, 100, mean = 1, sd = 1)$value

# Area under the likelihood in theta (y fixed)
integrate(lik, -100, 100)$value


# Plot of the likelihood for the sample of size n = 5
curve(lik(x), xlim = c(-2, 4), lwd = 2, col = "blue",
      xlab = expression(theta), ylab = expression(L(theta)),
      main = "Likelihood for a sample of size n = 5")
abline(v = mean(y), lty = 2, col = "red")

# Plots of likelihood, score, observed information matrix, and asymptotic distribution ----

# Maximum likelihood and score: N(mu, 1), n = 5 

set.seed(7)
y <- rnorm(5, mean = 1, sd = 1)
n <- length(y)
mu_hat <- mean(y)

loglik <- function(mu, y)
  sapply(mu,
         function(t) sum(dnorm(y, mean = t, sd = 1, log = TRUE)))
score  <- function(mu, y) length(y) * (mean(y) - mu)

# Log-likelihood (maximum and tangent lines) and score
theta_hat <- mean(y)
par(mfrow = c(1, 2), mar = c(3.5, 3.5, 2, 1), mgp = c(2, 0.6, 0))

# Left: log-likelihood, maximum and tangent lines at two values of theta
curve(loglik(x, y), xlim = c(-2, 4), lwd = 2, col = "blue",
      xlab = expression(theta), ylab = expression(l(theta)),
      main = "Log-likelihood")
abline(v = theta_hat, lty = 2, col = "red")
abline(h = loglik(theta_hat, y), lty = 3, col = "red")
points(theta_hat, loglik(theta_hat, y), pch = 19, col = "red", cex = 1.5)
text(theta_hat + 0.1, loglik(theta_hat, y) - 3, expression(hat(theta)), col = "red", pos = 4)
th0 <- c(-0.5, 2.2); cl <- c("darkgreen", "darkorange")
for (j in 1:2) {
  xs <- th0[j] + c(-0.7, 0.7)
  lines(xs, loglik(th0[j], y) + score(th0[j], y) * (xs - th0[j]), col = cl[j], lwd = 2)
  points(th0[j], loglik(th0[j], y), pch = 19, col = cl[j])
}

# Right: score function, zero at theta_hat
curve(score(x, y), xlim = c(-2, 4), lwd = 2, col = "blue",
      xlab = expression(theta), ylab = expression(l*"'"*(theta)),
      main = "Score")
abline(h = 0, lty = 3); abline(v = theta_hat, lty = 2, col = "red")
points(theta_hat, 0, pch = 19, col = "red", cex = 1.5)
points(th0, score(th0, y), pch = 19, col = cl)


# Observed information and (asymptotic) distribution 
set.seed(11)
y50 <- rnorm(50, mean = 1, sd = 1)

# observed information by numerical second derivative
h <- 1e-3
jhat <- function(y) {
  th <- mean(y)
  -(loglik(th + h, y) - 2 * loglik(th, y) + loglik(th - h, y)) / h^2
}
jhat(y); jhat(y50)   # equal to n

par(mfrow = c(1, 2), mar = c(3.5, 3.5, 2, 1), mgp = c(2, 0.6, 0))

# Left: relative log-likelihood, n = 5 vs n = 50
curve(loglik(mean(y) + x, y) - loglik(mean(y), y), xlim = c(-1.5, 1.5), ylim = c(-12, 0),
      lwd = 2, col = "blue", xlab = expression(theta - hat(theta)),
      ylab = expression(l(theta) - l(hat(theta))), main = "observed info")
curve(loglik(mean(y50) + x, y50) - loglik(mean(y50), y50), add = TRUE, lwd = 2, col = "darkred")
legend("bottomleft", legend = c("n = 5", "n = 50"), col = c("blue", "darkred"), lwd = 2, bty = "n")

# Right: asymptotic N(theta, 1/j) distribution of the estimator
curve(dnorm(x, mean = 1, sd = sqrt(1 / jhat(y))), xlim = c(-0.5, 2.5), lwd = 2, col = "blue",
      ylim = c(0, 3), xlab = expression(theta), ylab = "density",
      main = expression("Distribution of " * hat(theta)))
curve(dnorm(x, mean = 1, sd = sqrt(1 / jhat(y50))), add = TRUE, lwd = 2, col = "darkred")
abline(v = 1, lty = 3)
points(c(mean(y), mean(y50)), c(0, 0), pch = 19, col = c("blue", "darkred"), cex = 1.5)

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1, mgp = c(3, 1, 0))   # back to default


# numerical optimisation ----

# Example 
# To do togheter
# 1) write down the (negative) likelihood function 
# 2) Let's generate a sample of size $n = 50$ from a $\mathcal{N}(\mu = 10, \sigma^2 = 4)$
# 3) Get the MLE using the nlminb() optimizer
# 4) compare with the analytic solution

nllik <- function(par, data){
  mu <- par[1]
  sigma2 <- par[2]
  n <- length(data)
  out <- n/2 * log(sigma2) + 1/(2*sigma2) * sum((data-mu)^2)
  return(out)
}

set.seed(3)
n <- 50 
y <- rnorm(n, mean = 10, sd = 2) 

par0 <- c(1,2)

nllik(par0, y) # check 

mle1 <- nlminb(start  = par0, 
               objective = nllik, 
               data = y, 
               lower=c(-Inf, 1e-9))

str(mle1)

nllik(c(9.87,3.1), y)

mean(y)
var(y)*(length(y)-1)/length(y)

# Exercise ----
# 1) Adapt the negative log-likelihood function built above to the case of a simple linear regression model 
# 2) Let's generate some data: 
#    * set $n=100$ 
#    * generate $x$ from a uniform distribution in (0,10) 
#    * fix $\beta_1 = 2$, $\beta_2 = 0.8$, and $\sigma = 1.5$
# 3) fit the model using nlminb
# 4) then, compare the results with those obtained via lm function 

nllik_lm <- function(par, y, x){
  beta_1 <- par[1]
  beta_2 <- par[2]
  sigma2 <- par[3]
  mu <- beta_1 + beta_2*x
  n <- length(y)
  out <- n/2 * log(sigma2) + 1/(2*sigma2) * sum((y-mu)^2)
  return(out)
}

set.seed(3)
n <- 100
x <- runif(n, 0, 10)
beta_1 <- 2 
beta_2 <- 0.8
sigma2 <- (1.5)^2

y <- beta_1 + beta_2*x + rnorm(n, 0, sqrt(sigma2))

plot(x,y)
abline(c(2,0.8), col = "red")


fit_lm <- nlminb(start = c(0,0,1),
                 objective = nllik_lm ,
                 x = x, 
                 y = y, 
                 lower = c(-Inf, -Inf, 1e-9))
                 
str(fit_lm)             
               

lm1 <- lm(y ~ x)                    
summary(lm(y ~ x))

sum(residuals(lm1)^2)/length(y)


# Inference on a single coefficient ----

data("Wage")
lmFit <- lm(wage ~ age + year + education, data = Wage)
summary(lmFit)$coefficients

# Some useful quantities
X <- model.matrix(wage ~ age + year + education, data = Wage)
n <- nrow(X)
p <- ncol(X)
hat_s2 <- sum(residuals(lmFit)^2)/(nrow(Wage) - p)

# Confidence intervals: confint() and by hand
confint(lmFit)
cbind(lmFit$coefficients - sqrt(hat_s2 * diag(solve(t(X) %*% X)))* qt(0.975, df = n-p),
      lmFit$coefficients + sqrt(hat_s2 * diag(solve(t(X) %*% X)))* qt(0.975, df = n-p))

# Hypothesis tests: summary() and by hand
summary(lmFit)$coefficients[,3:4]

tval <- lmFit$coefficients/(sqrt(hat_s2 * diag(solve(t(X) %*% X))))

cbind(tval, 2*pt(-abs(tval), n - p))


# Inference on a pool of coefficients ---- 

hat_sigma2 <- sum(residuals(lmFit)^2)/n
tilde_sigma2 <- var(Wage$wage)*(n-1)/n

fval <- ((tilde_sigma2 - hat_sigma2)/(p-1))/((hat_sigma2)/(n-p))
fval
pf(fval, p-1, n-p, lower = FALSE)

summary(lmFit)
summary(lmFit)$fstatistic


# Inference on a subset of coefficients (education) ----

boxplot(wage ~ education, data = Wage,
        names = c("<HS", "HS", "Some coll.", "College", "Advanced"),
        xlab = "education", ylab = "wage", col = "lightblue")

# Full and reduced model
fitFull <- lm(wage ~ age + year + education, data = Wage)
fitRed  <- lm(wage ~ age + year, data = Wage)   # without education

round(summary(fitFull)$coefficients, 3)

# The test by hand
n_obs  <- nrow(Wage)
p_full <- length(coef(fitFull))
q      <- p_full - length(coef(fitRed))
c(n = n_obs, p = p_full, q = q)

hat_sigma2   <- sum(residuals(fitFull)^2)/n_obs
tilde_sigma2 <- sum(residuals(fitRed)^2)/n_obs

Fobs <- ((tilde_sigma2 - hat_sigma2)/q)/(hat_sigma2/(n_obs - p_full))
Fobs
pf(Fobs, q, n_obs - p_full, lower.tail = FALSE)

# The test with anova(): reduced model, full model
anova(fitRed, fitFull)