Tutorial explicativo em R para o cálculo do Mínimos Quadrados Ordinários (MQO) em uma regressão linear univariada.

O método mais utilizado em regressão linear para estimação dos parâmetros de regressão é o Método dos Mínimos Quadrados Ordinários (MQO), que tem como princípio a minimização da soma do quadrado dos desvios dos valores observados em relação à reta de regressão.

O objetivo deste tutorial é explicar de maneira objetiva, simplificada e replicável como utilizar esta técnica para produzir suas regressões lineares.

Ao final deste tutorial, mostro como fazer o calculo nativamente pelo R.

Minimização

Minimizar algo é atingir um valor cada vez menor. Se você é um pedreiro construindo uma casa, e seu contratante pediu para minimizar os custos da obra, significa que você encontrará o ajuste perfeito da quantidade de areia, cimento, pedra e água que compõem a massa para erguer as paredes. Você quer o mínimo de cimento que forneça uma parede ideal, forte e robusta.

Tomando a analogia de uma bola rolando em uma supercie tal qual uma montanha, esta irá rolar até encontrar o ponto mais baixo desta função. Este é o ponto de mínimo da superfície da montanha.

O mesmo funciona para funções matemáticas. Minimizar uma função é encontrar o menor valor que satisfaça determinadas condições.

Valores Observados

Valores Observados são dados coletados durante a realização de um experimento ou provindo de registros históricos. A cotação de ativos na bolsa de valores, por exemplo, é um exemplo de uma série de dados ordenada no tempo. Um pesquisador que sai na floresta da UFSCar para coletar dados do tamanho de flores de uma espécie também esta coletando dados.

Geralmente, quando se quer avaliar o efeito explicativo de uma variável na outra, utiliza-se técnicas estatisticas, a variável dependente é comumente chamada de Y e a explicativa de X. Um exemplo prático disto é fornecido abaixo, no qual a Nota da Prova (Y) depende, em certa medida, das horas de estudo (X) empregadas.

Espera-se, neste exemplo fictício, que quanto maior as horas de estudo, maior seja o impacto na nota obtida.

Os pares ordenados (Xi, Yi) representam as observações individuais no conjunto de dados e são usados para traçar o gráfico de dispersão.

# Dataframes genéricos contendo as variáveis observadas
horas_estudo <- c(2, 3, 4, 5, 6)
notas_prova <- c(60, 70, 80, 90, 95)
dados <- data.frame(Xi = horas_estudo, Yi = notas_prova)
dados
##   Xi Yi
## 1  2 60
## 2  3 70
## 3  4 80
## 4  5 90
## 5  6 95
plot(dados$Xi, dados$Yi, main = "Gráfico de dispersão", xlab = "horas de estudo", ylab = "notas obtidas")

Desvios

Em estatística, o termo desvio refere-se a uma subtração em relação a outra variável. O Desvio aconte em relalção a alguma medida estabelecida anterioremente. Existe o desvio padrão, o desvio absoluto e no caso do MQO estamos nos referindo ao desvio em relação à média das variáveis oservadas. Seguindo o exemplo do dataframe dado:

# Caluclando a média das horas de estudo e das notas da prova
MediaY <- mean(notas_prova)
MediaX <- mean(horas_estudo) 

# Calculando o desvio das horas de estudo 
DesvioX <- horas_estudo - MediaX
DesvioY <- notas_prova - MediaY

# Juntando os desvios calculados no dataframe
dados <- cbind(dados, DesvioY = DesvioY, DevioX = DesvioX)
dados
##   Xi Yi DesvioY DevioX
## 1  2 60     -19     -2
## 2  3 70      -9     -1
## 3  4 80       1      0
## 4  5 90      11      1
## 5  6 95      16      2

Aplicando o método MQO

A partir do método dos mínimos quadrados ordinários é possível encontrar os coeficientes B1 e β0 de um regressão linear.

β1 é o coeficiente angular da reta de regressão, que indica a mudança esperada na variável dependente (Y) para cada unidade de mudança na variável independente (X).

Em outras palavras é a inclinação da reta que se pretender traçar que melhor se ajusta a cada par de pontos Xi e Yi. Quanto maior o for a relação de dependência entre Y e X, maior será o ajuste da reta traçada em relação ao gráfico de dispersão. A geração desta reta nada mais é do que o modelo de regressão linear simples.

Matemáticamente,

β1 = Σ(produto dos desvios médios de XiYi) / Σ(dos quadrados dos desviosde Xi),

Em outras palavras:

Para obter o numerador, multiplica-se cada elmento da coluna DesvioX pelo seu par correspondente na coluna DevioY. Aplica-se o somatório aos produtos resultantes.

Para obter o denominador, elevam-se os desvios médios de X ao quadrado. Multiplica-se cada elemento da coluna DesvioX por ele mesmo. Aplica-se o somatório aos produtos resultantes.

# Obtendo o produto dos desvios médios
XiYi <- DesvioY *DesvioX

# Adicionando ao dataframe

dados <- cbind(dados, XiYI = XiYi)
dados
##   Xi Yi DesvioY DevioX XiYI
## 1  2 60     -19     -2   38
## 2  3 70      -9     -1    9
## 3  4 80       1      0    0
## 4  5 90      11      1   11
## 5  6 95      16      2   32
# Calculando o somatório = numerador 
Numerador <- sum(XiYi)
Numerador
## [1] 90
# Obtendo o quadrado dos desvios médios de X
desvioXiq <- DesvioX^2

# Adicionando ao dataframe
dados <- cbind(dados, DesvioXiq = desvioXiq)
dados
##   Xi Yi DesvioY DevioX XiYI DesvioXiq
## 1  2 60     -19     -2   38         4
## 2  3 70      -9     -1    9         1
## 3  4 80       1      0    0         0
## 4  5 90      11      1   11         1
## 5  6 95      16      2   32         4
# calculando o somatório = denominador
Denominador <- sum(desvioXiq)
Denominador
## [1] 10
β1 <- Numerador/Denominador
β1
## [1] 9
# pode-se sintetizar a expressão matemática como: 
beta1 <- sum(DesvioY *DesvioX)/sum(DesvioX^2)
beta1 # todavia não seria de fácil compreensão para o leitor.
## [1] 9

β0 é o valor da interceptação da linha de regressão com o eixo y. Em outras palavras, é o valor estimado da variável dependente (Y) quando a variável independente (X) é igual a zero.

O coeficiente linear é uma medida importante na análise de regressão, pois indica o ponto de partida da reta de regressão no eixo y.

Matemáticamente, podemos obtêlo da seguinte forma:

β0 = y(medio) - β1 * x(médio)

# Calculo do β0
β0 <- MediaY - β1 * MediaX
β0
## [1] 43

Por fim, podemos plotar a equação de reta que descreve a regressão linear obtida pelo método dos mínimos quadrados ordinários, sendo esta que melhor se adequa ao conjuto de dados obtidos.

Y = β0 + β1X

Y = 43 + 9X

Por fim, basta calcular os valores de Y esperados pas os mesmos valores de X passados no dataframe inicial. Em outras palavras, seremos capazes de prever a variável dependente Y de acordo com a variável explicativa X.

# Assumindo Xi como pontos iniciais do dataframe dados tal que Xi = horas_estudo

Ye <- 43 + 9*horas_estudo

# adicionando ao dataframe

dados <- cbind(dados, Ye = Ye)
dados
##   Xi Yi DesvioY DevioX XiYI DesvioXiq Ye
## 1  2 60     -19     -2   38         4 61
## 2  3 70      -9     -1    9         1 70
## 3  4 80       1      0    0         0 79
## 4  5 90      11      1   11         1 88
## 5  6 95      16      2   32         4 97
# Faz o gráfico da reta prevista pelo MQO junto com o gráfico orginal de dispersão. 
plot(dados$Xi, dados$Yi, main = "Gráfico de dispersão", xlab = "horas de estudo", ylab = "notas obtidas")
abline(a=β0, b = β1, col = "red", lwd = 2)

Utilizando o R para calcular MQO por você

Sabendo que os resultados de interesse do MQO são os coeficientes β0 e β1, podemos utilizar as funções nativas do R.

# em uma única linha fazemos todo o tutorial acima. 
reglin <-lm(Yi ~ Xi, data = dados)

# com a função summary percebemos o intercepto β0 = 43 e β1 = 9  
summary(reglin)
## 
## Call:
## lm(formula = Yi ~ Xi, data = dados)
## 
## Residuals:
##          1          2          3          4          5 
## -1.000e+00  8.937e-15  1.000e+00  2.000e+00 -2.000e+00 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  43.0000     2.4495   17.55 0.000403 ***
## Xi            9.0000     0.5774   15.59 0.000574 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.826 on 3 degrees of freedom
## Multiple R-squared:  0.9878, Adjusted R-squared:  0.9837 
## F-statistic:   243 on 1 and 3 DF,  p-value: 0.0005737
# outra forma de chamar os coeficientes descrito é por meio do objeto coef()
coef(reglin)
## (Intercept)          Xi 
##          43           9