Let’s see if we can roll our own weighted linear regression using the mle function and our own definition of likelihood.



Data to play with

set.seed(1)
data <- data.frame(
  x = 1:10,
  # Exponentiated error to thwart the normalcy 
  # assumption.
  y = 1:10 + exp(rnorm(10))  
)



R’s native linear regression model

model <- lm(y ~ x, 
            data=data,
            weights = 1:10)
coefficients(summary(model))
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    2.083      1.171    1.78 0.113085
x              0.923      0.158    5.85 0.000382



Homebrewed regression

library(stats4)
negative_log_likelihood <- function(i, x) {
  w <- 1:10 / sum(1:10)
  n <- nrow(data)
  df <- length(start)
  prediction <- x * data$x + i
  e <- data$y - prediction
  s2 <- sum(w * e^2)
  
  ((n - df)/2) * (log(2*pi) + log(s2) + 1)
}
start <- list(
  i = 0,
  x = 0
)
model <- mle(negative_log_likelihood, start)
summary(model)@coef
  Estimate Std. Error
i    2.083      1.171
x    0.923      0.158



Yes we can.