Importa a base de dados crabs.csv
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/crabs.csv"
dados <- read.table(url, header = TRUE, sep = ";", dec = ",")
str(dados)
## 'data.frame': 156 obs. of 7 variables:
## $ especie: Factor w/ 2 levels "azul","laranja": 1 1 1 1 1 1 1 1 1 1 ...
## $ sexo : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
## $ FL : num 8.1 8.8 9.2 9.6 10.8 11.6 11.8 12.3 12.6 12.8 ...
## $ RW : num 6.7 7.7 7.8 7.9 9 9.1 10.5 11 10 10.9 ...
## $ CL : num 16.1 18.1 19 20.1 23 24.5 25.2 26.8 27.7 27.4 ...
## $ CW : num 19 20.8 22.4 23.1 26.5 28.4 29.3 31.5 31.7 31.5 ...
## $ BD : num 7 7.4 7.7 8.2 9.8 10.4 10.3 11.4 11.4 11 ...
Vamos analisar a relação que existe entre CL e CW
plot(CW ~ CL, data = dados)
Um modelo linear entre duas variáveis \(X\) e \(Y\), é definido matematicamente como uma equação com dois parâmetros desconhecidos,
\[ Y = \beta_0 + \beta_1 X \]
A análise de regressão é a técnica estatística que analisa as relações existentes entre uma única variável dependente, e uma ou mais variáveis independentes O objetivo é estudar as relações entre as variáveis, a partir de um modelo matemático, permitindo estimar o valor de uma variável a partir da outra
O problema da análise de regressão consiste em definir a forma de relação existente entre as variáveis. Por exemplo, podemos ter as seguintes relações
\[ \begin{align} Y &= \beta_0 + \beta_1 X &\qquad \text{linear} \\ Y &= \beta_0 X^{\beta_1} &\qquad \text{potência} \\ Y &= \beta_0 e^{\beta_1 X} &\qquad \text{exponencial} \\ Y &= \beta_0 + \beta_1 X + \beta_2 X^2 &\qquad \text{polinomial} \\ \end{align} \]
Em todos os casos, a variável dependente é \(Y\), aquela que será predita a partir da relação e da variável independente \(X\).
Em uma análise de regressão linear consideraremos apenas as variáveis que possuem uma relação linear entre si. Uma análise de regressão linear múltipla pode associar \(k\) variáveis independentes (\(X\)) para “explicar” uma única variável dependente (\(Y\)),
\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k + e \]
Uma análise de regressão linear simples associa uma única variável independente (\(X\)) com uma variável dependente (\(Y\)),
\[ Y = \beta_0 + \beta_1 X + e \]
Assim, dados \(n\) pares de valores, \((X_1, Y_1), (X_2, Y_2), \ldots, (X_n, Y_n)\), se for admitido que \(Y\) é função linear de \(X\), pode-se estabelecer uma regressão linear simples, cujo modelo estatístico é
\[ Y_i = \beta_0 + \beta_1 X_i + e_i, \quad i = 1, 2, \ldots, n \]
onde:
Interpretação dos parâmetros:
O problema agora consiste em estimar os parâmetros \(\beta_0\) e \(\beta_1\).
Como através de uma amostra obtemos uma estimativa da verdadeira equação de regressão, denominamos
\[ \hat{Y}_i = \hat{\beta}_0 + \hat{\beta}_1 X_i \]
ou seja, \(\hat{Y}_i\) é o valor estimado de \(Y_i\), através das estimativas de \(\beta_0\) e \(\beta_1\), que chamaremos de \(\hat{\beta}_0\) e \(\hat{\beta}_1\). Para cada valor de \(Y_i\), temos um valor \(\hat{Y}_i\) estimado pela equação de regressão,
\[ Y_i = \hat{Y}_i + e_i \]
Portanto, o erro (ou desvio) de cada observação em relação ao modelo adotado será
\[ \begin{align} e_i &= Y_i - \hat{Y}_i \\ e_i &= Y_i - (\beta_0 + \beta_1 X_i) \end{align} \]
Devemos então adotar um modelo cujos parâmetros \(\beta_0\) e \(\beta_1\), tornem esse diferença a menor possível. Isso equivale a minimizar a soma de quadrados dos resíduos (\(SQR\)), ou do erro,
\[ SQR = \sum_{i=1}^{n} [Y_i - (\beta_0 + \beta_1 X_i)]^2 \]
O método de minimizar a soma de quadrados dos resíduos é denominado de método dos mínimos quadrados. Para se encontrar o ponto mínimo de uma função, temos que obter as derivadas parciais em relação a cada parâmetro,
\[ \begin{align} \frac{\partial SQR}{\partial \beta_0} &= 2 \sum_{i=1}^{n} [Y_i - \beta_0 - \beta_1 X_i] (-1) \\ \frac{\partial SQR}{\partial \beta_1} &= 2 \sum_{i=1}^{n} [Y_i - \beta_0 - \beta_1 X_i] (-X_i) \end{align} \]
e igualar os resultados a zero
\[ \hat{\beta}_0 = \frac{\partial SQR}{\partial \beta_0} = 0 \qquad \text{e} \qquad \hat{\beta}_1 = \frac{\partial SQR}{\partial \beta_1} = 0 \]
Dessa forma, chegamos às estimativas de mínimos quadrados para os parâmetros \(\beta_0\) e \(\beta_1\):
\[ \begin{align} \hat{\beta}_1 &= \frac{\sum_{i=1}^{n} X_iY_i - \frac{\sum_{i=1}^{n} X_i \sum_{i=1}^{n} Y_i}{n}}{\sum_{i=1}^{n}X_i^2 - \frac{(\sum_{i=1}^{n} X_i)^2}{n}} \\ & \\ \hat{\beta_0} &= \bar{Y} - \hat{\beta}_1 \bar{X} \end{align} \]
onde
\[ \begin{align} \bar{Y} = \frac{1}{n} \sum_{i=1}^{n} Y_i \qquad \text{e} \qquad \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i \end{align} \]
Portanto, podemos usar estas soluções para calcular os parâmetros \(\hat{\beta_0}\) e \(\hat{\beta_1}\) na relação linear entre CL e CW, ou seja,
\[ CW = \hat{\beta}_0 + \hat{\beta}_1 CL + e \]
Como vimos pelas soluções acima, primeiro calculamos \(\hat{\beta_1}\), e depois \(\hat{\beta_0}\). Para facilitar as contas, vamos criar objetos X
e Y
com as colunas CL e CW respectivamente, e n
que é o tamanho da amostra
X <- dados$CL
Y <- dados$CW
n <- length(X)
Agora calculamos \(\hat{\beta_1}\) com
beta1 <- (sum(X * Y) - (sum(X) * sum(Y))/n)/(sum(X^2) - (sum(X)^2/n))
beta1
## [1] 1.097451
E \(\hat{\beta_0}\) é calculado com
beta0 <- mean(Y) - beta1 * mean(X)
beta0
## [1] 1.18695
Uma parte importante em uma análise de regressão linear é a verificação dos resíduos do modelo, ou seja, os desvios de cada valor observado \(Y\) em relação aos valores preditos pelo modelo, \(\hat{Y}\).
Os valores preditos (ou estimados) de \(Y\) pelo modelo podem ser calculados usando os coeficientes estimados
Y.pred <- beta0 + beta1 * X
Portanto, os resíduos podem ser calculados como a diferença entre cada valor observado, e cada valor predito de \(Y\), ou seja,
res <- Y - Y.pred
Como vimos que a suposição do modelo é de que os resíduos possuam uma distribuição normal com média 0 e variância constante, \(e_i \sim \text{N}(0, \sigma^2)\), então podemos verificar essa suposição fazendo um histograma destes resíduos
hist(res)
Note que por esta inspeção visual, a média dos resíduos parece estar próxima de 0, e todos os valores estão entre -2 e 2, o que indica que não há valores extremos em relaçao à normal padrão.
Observação: Existem outros métodos para analisar os resíduos e verificar as suposições do modelo, mas essa não é a intenção aqui.
Para ajustar um modelo linear no R, usamos a função lm()
mod <- lm(CW ~ CL, data = dados)
mod
##
## Call:
## lm(formula = CW ~ CL, data = dados)
##
## Coefficients:
## (Intercept) CL
## 1.187 1.097
Entre os diversos usos de fórmulas, o mais importante deles é sem dúvida o fato que fórmulas são utilizadas na declaração de modelos estatísticos. Um aspecto particularmente importante da linguagem S, o portanto no programa R, é que adota-se uma abordagem unificada para modelagem, o que inclui a sintaxe para especificação de modelos. Variáveis respostas e covariáveis (variáveis explanatórias) são sempre especificadas de mesma forma básica, ou seja, na forma
resposta ~ covariavel
onde:
~
significa “é modelada por”No restante deste texto vamos, por simplicidade, considerar que há apenas uma variável resposta que poderá ser explicada por uma ou mais covariáveis. Considere, o conjunto de dados do exemplo acima. Na sintaxe
lm(CW ~ CL, data = dados)
deve-se ler: CW é modelada por CL, através de um modelo linear lm()
, o que implica no modelo acima.
Vimos acima que o objeto mod
contém o modelo linear, e ao chamar esse objeto, vemos os coeficientes estimados, os mesmo que calculamos manualmente. Note a classe desse objeto:
class(mod)
## [1] "lm"
Existem uma série de funções genéricas capazes de extrair informações que estão nesse objeto. Objetos da classe lm
possuem muito mais informações do que os coeficientes. Você pode conferir quais são estas informações com
names(mod)
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
e para acessar algum destes componentes, fazemos como se fosse uma lista. Por exemplo, para acessar os coeficientes usamos
mod$coefficients
## (Intercept) CL
## 1.186950 1.097451
No entanto, também podemos usar a função genérica coef()
coef(mod)
## (Intercept) CL
## 1.186950 1.097451
Outras funções úteis nesse caso são:
fitted()
para acessar os valores ajustados, ou preditos pelo modeloresiduals()
para acessar os resíduos do modeloOutra função genérica extremamente importante para modelos de qualquer classe é a função summary()
summary(mod)
##
## Call:
## lm(formula = CW ~ CL, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7762 -0.5699 0.1098 0.4629 1.8273
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.186950 0.285340 4.16 5.28e-05 ***
## CL 1.097451 0.008698 126.17 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7827 on 154 degrees of freedom
## Multiple R-squared: 0.9904, Adjusted R-squared: 0.9904
## F-statistic: 1.592e+04 on 1 and 154 DF, p-value: < 2.2e-16
Note que o resultado é um resumo completo sobre a distribuição dos resíduos do modelo, e uma tabela com os coeficientes estimados, incluindo um teste \(t\) para a hipótese nula de que os coeficientes são iguais a zero.
Veja também que podemos atribuir o sumário a um outro objeto, verificar sua classe e seus componentes
mod.sum <- summary(mod)
class(mod.sum)
## [1] "summary.lm"
names(mod.sum)
## [1] "call" "terms" "residuals" "coefficients"
## [5] "aliased" "sigma" "df" "r.squared"
## [9] "adj.r.squared" "fstatistic" "cov.unscaled"
Portanto, verificamos que outras informações estão disponíveis também para esse objeto, resultante do mesmo modelo. Por exemplo, para acessar o valor de \(R^2\) do modelo, usamos
mod.sum$r.squared
## [1] 0.9904188
Outra vantagem de se ajustar um modelo ocm a função lm()
é que podemos também avaliar a significância geral do modelo através de uma Análise de Variância (ANOVA) para a regressão. Como vimos na estimação dos parâmetros, o objetivo é encontrar parâmetros que façam com que a soma de quadrados dos resíduos seja mínima. Podemos particionar a soma de quadrados da seguinte forma:
\[ SQTot = SQMod + SQRes \]
Portanto, se um modelo é bem ajustado, esperamos que a soma de quadrados do modelo \(SQMod\) seja grande, e a \(SQRes\) sejá mínima. Um quadro de ANOVA para o modelo irá testar, através de um teste F, se a soma de quadrados do modelo é significativamente diferente de zero. Para fazer essa ANOVA, usamos a função anova()
anova(mod)
## Analysis of Variance Table
##
## Response: CW
## Df Sum Sq Mean Sq F value Pr(>F)
## CL 1 9752.6 9752.6 15919 < 2.2e-16 ***
## Residuals 154 94.3 0.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Com esse resultado, vemos que o modelo linear é significativo, ou seja, a variável explicativa CL é útil para explicar a relação com a variável resposta CW.
Podemos também associar a saída dessa ANOVA para um objeto, e verificar sua classe
mod.an <- anova(mod)
class(mod.an)
## [1] "anova" "data.frame"
Veja que este objeto possui duas classes: anova
e data.frame
. Fundamentamente, a estrutura do objeto é a de um data.frame
, portanto, podemos acessar essa tabela usando indexação conforme qualquer outro data frame. Veja os nomes das colunas e linhas desse objeto
names(mod.an)
## [1] "Df" "Sum Sq" "Mean Sq" "F value" "Pr(>F)"
row.names(mod.an)
## [1] "CL" "Residuals"
Por exemplo, veja que o Residual standard error: 0.7827
, como reportado no objeto mod.sum
é o estimador do desvio-padrão residual \(\hat{\sigma}_{e} = \sqrt{\frac{\text{SQRes}}{n-2}}\), ou seja,
sqrt(mod.an$Sum[2]/mod.an$Df[2])
## [1] 0.7827079
e que F-statistic: 1.592e+04
(~15920) é o mesmo valor de
mod.an$F[1]
## [1] 15919.11
que testa a mesma hipótese da ANOVA. De fato, o valor de \(t^2\) para \(\beta_1\) no sumário do modelo é
mod.sum$coef[2,3]^2
## [1] 15919.11
Ainda aproveitando o objeto mod
, podemos ajustar o modelo graficamente ao gráfico da relação entre CW e CL. Para isso, usamos uma função gráfica de baixo nível, abline()
, utilizada para inserir linhas em gráficos. Esta função aplicada à um objeto da classe lm
produz o seguinte resultado.
plot(CW ~ CL, data = dados)
abline(mod)
Apenas para ilustrar a interpretação dos parâmetros, vamos alterar a escala dos eixos X e Y para visualizar o intercepto do modelo
plot(CW ~ CL, data = dados, xlim = c(0,50), ylim = c(0,55))
abline(mod)
Assim como fizemos antes, vamos verificar a adequção do modelo através dos resíduos. Para isso, podemos extrair os resíduos diretamente do objeto mod
com a função residuals()
, e fazer o histograma diretamente
hist(residuals(mod))
Veja que os resíduos retornados pela função são exatamente os mesmos daqueles que calculamos anteriormente (objeto res
)
all.equal(res, residuals(mod), check.names = FALSE)
## [1] TRUE
Isto é obviamente consequência do fato de que os valores preditos pelos coeficientes que calculamos à mão (objeto Y.pred
) são os mesmos daqueles retornados pela função fitted()
all.equal(Y.pred, fitted(mod), check.names = FALSE)
## [1] TRUE
Uma forma mais completa de se analisar os resíduos de um modelo é através do método lm
para a função genérica plot()
. Ao utilizar essa função em um objeto da classe lm
, por padrão são gerados quatro gráficos com diferentes informações sobre os resíduos
## Deixa a janela gráfica com 2 linhas e 2 colunas
par(mfrow = c(2, 2))
plot(mod)
par(mfrow = c(1, 1))
Com as colunas BD e CL do objeto dados
: