© Ricardo Solar/UFMG - compartilhamento e utilização não-comercial livres. Não reproduzir sem autorização >http://doi.org/10.5281/zenodo.7392285

Vimos na aula anterior perguntas que continham variáveis x categóricas, ou seja, havia classes predeterminadas para resposta. Nesses casos, utilizamos análises como o teste T e ANOVA. Agora, nosso interesse é avaliar uma variável x contínua, que pode assumir quaisquer valores. Nesse caso, utilizaremos um tipo de análise conhecido como regressão linear. Mais pra frente, iremos misturar variáveis x dos dois tipos, análise conhecida como ANCOVA.

Vamos utilizar a pergunta inicial de se o peso dos estudantes é função de sua altura, seguindo a hipótese de que indivíduos mais altos devem ter peso maior. Para essa primeira pergunta, utilizaremos uma Regressão Linear simples, pressupondo que há uma relação de causa e efeito entre essas variáveis.

1 Regressão linear

1.1 Importação dos dados

setwd("~/Dropbox/UFMG/Disciplinas/R UFAC/Aulas/teste t e anova")
dados <- read.table("~/Dropbox/UFMG/Disciplinas/R UFAC/Aulas/teste t e anova/RTurmas.txt", h = T, stringsAsFactors = T)

summary(dados)
##      ID_Ind          Altura           Peso        Sexo    Turma   
##  Min.   : 1.00   Min.   :150.0   Min.   : 44.50   F:45   Grad:45  
##  1st Qu.:23.75   1st Qu.:163.0   1st Qu.: 56.77   M:47   Pos :47  
##  Median :46.50   Median :168.0   Median : 66.02                   
##  Mean   :46.50   Mean   :170.1   Mean   : 68.33                   
##  3rd Qu.:69.25   3rd Qu.:177.2   3rd Qu.: 76.03                   
##  Max.   :92.00   Max.   :195.0   Max.   :113.70                   
##       Ano      
##  Min.   :2015  
##  1st Qu.:2016  
##  Median :2017  
##  Mean   :2017  
##  3rd Qu.:2017  
##  Max.   :2019

1.2 Criação e avaliação do modelo

Ainda trabalhamos com a premissa de que os resíduos seguem distribuição normal, conforme estamos trabalhando. De toda forma, os modelos precisarão ser avaliados através da análise de resíduos. Diferentemente do que utilizamos antes, aov para análise de variância, utilizaremos o comando lm (linear models). Na verdade, seguindo em frente, o comando lm abarca já o que é feito pelo aov, podemos seguir usando somente essa função.

modelo <- lm(Peso~Altura, data=dados)

Criado o modelo, o próximo passo é avaliar a adequabilidade do mesmo. Vamos utilizar aquela função simulateResiduals do pacote library(DHARMa). Para efeitos de treinamento, não irei carregar o pacote, mas apenas chamar uma função dentro dele, isso é feito utilizando a sintaxe pacote::NomeDaFunção, veja abaixo.

DHARMa::simulateResiduals(modelo, plot=T)
## qu = 0.25, log(sigma) = -3.112632 : outer Newton did not converge fully.

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.056 0.132 0.196 0.164 0.2 0.3 0.628 0.616 0.472 0.952 0.48 0.388 0.48 0.264 0.332 0.524 0.44 0.744 0.944 0.98 ...

Observe que a função até detectou uma variação, mas mostra que a análise combinada não foi significativa, ou seja, o modelo está OK. Posse seguir avaliando a significância do mesmo. Essa significância é avaliada pelo teste F, iremos disctuir em sala sobre como isso se processa do ponto de vista estatístico.

anova(modelo, test="F")

Ou seja, temos um efeito singificativo da altura do estudante no peso do mesmo. O próximo passo é visualizar isso num gráfico.

1.3 Construção do gráfico

No pacote ggplot2, a função geom_smooth é capaz de fazer essa curva sem nada além disso. Todavia, é importante que saibamos recuperar a equação da curva, e para isso, utilizaremos a função coef(modelo). Conforme já sabemos, a equação da reta é expressa como: \[ y \sim intercepto + inclinação \times x\]

Onde o intercepto é aquele ponto onde a reta toca o eixo Y, e a inclinação é quanto cada incremento e x reflete em y.

A função coef nos dá exatamente isso! Ela nos dará um vetor de dois elementos, o primeiro sendo o intercepto e o segundo a inclinação

coef(modelo)
## (Intercept)      Altura 
## -135.779875    1.199711

Portanto, nossa equação desta reta fica:

\[ y \sim -113.77 + 1.19 \times x\]

E na sintaxe de R, posso criar essa função utilizando os conteúdos de coef(modelo), sem ter que copiar esses valores. Vamos criar essa função, e a sintaxe em si, explicarei em sala:

reta.grafico <- function(x){coef(modelo)[1]+coef(modelo)[2]*x}

E agora irei plotar o gráfico em ggplot2, mas utilizando as duas formas de plotar a reta, vamos lá:

library(ggplot2)

ggplot(dados, aes(x=Altura, y=Peso))+
  geom_point()+
  geom_smooth(method = "lm", se=F)+ ##veja que apenas preciso especificar que tenho um lm. Esse é o método mais simples de plotar a curva, que será feita de maneira automática
  stat_function(fun=reta.grafico, col=2, linetype=2, linewidth=2) ##nesse caso, eu informo a equação da reta, que fiz acima. para que possamos ver ambas no mesmo gráfico, colocarei com a cor 2 (vermelho) e tipo de linha 2 (pontilhado) e com espessura grossa para visualização

O gráfico possui duas curvas idênticas, na prática, só precisamos de uma única! Para evitar erros ou confusões, eu prefiro utilizar sempre a equação informada por mim.

2 ANCOVA

Agora, quero refinar mais o meu modelo. Nossa hipótese agora é que o peso do estudante varia em função da altura e do sexo biológico de cada um. Nosso modelo fica um pouco mais elaborado. O que eu espero é que indivíduos do sexo masculino tenham peso maior do que indivíduos do sexo feminino e também que haja um efeito positivo de altura em peso. Não espero, todavia, que haja diferença na tendência da relação (não espero interação estatística). Sendo assim, meu modelo fica:

\[ Peso \sim Sexo + Altura \] ## Construção e avaliação do modelo

modelo2 <- lm(Peso~Sexo+Altura, data=dados)

Verificando o modelo - análise de resíduos (ver interpretação no tópico anterior)

DHARMa::simulateResiduals(modelo2, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.064 0.144 0.212 0.196 0.228 0.332 0.644 0.632 0.516 0.956 0.424 0.344 0.436 0.24 0.308 0.488 0.4 0.712 0.94 0.972 ...

E agora podemos avaliar a significância das variáveis.

anova(modelo2, test="F")

Ambas as variáveis são significativas. Como fica o plot desse gráfico? Como tenho duas categorias (M e F) e uma variável contínua, preciso construir duas retas. Não há interação, então isso singifica que temos apenas variação do intercepto, a inclinação é a mesma entre as retas. Teríamos então: \[ y_{fem} \sim intercepto_{fem} + inclinação \times x \]

\[ y_{masc} \sim intercepto_{masc} + inclinação \times x \]

Vamos ver os coeficientes:

coef(modelo2)
## (Intercept)       SexoM      Altura 
## -116.590917    3.202501    1.077305

O R mostra esses coeficientes de maneira peculiar. Ele mostra uma sequência Intercept, SexoM e Altura (inclinação). O que nos é dado aqui é o intercepto para a primeira categoria (em ordem alfabética) e depois a diferença na inclinação para a segunda categoria, que deve ser somada ou subtraída do intercepto da categoria 1.

Portanto, as duas equações ficam assim:

\[ y_{fem} \sim -116.59 + 1.07 \times x \]

\[ y_{masc} \sim (-116.59 + 3.20) + 1.07 \times x \] Entendido, isso, podemos partir para fazer as equações em forma de função do R, para podermos plotar isso no gráfico.

2.1 Confecção do gráfico

reta.fem <- function(x){coef(modelo2)[1]+coef(modelo2)[3]*x}
reta.mas <- function(x){(coef(modelo2)[1]+coef(modelo2)[2])+coef(modelo2)[3]*x}
ggplot(dados, aes(x=Altura, y=Peso, colour=Sexo))+
  geom_point(alpha=.8)+
  stat_function(fun=reta.fem, col="#6000CC", xlim = c(150, 180))+
  stat_function(fun=reta.mas, col="#FF8000", , xlim = c(160, 195))+
  scale_color_manual(values = c("#6000CC","#FF8000"))

Esse é o princípio básico das análises de regressão em ANCOVA. Iremos explorar isso bem mais nas próximas aulas.