Três Otimizações para Regressões Lineares

Este projeto visa demonstrar o cálculo fundamental por trás de uma função de regressão, utilizando três métodos distintos de otimização. Cada abordagem oferece uma perspectiva única sobre como os coeficientes da regressão são determinados, permitindo uma análise comparativa de suas eficiências e aplicabilidades em diferentes cenários.

Método de Newton: Este método encontra os coeficientes que minimizam a função de custo por meio da multiplicação da matriz Hessiana (segunda derivada da função de custo) e da matriz dos gradientes (primeira derivada da função de custo). Ele é conhecido pela sua rapidez e consistência, porém, apresenta limitações em problemas de alta dimensionalidade devido ao elevado custo computacional.
Equação Normal: Trata-se de um método que utiliza uma fórmula fechada, aplicável exclusivamente a regressões lineares com fortes correlações. Contudo, requer que a matriz de variáveis independentes seja invertível. Em situações onde duas ou mais variáveis apresentam alta correlação (correlação igual a 1), o uso desse método torna-se inviável.
Gradiente Descendente: O terceiro método de otimização utilizado é o gradiente descendente, que calcula as derivadas parciais em relação à função de custo para cada parâmetro beta e ajusta esses parâmetros na direção oposta ao gradiente. O ajuste é controlado por um parâmetro de aprendizado, o alpha, que determina a velocidade de convergência.

Cada método apresenta suas vantagens e desvantagens, e a escolha do mais adequado depende da natureza do problema e dos recursos computacionais disponíveis.

Pressupostos importantes

Para este trabalho, assumo que os dados apresentam correlação linear e trabalho em dimensões reduzidas (apenas uma feature \(x\)), além de considerar que a variável \(y\) segue uma distribuição normal. Essa suposição é comum em regressões lineares tradicionais, pois a normalidade dos resíduos garante que as estimativas dos coeficientes sejam eficientes e não viesadas.

No entanto, se a variável \(y\) não apresentar uma distribuição normal, seria necessário utilizar uma abordagem mais adequada, como as regressões do tipo GLM (Modelos Lineares Generalizados), que permitem modelar diferentes distribuições de probabilidade para a variável dependente. Além disso, em cenários com dados mais complexos ou alta variabilidade, poderíamos recorrer a técnicas de regressão mais robustas, como árvores de decisão ou Support Vector Machines (SVM), que lidam melhor com relações não lineares e distribuições mais diversas.

Lembrando que para GLM é preciso otimizar a função de log verossimilhança adaptada a distribuição esperada de \(y\) e ao final da previsão, utilizar uma link function para transformar \(y\) para a distribuição desejada. Entretando, nesses modelos, otimizar por Gradiente e Newton Method pode ser complicado, dado que a função de custo não se comporta de forma côncava ou convexa, aprendetando diversos mínimos locais.

Esses métodos avançados são úteis especialmente quando as suposições do modelo linear tradicional são violadas, proporcionando maior flexibilidade e poder preditivo em contextos mais desafiadores.

dados <- read.csv('alturas-pesos.csv')
dados2 <- read.csv('Student_Performance.csv')
dados3 <- flights

# Função para calcular RMSE - visto que queremos avaliar modelo pela escala original
rmse <- function(y_real, y_predito) {
  return(sqrt(mean((y_predito - y_real)^2)))
}

# Função para calcular R²
r_squared <- function(y_real, y_predito) {
  ss_res <- sum((y_real - y_predito)^2)
  ss_tot <- sum((y_real - mean(y_real))^2)
  return(1 - (ss_res / ss_tot))
}

# Função para calcular as previsões
y_predictions <- function(b0, b1, x1) {
  return(b0 + b1 * x1)
}

Introdução Método de Newton

O método de Newton-Raphson é amplamente utilizado para maximizar a função de log-verossimilhança, minimizando a log-verossimilhança negativa. Ele ajusta iterativamente os parâmetros ( \(\theta\) ), utilizando o gradiente e a matriz Hessiana da função de custo, até que se atinja a convergência.

Função de log-verossimilhança negativa

A função de log-verossimilhança negativa, proporcional ao erro quadrático, é definida como:

\[ \mathcal{L}(\theta) = \frac{1}{2 \sigma^2} \sum_{i=1}^{n} \left( y_i - X_i \theta \right)^2 \]

Onde: - ( \(\theta\) ) são os parâmetros; - ( \(X\) ) são as variáveis preditoras; - ( \(y\) ) é a variável dependente; - ( \(\sigma^2\) ) é a variância dos resíduos.

Gradiente da log-verossimilhança

O gradiente da função de log-verossimilhança, que indica a direção de maior crescimento, é dado por:

\[ \nabla \mathcal{L}(\theta) = -\frac{X^\top (y - X \theta)}{\sigma^2} \]

Hessiana da log-verossimilhança

A matriz Hessiana, que representa a segunda derivada da função de log-verossimilhança, é calculada como:

\[ \mathcal{H}(\theta) = \frac{X^\top X}{\sigma^2} \]

Atualização dos parâmetros

No método de Newton-Raphson, a atualização dos parâmetros ( \(\theta\) ) é feita da seguinte maneira:

\[ \theta_{\text{novo}} = \theta - \mathcal{H}^{-1}(\theta) \nabla \mathcal{L}(\theta) \]

# Função de log-verossimilhança negativa (proporcional ao erro quadrático)
log_likelihood <- function(theta, X, y, sigma2) {
  residuals <- y - X %*% theta
  return(0.5 * sum(residuals^2) / sigma2)
}

# Gradiente da log-verossimilhança
gradient <- function(theta, X, y, sigma2) {
  residuals <- y - X %*% theta
  return(- (t(X) %*% residuals) / sigma2)
}

# Hessiana da log-verossimilhança
hessian <- function(X, sigma2) {
  return(t(X) %*% X / sigma2)
}

# Estimativa de sigma^2 (variância dos resíduos)
estimate_sigma2 <- function(theta, X, y) {
  residuals <- y - X %*% theta
  return(sum(residuals^2) / length(y))
}

# Método de Newton-Raphson para máxima verossimilhança
newton_raphson_ml <- function(theta, X, y, tol = 1e-6, max_iter = 100) {
  for (i in 1:max_iter) {
    sigma2 <- estimate_sigma2(theta, X, y)
    
    grad <- gradient(theta, X, y, sigma2)
    hess <- hessian(X, sigma2)
    
    theta_new <- theta - solve(hess) %*% grad  # Atualização
    
    if (max(abs(theta_new - theta)) < tol) {
      cat("Convergiu em", i, "iterações\n")
      break
    }
    
    theta <- theta_new
  }
  
  return(theta)
}


# Exemplo de uso

y_newton <- dados2$Performance.Index
X_newton <- cbind(1, dados2$Previous.Scores)

theta_init <- c(0,0)  # Chute inicial para os parâmetros

# Executando a otimização
theta_final <- newton_raphson_ml(theta_init, X_newton, y_newton)

## Convergiu em 2 iterações

print(theta_final)

##            [,1]
## [1,] -15.181799
## [2,]   1.013837

y_pred_newton <- y_predictions(theta_final[1], theta_final[2], dados2$Previous.Scores)
rmse_1 <- rmse(dados2$Performance.Index, y_pred_newton)
r2 <- r_squared(dados2$Performance.Index, y_pred_newton)*100

ggplot(dados2, aes(x = Previous.Scores, y = Performance.Index))+
  geom_hex(col=alpha(colour = 'grey90'))+
  geom_abline(intercept =theta_final[1], slope = theta_final[2], lwd = 1.5, color = 'red2')+
  annotate("text", x = 50, y = 100, 
           label = paste("y = ", round(theta_final[1], 2), 
                         " + ", round(theta_final[2], 2), " * x  | R² = ", round(r2, 2), "%"), 
           size = 3, 
           color = "black")+
  scale_fill_gradientn(colours = c("#f7fbff", "#deebf7", "#9ecae1", "#3182bd", "#08519c"))+theme_classic2()+
  labs(
    title = "Calculation of the Performance Studies by Previous Scores - Newton Method",
    subtitle = "Using the Newton method by optimizing the negative log-likelihood" ,
    x = "Previous Scores",
    y = "Performance Studies"
  )+
  theme(
    plot.title = 
      element_text(family = "sans", size = 12, face='bold'),
    plot.subtitle = 
      element_text(family = "sans", size = 8),              
    axis.title.x = 
      element_text(family = "sans", size = 9),                 
    axis.title.y = 
      element_text(family = "sans", size = 9),                 
    axis.text = 
      element_text(family = "sans", size = 10),
    legend.title = 
      element_text(size = 0),    
    legend.text = 
      element_text(family = "sans", size = 7),      
    legend.position = "bottom",                                    
    legend.box = "horizontal",
    axis.line = element_line(linewidth = 1),  # Espessura das linhas dos eixos
    axis.ticks = element_line(linewidth = 1)
  )

Introdução a Equação Normal

A equação normal é um método utilizado para calcular os coeficientes de uma regressão linear de forma direta, por meio de uma solução analítica. Esse método é especialmente eficaz quando há fortes correlações entre as variáveis preditoras, pois utiliza uma fórmula fechada para determinar os parâmetros ( \(\theta\) ).

Fórmula da Equação Normal

Os coeficientes ( \(\theta\) ) podem ser obtidos pela seguinte fórmula:

\[ \theta = (X^T X)^{-1} X^T y \]

Onde: - ( \(X\) ) é a matriz das variáveis independentes (ou preditores); - ( \(y\) ) é o vetor da variável dependente; - ( \(X^T\) ) representa a transposta da matriz ( \(X\) ); - ( \((X^T X)^{-1}\) ) é a inversa do produto ( \(X^T X\) ).

Condição de Inversibilidade

É importante notar que a inversa de ( \((X^T X)\) ) só existe se a matriz for invertível, o que implica que não deve haver correlações perfeitas entre as variáveis independentes. Se pelo menos duas variáveis são perfeitamente correlacionadas (correlação igual a 1), o método se torna inviável.

cor(flights$dep_delay, flights$arr_delay, use = "complete.obs")

## [1] 0.9148028

sapply(dados3, function(x) sum(is.nan(x)))

##           year          month            day       dep_time sched_dep_time 
##              0              0              0              0              0 
##      dep_delay       arr_time sched_arr_time      arr_delay        carrier 
##              0              0              0              0              0 
##         flight        tailnum         origin           dest       air_time 
##              0              0              0              0              0 
##       distance           hour         minute      time_hour 
##              0              0              0              0

dados3 <- dados3[complete.cases(dados3), ]
  
  # Ler os dados
y <- dados3$arr_delay
X <- cbind(1, dados3$dep_delay)  # Adiciona a coluna de 1s para o intercepto

# Calcular os coeficientes usando a equação fechada
# X_t é a transposta de X
X_t <- t(X)
# Calcula (X^T X)
X_t_X <- X_t %*% X
# Calcula (X^T X)^-1
X_t_X_inv <- solve(X_t_X)
# Calcula (X^T y)
X_t_y <- X_t %*% y
# Calcula os coeficientes
theta <- X_t_X_inv %*% X_t_y

# Exibir os resultados
intercepto <- theta[1]
coeficiente <- theta[2]

cat("Intercepto:", intercepto, "\n")

## Intercepto: -5.899493

cat("Coeficiente:", coeficiente, "\n")

## Coeficiente: 1.019093

y_pred_normal <- y_predictions(intercepto, coeficiente, dados3$dep_delay)
rmse_2 <- rmse(dados3$arr_delay, y_pred_normal)
r2 <- r_squared(dados3$arr_delay, y_pred_normal)*100

ggplot(dados3, aes(x = dep_delay, y = arr_delay))+
  geom_hex(col=alpha(colour = 'grey90'))+
  geom_abline(intercept =intercepto, slope = coeficiente, lwd = 1.5, color = 'red2')+
  scale_fill_gradientn(colours = c("#9ecae1", "#08519c"))+theme_classic2()+
  annotate("text", x = 180, y = 1500, 
           label = paste("y = ", round(intercepto, 2), " + ", 
                         round(coeficiente, 2), " * x  | R² = ", round(r2,2),"%"), 
           size = 3, 
           color = "black")+
  labs(
    title = "Calculation of the Delay in Takeoff x Landing Delay Regression  ",
    subtitle = "Using the Normal Equation that automatically minimizes the MSE." ,
    x = "Delay in Takeoff",
    y = "Landing Delay",
    color = "Frequency"
  )+
  theme(
    plot.title = 
      element_text(family = "sans", size = 12, face="bold"),
    plot.subtitle = 
      element_text(family = "sans", size = 8),              
    axis.title.x = 
      element_text(family = "sans", size = 9),                 
    axis.title.y = 
      element_text(family = "sans", size = 9),                 
    axis.text = 
      element_text(family = "sans", size = 10),
    legend.title = 
      element_text(size = 0),    
    legend.text = 
      element_text(family = "sans", size = 6),      
    legend.position = "bottom",                                    
    legend.box = "horizontal",
    axis.line = element_line(linewidth = 1),  # Espessura das linhas dos eixos
    axis.ticks = element_line(linewidth = 1)
  )

Introdução ao Gradiente Descendente

O Gradiente Descendente é um algoritmo de otimização utilizado para minimizar a função de custo em problemas de regressão. Ele ajusta iterativamente os parâmetros, calculando os gradientes da função de custo e atualizando os coeficientes na direção oposta ao gradiente, controlada por um fator de aprendizado.

Função de Custo

A função de custo, que mede a qualidade das previsões, é geralmente representada pelo Erro Quadrático Médio (MSE):

\[ J(b_0, b_1) = \frac{1}{2n} \sum_{i=1}^{n} (y_i - (b_0 + b_1 x_{1i}))^2 \]

Onde: - ( \(J\) ) é a função de custo; - ( \(b_0\) ) é o intercepto; - ( \(b_1\) ) é o coeficiente para a variável ( \(x_1\) ); - ( \(n\) ) é o número de exemplos; - ( \(y_i\) ) são os valores reais.

Cálculo dos Gradientes

Os gradientes da função de custo em relação aos parâmetros ( \(b_0\) ) e ( \(b_1\) ) são dados por:

\[ \frac{\partial J}{\partial b_0} = \frac{1}{n} \sum_{i=1}^{n} (y_i - (b_0 + b_1 x_{1i})) \]

\[ \frac{\partial J}{\partial b_1} = \frac{1}{n} \sum_{i=1}^{n} (y_i - (b_0 + b_1 x_{1i})) x_{1i} \]

Atualização dos Parâmetros

Os parâmetros são atualizados da seguinte forma:

\[ b_0 := b_0 - \alpha \frac{\partial J}{\partial b_0} \]

\[ b_1 := b_1 - \alpha \frac{\partial J}{\partial b_1} \]

Onde ( \(\alpha\) ) é a taxa de aprendizado.

Vamos aplicar agora o Gradiente Descendente para encontrar a regressão que minimize MSE pela derivada parcial a mesma, para relaciona Peso e Altura

-> Uma obervação, eu iniciei \(b0\) como -158 e \(b1\) como 1.37 dado que por conta da taxa de aprendizado \(a\) demorava muito e exigia iterações exponenciais para encontrar os parâmetros corretos que minimizem a função de regressão

# Função de Custo - MSE
cost_function <- function(b0, b1, x1, y) {
  n <- length(y) 
  predictions <- b0 + b1 * x1  # Previsões
  mse <- sum((predictions - y)^2) / (2 * n)  # Erro quadrático médio
  return(mse)
}

# Implementando o Gradiente Descendente com critério de parada
gradient_descent_l1 <- function(x1, y, alpha, iterations, lambda, epsilon) {
  b0 <- -158  # Inicializa o intercepto
  b1 <- 1.37 # Inicializa o coeficiente para x1
  n <- length(y)  # Número de exemplos
  cost_history <- numeric(iterations)
  
  # Criando vetores para armazenar o histórico de coeficientes
  b0_history <- numeric(iterations)
  b1_history <- numeric(iterations)
  cost_history <- numeric(iterations)
  
  for (i in 1:iterations) {
    predictions <- b0 + b1 * x1
    
    # Calculando os gradientes
    db0 <- (1 / n) * sum(predictions - y)
    db1 <- (1 / n) * sum((predictions - y) * x1)
    
    # Atualizando os parâmetros
    b0 <- b0 - alpha * db0
    b1 <- b1 - alpha * db1
    
    # Guardando dados em cada iteração
    b0_history[i] <- b0
    b1_history[i] <- b1
    
    # Recalculando a função de custo com os novos b0 e b1
    cost_history[i] <- cost_function(b0, b1, x1, y)
    
    # Verificando o critério de parada baseado nos gradientes
    if (abs(db0) < epsilon && abs(db1) < epsilon) {
      print(paste("Convergiu na iteração:", i))
      break
    }
    
  }
  return(list(b0 = b0, b1 = b1, cost_history = cost_history, iterations = i,
              b0_history = b0_history, b1_history = b1_history))
}

# Configurando parâmetros
alpha <- 0.00001  # Taxa de aprendizado
iterations <- 10000  # Número máximo de iterações
epsilon <- 0.01  # Critério de parada baseado nos gradientes


# Carregando e normalizando os dados
y <- dados$Peso
X <- dados$Altura


# Executando o Gradiente Descendente Matricial
gd_result <- gradient_descent_l1(X, y, alpha, iterations, lambda, epsilon)

## [1] "Convergiu na iteração: 28"

# Resultados
cat("Coeficientes encontrados (theta):", gd_result$b1, "\n")

## Coeficientes encontrados (theta): 1.372833

cat("Coeficientes encontrados (beta):", gd_result$b0, "\n")

## Coeficientes encontrados (beta): -158

y_pred_gd <- y_predictions(gd_result$b0, gd_result$b1, dados$Altura)
rmse_3 <- rmse(dados$Peso, y_pred_gd)
r3 <- r_squared(dados$Peso, y_pred_gd)*100


ggplot(dados, aes(x = Altura, y = Peso))+
  geom_hex(col=alpha(colour = 'grey90'))+
  geom_abline(intercept =gd_result$b0, slope = gd_result$b1, lwd = 1.5, color = 'red2')+
  scale_fill_gradientn(colours = c("#f7fbff", "#deebf7", "#9ecae1", "#3182bd", "#08519c"))+theme_classic2()+
  annotate("text", x = 150, y = 125, 
           label = paste("y = ", round(gd_result$b0, 2), " + ", 
                         round(gd_result$b1, 2), " * x  | R² = ", round(r2,2),"%"), 
           size = 3, 
           color = "black")+
  labs(
    title = "Calculation of the Weight [Kg] x Height [cm] Regression  ",
    subtitle = "Using the Gradient Descent that automatically minimizes the MSE by n iterations." ,
    x = "Height [cm]",
    y = "Weight [Kg]"
  )+
  theme(
    plot.title = 
      element_text(family = "sans", size = 12, face="bold"),
    plot.subtitle = 
      element_text(family = "sans", size = 8),              
    axis.title.x = 
      element_text(family = "sans", size = 9),                 
    axis.title.y = 
      element_text(family = "sans", size = 9),                 
    axis.text = 
      element_text(family = "sans", size = 10),
    legend.title = 
      element_text(size = 0),    
    legend.text = 
      element_text(family = "sans", size = 6),      
    legend.position = "bottom",                                    
    legend.box = "horizontal",
    axis.line = element_line(linewidth = 1),  # Espessura das linhas dos eixos
    axis.ticks = element_line(linewidth = 1)
  )

Avaliação RMSE dos modelos

# Cálculo das médias das colunas especificadas
media_1 <- mean(dados2$Performance.Index)  # Exemplo de cálculo da média
media_2 <- mean(dados3$arr_time)
media_3 <- mean(dados$Peso)

# Criando a tabela de métricas com colunas separadas
tabela_metricas <- data.frame(
  Modelo = c("Newton Method", "Equação Normal", "Gradiente Descendente"),
  RMSE = c(rmse_1, rmse_2, rmse_3),
  Média = c(media_1, media_2, media_3)
)

# Exibindo a tabela formatada com colunas separadas
kable(tabela_metricas, col.names = c("Modelo", "RMSE", "Média"),
      caption = "Métricas de Performance dos Modelos")

Métricas de Performance dos Modelos
Modelo	RMSE	Média
Newton Method	7.742746	55.2248
Equação Normal	18.027398	1501.9082
Gradiente Descendente	5.568620	72.7339

Dada a forte linearidade entre as variáveis e a simplicidade do problema (com apenas uma dimensão), os modelos convergiram eficientemente, resultando em boas soluções. Os coeficientes de determinação (\(R²\)) também foram satisfatórios, confirmando a qualidade dos modelos ajustados.

No entanto, é importante destacar que, à medida que a dimensionalidade do problema aumenta, torna-se necessário considerar o uso de regularizações, como \(L2\) (Ridge) ou \(L1\) (Lasso), nas funções de custo. Lembre-se de que, ao calcular o gradiente, o fazemos com base na função de custo, o que implica que a regularização também será derivada e incorporada nas iterações do processo de otimização.

Por fim, vale mencionar que este estudo não tem como objetivo realizar tratamentos extensivos de dados ou gerar insights profundos. O foco aqui é demonstrar diferentes abordagens para regressões lineares, além do tradicional método de Mínimos Quadrados Ordinários (OLS), ao qual estamos mais habituados.

Manual Regression Calculate

Luis Henrique

2024-09-20