Simple linear models
A regressão linear impõe ou presume que exista uma relação linear entre uma ou mais variáveis preditoras \(x_i\) e uma variável resposta \(y\).
\[y_i \sim Normal(\mu_i, \sigma) \\ \mu_i = \alpha + \beta x_i\]
A primeira linha desse modelo descreve a distribuição de probabilidade da variável \(y\). Ela pode ser modelada como uma variável com uma distribuição normal, com uma média \(\mu\) e um desvio-padrão \(\sigma\).
A segunda linha do modelo relaciona as variáveis preditoras X com a média esperada μ (mu), baseada numa relação causal hipotetizada entre X e Y. Nesse caso simples, a relação é linear e pode ser representada por uma linha que é dada pelos parâmetros de intercepto α (alfa) e de inclinação β (beta).
Gerando um conjunto de dados.
Para esse tutorial introdutório, vamos utilizar dados simulados, de forma que conheçamos a relação real entre as variáveis X e Y.
N = 30
x <- runif(N, 0, 5)
n <- length(x)
mu <- 1.2 + 3.5 * x # α = 1.2 e β = 3.5.
y <- rnorm(N, mean = mu, sd = 3)
lm_data <- data.frame(x, mu, y, res = y - mu)Um gráfico simples para os dados:
plot(y ~ x, pch = 19)
Como podemos ver no gráfico, parece haver uma relação linear simples entre as duas variáveis, como previsto pela simulação.
Coeficiente de regressão
Podemos estimar o valor dos parâmetros alfa e beta usando muitos métodos diferentes. O método mais comum é o Ordinary Least Squares (OLS), que produz soluções analíticas para o valor desses dois coeficientes.
O valor do coeficiente beta é dado pela covariância entre X e Y, dividido pela variância de X. E o intercepto alpha é dado pela média de Y menos o produto do valor de beta pela média de X.
\[\widehat{\beta}= \frac{Cov(x, y)}{Var(x)} = \frac{\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\]
and
\[ \hat{\alpha} = \bar{y} - \hat{\beta}\bar{x} \]
[1] 3.601063[1] 3.601063[1] 1.284795Usando esse valor dos parâmetros, podemos incluir a retro-regressão no gráfico que fizemos acima.
Usando a função lm() para ajustar o modelo
Claro que não precisamos calcular os parâmetros diretamente. Podemos usar uma função do R para ajustar esse modelo linear. Essa função é a função “lm” e pode ser chamada usando a sintaxe de fórmula do R.
m1 = lm(y ~ x)
m1
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
1.285 3.601 Com o modelo ajustado, podemos usar a função “summary” para ver as estimativas dos coeficientes, que devem ser idênticas às que obtivemos acima, e também a sua significância e várias outras informações úteis.
summary(m1)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-6.6792 -2.9851 0.2627 2.6821 5.4490
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2848 1.2934 0.993 0.329
x 3.6011 0.4267 8.438 3.54e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.414 on 28 degrees of freedom
Multiple R-squared: 0.7178, Adjusted R-squared: 0.7077
F-statistic: 71.21 on 1 and 28 DF, p-value: 3.545e-09Resíduos da regressão
Uma vez que o modelo está ajustado, podemos calcular os seus resíduos (\(e_i\)). Ao invés de tratarmos os resíduos simplesmente como “erros”, vamos tratar esses valores como a variação residual não explicada pela variável preditora.
Matematicamente, o resíduo para a observação \(i\) é a diferença entre o valor real observado (\(y_i\)) e o valor predito pelo nosso modelo linear (\(\mu_i\)):
\[e_i = y_i - \mu_i\]
Podemos extrair essa variação não explicada diretamente do nosso modelo ajustado no R usando funções nativas:
Para entender o poder dessa extração, precisamos visualizar os resíduos de duas perspectivas complementares. A primeira delas é a perspectiva geométrica, pela qual enxergamos o resíduo como a distância vertical entre o dado observado e a reta de regressão.
# 1. Visualizando os resíduos como distância da reta
plot(y ~ x, pch = 19, main = "Perspectiva 1: Distância da reta de regressão")
abline(m1, col = "red", lwd = 2)
# Desenhando linhas conectando os pontos reais aos valores preditos
segments(x0 = x, y0 = y, x1 = x, y1 = mu, col = "blue", lty = 2)
A segunda perspectiva é muito mais poderosa para nossos usos: podemos isolar essas distâncias (as linhas azuis tracejadas do gráfico acima) e tratá-las como uma nova variável: os valores de y descontados os efeitos de x.
# 2. Visualizando os resíduos como uma variável independente
plot(res ~ x, pch = 19, col = "blue",
main = "Perspectiva 2: A variação residual como variável",
ylab = "Variação residual não explicada")
abline(h = 0, col = "red", lwd = 2, lty = 2)
Compare os dois gráficos. O segundo gráfico é essencialmente o primeiro, mas nós “rotacionamos” a linha de regressão até que ela ficasse horizontal, centralizada no zero.
Entender os resíduos como uma variável própria, limpa da influência de um preditor específico, é o que nos permitirá usar a regressão como uma ferramenta de controle estatístico. Mais à frente neste tutorial, usaremos essa propriedade para descontar nossas variáveis de interesse de fontes de variação espúria.
Variáveis preditoras discretas
Caso a variável preditora X não seja contínua, o modelo linear pode ser usado da mesma maneira. Vamos simular dados em que a variável resposta Y represente o peso de um indivíduo e a variável X seja o estado nutricional: bem alimentado ou mal alimentado.
set.seed(1)
# 1. Criar a variável preditora X (Categórica)
n <- 100
estado_nutricional <- sample(c("bem_alimentado", "mal_alimentado"), n, replace = TRUE)
peso_medio_por_classe = c("bem_alimentado" = 75, "mal_alimentado" = 60)
peso <- peso_medio_por_classe[estado_nutricional] + rnorm(n, mean = 0, sd = 5)
# Criar o dataframe
dados <- data.frame(peso, estado_nutricional)
head(dados) peso estado_nutricional
1 76.99053 bem_alimentado
2 56.93987 mal_alimentado
3 76.70560 bem_alimentado
4 69.35318 bem_alimentado
5 67.16512 mal_alimentado
6 84.90200 bem_alimentadoAgora vamos ajustar um modelo linear e tentar capturar a diferença nas médias nessas duas classes.
Call:
lm(formula = peso ~ estado_nutricional)
Residuals:
Min 1Q Median 3Q Max
-9.7171 -2.7586 -0.3636 2.7933 11.8628
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 75.1453 0.6741 111.48 <2e-16
estado_nutricionalmal_alimentado -15.4573 0.9439 -16.38 <2e-16
(Intercept) ***
estado_nutricionalmal_alimentado ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.718 on 98 degrees of freedom
Multiple R-squared: 0.7324, Adjusted R-squared: 0.7297
F-statistic: 268.2 on 1 and 98 DF, p-value: < 2.2e-16Nesse caso, o modelo retorna um coeficiente que nós chamamos de um contraste, mostrando a diferença na média entre os indivíduos bem alimentados e mal alimentados. Os indivíduos mal alimentados tendem a ser aproximadamente 15 kg mais leves que os bem alimentados.
# install.packages( "ggplot2") # se necessario...
library(ggplot2)
ggplot(dados, aes(estado_nutricional, peso, color = estado_nutricional)) + geom_boxplot() + geom_jitter(width = 0.1) 
O coeficiente beta representa a diferença na média desses dois conjuntos de observações.
Pergunta: Qual seria a estrutura dos resíduos dessa regressão?
Regressão múltipla
Nós podemos estender a regressão linear para mais de uma variável preditora. Nesse caso, a média esperada de Y é uma função linear de mais de uma variável X.
\[y_i \sim Normal(\mu_i, \sigma) \\ \mu_i = \alpha + \sum \beta_j x_{ij}\]
Relação entre regressão múltipla e simples
A relação entre os coeficientes da regressão múltipla e da regressão simples depende da relação entre as variáveis X. Quando as variáveis X são não correlacionadas, os coeficientes da regressão múltipla e da regressão simples são idênticos.
Vamos fazer um exemplo simulado para entender isso. Primeiro, vamos simular duas variáveis preditoras não correlacionadas, chamadas x e z.
Primeiro, vamos olhar para os coeficientes das regressões simples, de y em x, depois de y em z.
lm(y ~ x) |> coefficients() |> round(2)(Intercept) x
0.90 1.37 lm(y ~ z) |> coefficients() |> round(2)(Intercept) z
1.17 3.45 Os valores dos coeficientes de X e Z são muito parecidos com os valores simulados.
Finalmente, vamos olhar para os coeficientes da regressão múltipla, ajustando coeficientes para X e Y simultaneamente.
lm(y ~ x + z) |> coefficients() |> round(2)(Intercept) x z
1.03 1.37 3.45 Note como os valores dos coeficientes calculados pela regressão simples e pela regressão múltipla são idênticos. Isso acontece sempre que as variáveis preditoras são não correlacionadas.
Preditores correlacionados na regressão múltipla
A situação muda quando os preditores são correlacionados. Vamos simular um exemplo em que X e Z são ambos causados por uma mesma variável não observada U.
Agora, vamos novamente olhar para os coeficientes da regressão simples de Y em X e de Y em Z.
lm(y ~ x) |> coefficients() |> round(2)(Intercept) x
8.76 4.65 lm(y ~ z) |> coefficients() |> round(2)(Intercept) z
1.88 3.76 Veja como agora os coeficientes da regressão simples são completamente diferentes dos coeficientes que usávamos para fazer a simulação (\(b_x = 1.35\); \(b_z = 3.5\)).
Mas, se incluirmos tanto o X quanto o Z na regressão, conseguimos recuperar esses valores simulados.
lm(y ~ x + z) |> coefficients() |> round(2)(Intercept) x z
1.27 1.35 3.51 Recuperando os coeficientes da regressão múltipla usando residuos
Lembra-se de quando definimos os resíduos como a variação residual não explicada e vimos que poderíamos usá-los como uma nova variável independente?
Podemos reproduzir os coeficientes da regressão múltipla usando apenas regressões simples, aplicando o conceito de resíduos. Para descobrir o efeito isolado de \(x\) sobre \(y\) (descontando \(z\)), precisamos primeiro descontar o efeito de \(z\) na variável \(x\).
Podemos fazer isso em três passos:
- Ajustamos um modelo em que \(x\) é a variável resposta e \(z\) é a preditora.
- Extraímos os resíduos dessa regressão. Essa nova variável representará a variação de \(x\) que é independente de \(z\).
- Ajustamos um novo modelo simples de \(y\) contra esses resíduos de \(x\).
Vamos ver isso em código:
# 1. Regredimos a variável X contra a variável Z
m_xz <- lm(x ~ z)
# 2. Extraímos os resíduos (a variação de X que não tem relação com Z)
res_x <- residuals(m_xz)
# 3. Regredimos Y contra esse X residual
lm(y ~ res_x) |> coefficients() |> round(2)(Intercept) res_x
13.47 1.35 Observe a estimativa do coeficiente para res_x. Ela é idêntica ao coeficiente de \(x\) que obtivemos na regressão múltipla lm(y ~ x + z)!
Podemos fazer exatamente o mesmo procedimento para isolar o efeito de \(z\) sem a influência de \(x\):
# Isolando Z: regredimos Z em função de X e extraímos os resíduos
res_z <- residuals(lm(z ~ x))
# Regredimos Y contra o Z residual
lm(y ~ res_z) |> coefficients() |> round(2)(Intercept) res_z
13.47 3.51 Novamente, o coeficiente obtido na regressão simples bate com o coeficiente de \(z\) no modelo múltiplo.
Esse processo ilustra o Teorema de Frisch-Waugh-Lovell: a regressão múltipla não é nada além de uma série de regressões simples nas quais o efeito de cada variável preditora é estimado usando apenas a parte da sua variação que é única e não compartilhada com os demais preditores do modelo.
Controlando efeitos espúrios
Essas propriedades da regressão múltipla são usadas no contexto de genética quantitativa para remover fontes de variação que não são relevantes para o estudo em questão.
Vamos usar um exemplo com dados reais: o conjunto de dados ratones, que contém medidas cranianas de 5 populações experimentais de camundongos (https://onlinelibrary.wiley.com/doi/full/10.1111/evo.13304).
ratones <- read.table("https://raw.githubusercontent.com/diogro/evofencom/refs/heads/main/Tutoriais/ratones.tsv", header = TRUE)Esse é um conjunto de dados relativamente grande.
As primeiras colunas representam variáveis categóricas, como:
- ID de cada indivíduo
- O sexo de cada indivíduo
- A qual geração o indivíduo pertence
- A qual linhagem ele pertence
- O ID da mãe e do pai do indivíduo
- A idade do indivíduo em dias
- Qual linhagem ele pertence
- Qual o peso dele aos 49 dias
E a partir da coluna IS_PM até a coluna BA_OPI, temos as distâncias cranianas em milímetros.
str(ratones)'data.frame': 329 obs. of 46 variables:
$ .id : chr "control.t" "control.t" "control.t" "control.t" ...
$ ID : int 270158 270162 270164 270166 270189 270191 270193 270197 270199 270264 ...
$ SEX : chr "F" "F" "F" "F" ...
$ GER : int 55 55 55 55 55 55 55 55 55 55 ...
$ LIN : chr "t" "t" "t" "t" ...
$ MADRE : int 260902 260902 260918 260918 260902 260902 260918 260918 260918 260882 ...
$ PADRE : int 261089 261089 261065 261065 261089 261089 261065 261065 261065 261071 ...
$ AGE : num 170 118 117 169 113 113 154 112 112 168 ...
$ P49 : num 14.7 17.4 16.7 17 24.9 19.8 18.9 23.5 20.7 35.7 ...
$ line : chr "t" "t" "t" "t" ...
$ selection: chr "control" "control" "control" "control" ...
$ IS_PM : num 4.2 4.71 4.81 4.82 4.66 ...
$ IS_NSL : num 2.99 2.8 2.67 3.11 2.81 ...
$ IS_PNS : num 11.6 11.8 12.2 12.6 11.9 ...
$ PM_ZS : num 7.97 7.94 8.09 8.3 8.01 ...
$ PM_ZI : num 9.36 8.77 9.7 9.7 9.12 ...
$ PM_MT : num 6.99 6.91 7.14 7.21 7.03 ...
$ NSL_NA : num 8.49 9.07 9.23 9.31 9.4 ...
$ NSL_ZS : num 12 12.2 12.5 12.8 12.5 ...
$ NSL_ZI : num 14.1 14 14.8 15 14.5 ...
$ NA_BR : num 6.92 6.68 6.92 7.05 6.48 ...
$ NA_PNS : num 8.07 7.87 8.26 8.7 8.23 ...
$ BR_PT : num 4.47 3.47 4.11 3.69 5.16 ...
$ BR_APET : num 8.28 7.7 8.09 8.38 8.08 ...
$ PT_APET : num 8.04 8.29 8.54 8.84 8.3 ...
$ PT_BA : num 9.77 10.67 11.25 11.34 10.45 ...
$ PT_EAM : num 7.47 8.11 8.42 8.39 8.05 ...
$ PT_ZYGO : num 5.45 6.31 6.05 6.27 5.82 ...
$ PT_TSP : num 2.85 3.2 2.87 3.28 2.95 ...
$ ZS_ZI : num 2.52 2.59 2.66 2.84 2.72 ...
$ ZI_MT : num 4.81 4.49 5 5.06 4.5 ...
$ ZI_ZYGO : num 3 2.16 2.28 2.46 2.66 ...
$ ZI_TSP : num 4.08 4.22 3.74 4.06 3.99 ...
$ MT_PNS : num 1.85 1.78 1.9 1.92 1.95 ...
$ PNS_APET : num 3.43 4.72 4.49 4.61 4.68 ...
$ APET_BA : num 2.79 3.57 3.92 4.1 3.66 ...
$ APET_TS : num 3.83 3.66 3.9 4.05 3.97 ...
$ BA_EAM : num 5.51 5.42 5.69 5.76 5.57 ...
$ EAM_ZYGO : num 3.39 5.27 4.39 4.75 4.53 ...
$ ZYGO_TSP : num 4.18 4.18 3.89 4.04 4 ...
$ LD_AS : num 6.78 6.54 6.59 6.8 7.73 ...
$ BR_LD : num 7.89 8.27 7.95 7.02 9.65 ...
$ OPI_LD : num 3.64 3.68 3.64 4.92 4.17 ...
$ PT_AS : num 6.81 7.62 7.83 7.31 7.5 ...
$ JP_AS : num 4.76 4.14 4.45 4.35 4.14 ...
$ BA_OPI : num 4.48 3.82 4.12 4.05 3.78 ...Nosso objetivo é estimar a correlação fenotípica entre um par de variáveis cranianas: o comprimento do osso nasal e o comprimento do palato, dados pelas variáveis NSL_NA e IS_PNS.
Figure 1: Vista ventral e dorsal de um crânio de roedor, com os pontos envolvidos no cálculo das distâncias de interesse circulados em vermelho
Vamos começar onde devemos começar sempre, com uma visualização dos dados. Vamos fazer um gráfico da relação entre essas duas variáveis, usando cores e símbolos para visualizar os fatores de confusão: sexo e linhagem.
ggplot(ratones, aes(NSL_NA, IS_PNS, color = SEX, shape = line, group = line)) +
geom_point() + stat_ellipse() + scale_color_manual(values = 1:2) + theme_classic()
Parece haver uma relação entre as duas variáveis; porém, diferenças nas médias entre sexos e linhagens também são aparentes. As elipses estão desenhadas ao redor de cada linhagem, e podemos ver que a média das duas variáveis parece ser diferente em muitas das linhagens.
Antes de usar o modelo linear para controlar esses efeitos espúrios, vamos calcular a correlação bruta, sem nenhuma correção, entre as duas variáveis.
x = ratones$NSL_NA
y = ratones$IS_PNS
cor(x, y)[1] 0.3300104Agora, vamos usar um modelo linear para criar um conjunto de dados residual, removendo os efeitos de sexo e linhagem.
traits = names(ratones)[12:46]
f = paste0("cbind(", paste(traits, collapse = ","), ") ~ SEX + line")
# Replicando o formato original dos dados.
ratones_res = ratones
# Substituindo as medidas pelos resíduos de um modelo linear, incluindo sexo e linhagem.
ratones_res[, traits] = lm(as.formula(f), data = ratones) |> residuals()Esse novo conjunto de dados residual pode ser usado para fazer o mesmo gráfico que fizemos anteriormente, e calcular a mesma correlação.
ggplot(ratones_res, aes(NSL_NA, IS_PNS, color = SEX, shape = line, group = line)) +
geom_point() + stat_ellipse() + scale_color_manual(values = 1:2) + theme_classic()
Veja como agora as elipses estão todas sobrepostas: as diferenças nas médias foram removidas pelo modelo linear.
Podemos então calcular a correlação entre essas duas medidas, descontando os efeitos dessas variáveis de confusão.
x = ratones_res$NSL_NA
y = ratones_res$IS_PNS
cor(x, y)[1] 0.5596727A correlação, descontadas as variáveis de confusão, é sensivelmente maior do que o valor que havíamos calculado anteriormente. O efeito das variáveis espúrias (sexo e linhagem) estava reduzindo o valor da correlação bruta nos dados. Sem essa correção, nós estaríamos calculando uma correlação que não corresponde às relações fenotípicas encontradas dentro dessas populações.