Modelo Linear (Regressão Linear Múltipla) aplicada à fibropapilomatose no R

Tartarugas-marinhas são afligidas por tumores chamados fibropapilomatose (FP). Antigamente, alguns cientistas imaginavam que animais doentes (com FP) poderiam não conseguir se alimentar direito e portanto seriam mais magros do que animais saudáveis (sem FP). Hoje, muitos estudos ao redor do mundo indicam que animais com tumores são na verdade mais pesados, possivelmente devido ao peso que os tumores acrescentam na massa total da tartaruga. 

Imagine que você esteve no Projeto TAMAR 🐢 e conseguiu ter acesso aos dados de 29 tartarugas-verdes (Chelonia mydas) para essas variáveis: massa, comprimento curvilíneo da carapaça, largura curvilínea da carapaça, idade e presença/ausência de FP. Como um jovem ecólogo curioso que é, você faz o seguinte questionamento: 

• Pergunta:“Como o comprimento, largura, idade e FP podem influenciar a massa do animal?”

• Hipótese: Indivíduos com maior comprimento, maior largura, mais velhos e com FP terão mais massa."

Você então procura no Google alguma maneira de testar sua hipótese no R e lê algo sobre regressão linear múltipla (= modelo linear, na prática). Com ela você pode testar se múltiplas variáveis preditoras quantitativas (no caso, comprimento, largura e idade) ou categórica (no caso, presença/ausência de FP) podem influenciar uma variável resposta (massa). O modelo de regressão linear múltipla pode ser definido como

y = β0 ​+ β1​*x1 ​+ β2*​x2 ​+ ⋯ + βk*​xk

na qual x são as variáveis preditoras, y é a variável resposta, β0 é o ponto de intercepto da linha com o eixo y e βk​ é coeficiente angular da reta. Vamos agora rodar isso no R! 

1. Regressão Linear Múltipla 

Escreva o código abaixo no console do R: 

#Carrega os dados
> tartuguitalinear <- read.csv ("https://drive.google.com/u/0/uc?id=1J7dtTpKD2Ut3HKWTBKuyW0IKLDPKSFfC&export=download", sep=";")
> View (tartuguitalinear)

#Cria o objeto 'tortuguita' e usa a função 'lm' para modelar y ~ x1 + ... + xk
> tortuguita <- lm (massa ~ comprimento + largura + idade + fp, data = tartuguitalinear)

#Usa a função summary para apresentar os resultados
>summary (tortuguita)

Como resultado, temos então o seguinte diagnóstico:

Call:
lm(formula = massa ~ comprimento + largura + idade + fp, data = tartuguitalinear)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.226  -6.843   2.131   4.967  10.272 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -7.04322   18.45026  -0.382   0.7060  
comprimento  0.01685    1.28449   0.013   0.9896  
largura      0.52271    1.37346   0.381   0.7069  
idade        0.48036    0.26361   1.822   0.0809 .
fpyes        2.52356    3.79900   0.664   0.5128  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.878 on 24 degrees of freedom
Multiple R-squared:  0.7599,	Adjusted R-squared:  0.7199 
F-statistic: 18.99 on 4 and 24 DF,  p-value: 3.715e-07

Como interpretar esses resultados? Observe os Coefficients. Farei perguntas, tentem responder e depois olhem as respostas:

• Qual o significado de 0.01685, 0.52271 e 0.48036?

São os coeficientes angulares das variáveis preditoras comprimento (cm), largura (cm) e idade (anos de vida), respectivamente, em relação à variável resposta massa (kg). Por exemplo, podemos dizer que se idade, largura e FP forem constantes, espera-se que a cada aumento de um centímetro de comprimento (variável preditora) haja um aumento de 0.01685 kg na massa (variável resposta).

• Como interpretar o valor 2.52356 de fpyes?

No caso da variável categórica binária, fpyes significa que o modelo usou como referência (= 0) o estado de caráter no. Basicamente, ele codificou numericamente yes = 1 e no = 0. É assim que é possível misturar variáveis contínuas e categóricas em um mesmo modelo, sem precisar decorar aquela clássica tabela de tipos de testes de inferência estatística (teste t, ANOVA, Mann-Whitney, etc). Assim, se quisermos modelar para quando não há FP (fpno), devemos substituir x4 por zero na fórmula: 

y = β0 ​+ β1​*x1 ​+ β2*​x2 ​+ β3*​x3 + β4*​x4
massa = -7.04322 ​+ ​0.01685*comprimento ​+ 0.52271*​largura ​+ 0.48036*idade + 2.52356*0

Se quisermos modelar para quando há FP (fpyes), devemos substituir x4 por 1 na fórmula:

y = β0 ​+ β1​*x1 ​+ β2*​x2 ​+ β3*​x3 + β4*​x4
massa = -7.04322 ​+ ​0.01685*comprimento ​+ 0.52271*​largura ​+ 0.48036*idade + 2.52356*1

Com comprimento, largura e idade constantes, podemos concluir que a massa de tartarugas com FP em média é 2.52356 maior do que de tartarugas sem FP.

• Qual o significado do intercepto?

-7.04322 é a estimativa de massa quando o comprimento, largura e idade forem zero, e também não haver FP. Se quisermos estimar um indivíduo com comprimento, largura e idade = 0, mas com FP, isso resultaria em massa = -7.04322 + 1*2.52356. Obviamente não faz nenhum sentido uma tartaruga de massa negativa, o que demonstra a importância de se conhecer a biologia dos bichos e não só a matemática por trás dos modelos. 

• Como interpretar o R2?

Na regressão linear simples, geralmente interpretamos o Multiple R-squared como a porcentagem que a variação de x explica a variação de y. Contudo, a cada variável preditora adicionada na regressão linear múltipla, o R2 aumenta quando nem sempre deveria e passa a não ser confiável. Assim, devemos usar o Adjusted R-squared na regressão linear múltipla, pois ele penaliza o número de variáveis preditoras em x. O R2 ajustado apenas aumenta se a variável preditora de fato contribui para a explicação da variação em y.

• Como calcular o número de graus de liberdade? É útil saber?

Obtivemos 24 graus de liberdade, pois o nosso número amostral é 29 e as variáveis consideradas foram 5. Assim, 29 - 5 = 24. O número de graus de liberdade seria útil para calcular o intervalo de confiança para o coeficiente de cada variável preditora, mas o R já faz isso pra você.

• Vamos pensar em um teste de hipóteses: nossa hipótese nula (Ho) considera todos os coeficientes das variáveis preditoras como zero (nesse caso, a massa seria o intercepto), já nossa hipótese alternativa (Ha) considera pelo menos um coeficiente diferente de zero. O que o teste F nos informa para nos ajudar neste teste?

No teste F, note que p = 3.715e-07. Assim, o p < 0.05 do teste F nos diz que podemos considerar nosso modelo como um todo significativo. Contudo, não significa que o modelo realmente se ajusta 100% bem aos nossos dados, apenas nos informa que pelo menos um dos coeficientes das variáveis preditoras é diferente de 0, rejeitando a hipótese nula. Se p do teste F fosse maior que 0,05, significaria que a combinação das variáveis preditoras não foi boa para modelar os nossos dados.

Se quisermos fazer o teste de hipóteses para cada coeficientes das variáveis preditoras, basta ver se o p<0,05 no diagnóstico. Não houve valores de p menores que 0,05, então provavelmente não podemos confiar nos valores dos coeficientes de nenhuma das variáveis hauhauhauhuahahauhahua! 

2. Colinearidade e Parcimônia

Colinearidade deve ser evitada ao máximo! A colinearidade é quando uma variável preditora está altamente associada a outra variável preditora. Basicamente, uma delas não traz nada de informativo para o modelo (redundância) e só o torna desnecessariamente mais complexo. Por isso devemos adotar a filosofia da parcimônia, ou seja, o mínimo número de variáveis preditoras informativas para o modelo.

Como podemos saber se existe colinearidade? Basicamente podemos realizar uma seleção de modelo. Há 2 maneiras de fazer isso Backward (baseado em p ou R2 ajustado) ou Forward (baseado em p ou R2 ajustado também).

Backward

Nas seleções backward, começamos com o modelo com todas as variáveis preditoras e vamos dropando variáveis até chegarmos no modelo mais parcimonioso. Há muitos critérios sofisticados (como AIC e BIC), mas veremos critérios baseados em p e R2 ajustado:

(I) Baseado nos valores de p, devemos apagar a variável preditora com o maior p (> 0,05) e rodar o código de novo sem a variável apagada até que p < 0,05 para todas as preditoras. 

#Cria o objeto 'remove_comprimento', pois comprimento teve o maior p (0.9896) no diagnóstico
> remove_comprimento <- lm (massa ~  largura + idade + fp, data = tartuguitalinear)
> summary (remove_comprimento)

Como resultado, temos então o seguinte diagnóstico removendo comprimento das variáveis preditoras:

Call:
lm(formula = massa ~ largura + idade + fp, data = tartuguitalinear)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.232  -6.844   2.133   4.915  10.297 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -6.9140    15.2862  -0.452   0.6549  
largura       0.5396     0.4689   1.151   0.2607  
idade         0.4782     0.2001   2.390   0.0247 *
fpyes         2.5212     3.7181   0.678   0.5039  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.719 on 25 degrees of freedom
Multiple R-squared:  0.7599,	Adjusted R-squared:  0.7311 
F-statistic: 26.37 on 3 and 25 DF,  p-value: 6.518e-08

Observe que agora idade passou a ter um p significativo (0.0247). Significa que a colinearidade do comprimento com alguma outra variável preditora havia mudado nosso modelo para pior do que realmente é no diagnóstico original que fizemos no começo desse post. Poderíamos continuar removendo variáveis com valor de p > 0.05 até todos os valores de p se tornarem menores que 0.05, mas já deu para entender a lógica.

(II) Baseado nos valores de R2 ajustado, devemos repetir o código omitindo uma variável preditora de cada vez; no final escolhemos o modelo com maior R2. 

#Cria o objeto 'remove_comprimento' e usa a função 'lm' para modelar y ~ x1 + ... + xk
> remove_comprimento <- lm (massa ~  largura + idade + fp, data = tartuguitalinear)
> summary (remove_comprimento)
#Anote o R2 (=0.7311)

#Cria o objeto 'remove_largura' e usa a função 'lm' para modelar y ~ x1 + ... + xk
> remove_largura <- lm (massa ~  comprimento + idade + fp, data = tartuguitalinear)
> summary (remove_largura)
#Anote o R2 (=0.7295)

#Cria o objeto 'remove_idade' e usa a função 'lm' para modelar y ~ x1 + ... + xk 
> remove_idade <- lm (massa ~  comprimento + largura + fp, data = tartuguitalinear) 
> summary (remove_idade) 
#Anote o R2 (=0.6939 )

#Cria o objeto 'remove_fp' e usa a função 'lm' para modelar y ~ x1 + ... + xk
> remove_fp <- lm (massa ~  comprimento + largura + idade, data = tartuguitalinear) 
> summary (remove_fp) 
#Anote o R2 (=0.7261 )

Lembra que no modelo original com todas as variáveis preditoras obtivemos um R2 ajustado= 0.7199? Após começarmos a seleção de modelo backward, já encontramos um R2 maior que do modelo original. Um modelo sem comprimento parece mais parcimonioso comparado aos demais, então vamos explorar dropar comprimento + outra variável:

#Cria o objeto 'remove_comprimento.largura' e usa a função 'lm' para modelar y ~ x1 + ... + xk
> remove_comprimento.largura <- lm (massa ~  idade + fp, data = tartuguitalinear) 
> summary (remove_comprimento.largura) 
#Anote o R2 (= 0.7277)

#Cria  o objeto 'remove_comprimento.idade =sa a função 'lm' para modelar y ~ x1 + ... + xk
> remove_comprimento.idade <- lm (massa ~ largura + fp, data = tartuguitalinear)
> summary (remove_comprimento.idade)
#Anote o R2 (= 0.6824)

#Cria o objeto'remove_comprimento.fp" = a a função 'lm' para modelar y ~ x1 + ... + xk
> remove_comprimento.fp <- lm (massa ~ largura + idade, data = tartuguitalinear)
>  summary (remove_comprimento.fp)
#Anote o R2 (=0.7367 )

Encontramos um modelo ainda melhor removendo comprimento e FP. Em teoria, deveríamos remover o comprimento e FP (e talvez mais variáveis até que não encontremos um R2 mais alto), mas pouparei a nós todos para esse post não se alongar. Claro que se estou interessado em responder perguntas sobre a FP (mais do que em relação às outras variáveis), posso querer não excluir FP do modelo. Bom, já deu para entender a lógica.

PS.: Qual o melhor critério de eliminação backward, afinal? Depende... Backward por p permite criar modelos com variáveis preditoras mais significativas, enquanto backward por R2 permite criar modelos com variáveis preditoras mais explicativas. O mais comum é usar backward por p, mas é importante lembrar que ele utiliza níveis de significância arbitrários (0,05).

Forward

Nas seleções forward, começamos com um modelo com 1 variável preditora e vamos adicionando preditoras de interesse até alcançarmos o modelo mais parcimonioso (com maior R2 ou menor valor de p). Repetimos o processo até que não encontremos mais modelos com R2 maior ou com valores de p menores.

3. Premissas do Modelo de Regressão Linear Múltipla

• Relação linear entre variáveis preditoras e resposta: pode ser checado num gráfico de dispersão de resíduos, os quais devem estar aleatoriamente ao redor de 0. Abaixo, digamos que removendo comprimento e FP tenhamos encontrado nosso modelo mais parcimonioso. Checamos abaixo apenas para largura, mas devemos fazer para idade também.

> remove_comprimento.fp <- lm (massa ~ largura + idade, data = tartuguitalinear)
> plot(remove_comprimento.fp$residuals ~ tartuguitalinear$largura)

Não parece que o modelo passou nessa premissa, né?

• Resíduos devem ter distribuição normal ao redor de 0: pode ser checado usando plot de probabilidade normal e histograma dos resíduos (que deve estar ao redor de 0)

> hist (remove_comprimento.fp$residuals)
> qqnorm (remove_comprimento.fp$residuals)
> qqline (remove_comprimento.fp$residuals)

• Variabilidade constante dos resíduos: pode ser checado usando plot dos resíduos

> plot (remove_comprimento.fp$residuals ~ remove_comprimento$fitted)
> plot (abs(remove_comprimento.fp$residuals) ~ remove_comprimento$fitted)

• Independência: checado usando gráfico de dispersão de resíduos

> plot (remove_comprimento.fp$residuals)

Quebramos totalmente as premissas do modelo hahahahah bora testar com mais dados

LINKS
Além de Ecologia e Estatística, curte artes visuais? Dá uma olhada no meu Instagram!

E-MAIL

Mande suas críticas, elogios, sugestões e propostas para dani_ymn@hotmail.com