###############################################################################
# Title: Optimisation of a Linear Regression Model in R
# Date:  2020-05-26
#
# Author: Joshua P. Entrop
# Website: joshua-entrop.com
###############################################################################

# 1. Prefix -------------------------------------------------------------------

# 1.1 Clear working space
rm(list = ls())

# 1.2 Loading the example data set icu from the package aplore3 ---------------
library(aplore3)
y = icu$sys                 #Set our depended variable
x1 = icu$age                #Set our fist independed variable
x2 = as.numeric(icu$gender) #Set our second independed variable

# 2. Define the liklihood function we like to optimise ------------------------
ll_lm <- function(par, y, x1, x2){
  
  alpha <- par[1]
  beta1 <- par[2]
  beta2 <- par[3]
  sigma <- par[4]
  
  R = y - alpha - beta1 * x1 - beta2 * x2
  
  -sum(dnorm(R, mean = 0, sigma, log = TRUE))
}

# 3. Define initials values ---------------------------------------------------
est_alpha <- mean(icu$sys)
est_beta1 <- mean(icu$sys[icu$age >= 40 & icu$age <= 41]) - mean(icu$sys[icu$age >= 41 & icu$age <= 52])
est_beta2 <- mean(icu$sys[icu$gender == "Male"]) - mean(icu$sys[icu$gender == "Female"])
est_sigma <- sd(icu$sys)

# 4. Use optim() minimise the likelihood function -----------------------------
mle_par <- optim(fn = ll_lm,                #Function to be optimised
                 par = c(alpha = est_alpha, #Initial values
                         beta1 = est_beta1, 
                         beta2 = est_beta2, 
                         sigma = est_sigma), 
                 y = icu$sys,               #Empirical data from the data set icu
                 x1 = icu$age,
                 x2 = as.numeric(icu$gender))

mle_par$par                                 #Showing estimates for unknown parameters

# 5. Compare results to estimates from lm() -----------------------------------
summary(lm(sys ~ age + as.numeric(gender), data = icu))

#' ////////////////////////////////////////////////////////////////////////////
#' END OF R-FILE