Revisão de Modelos Lineares

1 Regressão linear simples

O modelo de regressão linear simples assume que \[Y_i=\beta_0 + \beta_1 x_i + \epsilon_i\] com \(\epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\).

Para cada valor de \(x\), temos uma variável aleatória \(Y\) que é normalmente distribuída com a média dependendo do valor de \(x\).

Simule um conjunto de dados de tamanho \(n=1000\) e supondo \(\beta_0 = 100, \beta_1 = 2, \sigma^2 = 4, x \sim N(0,1)\).

n       = 1000
beta0   = 10
beta1   = 2
sigma2  = 0.2
x       = rnorm(n, 0, 1)
e       = rnorm(n, 0, sqrt(sigma2))
y       = beta0+beta1*x + e

#Analisando os dados gerados graficamente
plot(x,y,lwd=2, bty="n")

par(mfrow=c(1,2))
hist(y, freq=F, main="Histograma", xlab="Y", ylab="densidade")
curve(dnorm(x, mean(y), sd(y)), add=T, lwd=2)
hist(e, freq=F, main="Histograma", xlab="e", ylab="densidade")
curve(dnorm(x, mean(e), sd(e)), add=T, lwd=2)

#Analisando os dados gerados numericamente
ks.test(y, pnorm, mean(y), sd(y))$p.value
## [1] 0.9406667
shapiro.test(y)$p.value
## [1] 0.256065
ks.test(e, pnorm, mean(e), sd(e))$p.value
## [1] 0.8023007
shapiro.test(e)$p.value
## [1] 0.8562613

Note o que acontece quando geramos as covariáveis de outra forma:

x2       = rgamma(n, 1, 1)
y2       = beta0+beta1*x2 + e

#Analisando os dados gerados graficamente
plot(x2,y2,lwd=2, bty="n")

par(mfrow=c(1,2))
hist(y2, freq=F, main="Histograma", xlab="Y", ylab="densidade")
curve(dnorm(x, mean(y2), sd(y2)), add=T, lwd=2)
hist(e, freq=F, main="Histograma", xlab="e", ylab="densidade")
curve(dnorm(x, mean(e), sd(e)), add=T, lwd=2)

#Analisando os dados gerados numericamente
ks.test(y2, pnorm, mean(y2), sd(y2))$p.value
## [1] 8.926193e-14
shapiro.test(y2)$p.value
## [1] 4.153272e-29
ks.test(e, pnorm, mean(e), sd(e))$p.value
## [1] 0.8023007
shapiro.test(e)$p.value
## [1] 0.8562613

1.1 Estimação por máxima verossimilhança

Assuma que foram observados \(\textbf{y}=(y_1,\ldots,y_n)\). Seguindo o modelo, a função de verossimilhança é dada por \[\begin{align} L(\beta_0,\beta_1,\sigma^2;\textbf{y}) & =\prod_{i=1}^{n} f(y_i | \beta_0,\beta_1,\sigma^2) \nonumber \\ & = \prod_{i=1}^{n} (2\pi\sigma^2)^{-1/2}\exp\left\{-\dfrac{1}{2\sigma^2}(y_i - \beta_0 - \beta_1x_i)^2 \right \} \nonumber \\ & = (2\pi\sigma^2)^{-n/2}\exp\left\{-\dfrac{1}{2\sigma^2}\sum_{i=1}^{n}(y_i - \beta_0 - \beta_1x_i)^2 \right \} \nonumber \end{align}\]

O estimador de máxima verossimilhança é dado pelos valores de \(\beta_0\), \(\beta_1\) e \(\sigma^2\) que maximizam a função de verossimilhança. De forma similar, são os valores que maximizam o logaritmo de \(L(\beta_0,\beta_1,\sigma^2;\textbf{y})\)

Logaritmo da função de verossimilhança: \[ \log L(\beta_0,\beta_1,\sigma^2;\textbf{y}) = -\dfrac{n}{2} \log (2\pi\sigma^2)-\dfrac{1}{2\sigma^2}\sum_{i=1}^{n}(y_i - \beta_0 - \beta_1x_i)^2 \]

logvero = function(beta0, beta1, sigma2){
  soma = 0
  for(i in 1:n){
    soma = soma - 0.5*log(2*pi*sigma2) - 1/(2*sigma2)*(y[i]-beta0-beta1*x[i])^2
  }
  return(soma)
}

logvero2 = function(beta0, beta1, sigma2){
  return(sum(dnorm(y, beta0+beta1*x, sqrt(sigma2), log=T)))
}

logvero(beta0, beta1, sigma2)
## [1] -602.4624
logvero2(beta0, beta1, sigma2)
## [1] -602.4624
#Supondo beta1 e sigma2 conhecidos e
#beta0 desconhecido
b0 = seq(0,20,l=1000)
plot(b0, logvero(b0, beta1, sigma2), type="l", lwd=2, xlab=expression(beta[0]), ylab="log verossimilhança")
abline(v=beta0, lty=3, col="red", lwd=2)

Podemos mostrar que os estimadores de máxima verossimilhança são dados por

\[\begin{eqnarray*} \widehat{\beta}_1 &=& \dfrac{\sum_{i=1}^{n}x_iy_i - n\bar{x}\bar{y}}{\sum_{i=1}^{n}x_i^2 - n\bar{x}^2}\\ \widehat{\beta}_0 &=& \bar{y} - \widehat{\beta}_1 \bar{x}\\ \widehat{\sigma}^2 &=& \dfrac{1}{n}\sum_{i=1}^{n}(y_i - \widehat{\beta}_0 - \widehat{\beta}_1x_i)^2 \end{eqnarray*}\]

Lembre que \[ \sum_{i=1}^{n}\epsilon_i^2 = \mbox{soma dos quadrados dos resíduos (erros) (SQE)} \] Então, um outro estimador pontual de \(\sigma^2\) pode ser descrito pelo { erro quadrático médio} dado por \[ S^2 = MQE = \dfrac{SQE}{n-2}. \] Este é um estimador não viesado para \(\sigma^2\).

#Suponha que os parâmetros beta0, beta1 e sigma2 sejam desconhecidos.
#Cálculo do EMV destes parâmetros:

beta1_chapeu    = (sum(x*y) - n*mean(x)*mean(y))/(sum(x^2)-n*mean(x)^2)
beta0_chapeu    = mean(y) - beta1_chapeu*mean(x)
sigma2_chapeu   = sum((y-beta0_chapeu-beta1_chapeu*x)^2)/n
s2              = sum((y-beta0_chapeu-beta1_chapeu*x)^2)/(n-2)
c(beta1_chapeu, beta1)
## [1] 1.996359 2.000000
c(beta0_chapeu, beta0)
## [1]  9.986715 10.000000
c(sigma2_chapeu, s2, sigma2)
## [1] 0.1951063 0.1954973 0.2000000
#Usando uma função do R para isso.
ajuste      = lm(y~x)
ajuste_resu = summary(ajuste)

c(beta1_chapeu, beta1, ajuste$coefficients[2])
##                          x 
## 1.996359 2.000000 1.996359
c(beta0_chapeu, beta0, ajuste$coefficients[1])
##                         (Intercept) 
##    9.986715   10.000000    9.986715
c(sigma2_chapeu, s2, sigma2, ajuste_resu$sigma^2)
## [1] 0.1951063 0.1954973 0.2000000 0.1954973

1.2 Estimação intervalar dos parâmetros

  • Sabemos que \(\widehat{\beta}_0 \sim N(\beta_0,\sigma_0^2), \quad \sigma_0 = \sigma \sqrt{\dfrac{1}{n} + \dfrac{\bar{x}}{\sum_{i=1}^{n} (x_i-\bar{x})^2}}\).

Seja \(S_0=S \sqrt{\dfrac{1}{n} + \dfrac{\bar{x}^2}{\sum_{i=1}^{n} (x_i-\bar{x})^2}}\). Então, temos que \(T_0 = \dfrac{\widehat{\beta}_0-\beta_0}{S_0} \sim t_{n-2}\) e o intervalo de \(100(1-\alpha)\%\) de confiança para \(\beta_0\) é dado por \[ \widehat{\beta}_0 \pm t_{\alpha/n,n-2}S_0. \]

  • Sabemos que \(\widehat{\beta}_1 \sim N(\beta_1,\sigma_1^2), \quad \sigma_1 = \dfrac{\sigma}{\sqrt{\sum_{i=1}^{n} (x_i-\bar{x})^2}}\).

Seja \(S_1= \dfrac{S}{\sqrt{\sum_{i=1}^{n} (x_i-\bar{x})^2}}\). Como \(\widehat{\beta}_0\), \(\widehat{\beta}_1\) e \(S\) são mutuamente independentes, a quantidade \(T_1 = \dfrac{\widehat{\beta}_1-\beta_1}{S_1} \sim t_{n-2}\) e o intervalo de \(100(1-\alpha)\%\) de confiança para \(\beta_1\) é dado por \[ \widehat{\beta}_1 \pm t_{\alpha/n,n-2}S_1. \]

Mostrando empiricamente as distribuições acima:

beta0_chapeu_amostra = beta1_chapeu_amostra = sigma2_chapeu_amostra = s2_amostra = NULL

for(j in 1:1500){
y3       = beta0+beta1*x + rnorm(n, 0, sqrt(sigma2))
beta1_chapeu_amostra[j]    = (sum(x*y3) - n*mean(x)*mean(y3))/(sum(x^2)-n*mean(x)^2)
beta0_chapeu_amostra[j]    = mean(y3) - beta1_chapeu_amostra[j]*mean(x)
sigma2_chapeu_amostra[j]   = sum((y3-beta0_chapeu_amostra[j]-beta1_chapeu_amostra[j]*x)^2)/n
s2_amostra[j]              = sum((y3-beta0_chapeu_amostra[j]-beta1_chapeu_amostra[j]*x)^2)/(n-2)

}

Sxx = sum(  (x-mean(x))^2   )

par(mfrow=c(2,2))

hist(beta1_chapeu_amostra, freq=F, main="Histograma", xlab=expression(beta[1]), ylab="densidade")
curve(dnorm(x, beta1, sqrt(sigma2/Sxx)), add=T, lwd=2)

S1 = sqrt(s2_amostra) / sqrt( Sxx )
hist((beta1_chapeu_amostra-beta1)/S1, freq=F, main="Histograma", xlab="beta1 transformado", ylab="densidade")
curve(dt(x, n-2), add=T, lwd=2)

hist(beta0_chapeu_amostra, freq=F, main="Histograma", xlab=expression(beta[0]), ylab="densidade")
aux = sqrt(sigma2) * sqrt(1/n + mean(x)^2/Sxx)
curve(dnorm(x, beta0, aux ), add=T, lwd=2)

S0 = sqrt(s2_amostra) * sqrt(1/n + mean(x)^2/Sxx )
hist((beta0_chapeu_amostra-beta0)/S0, freq=F, main="Histograma", xlab="beta0 transformado", ylab="densidade")
curve(dt(x, n-2), add=T, lwd=2)

Calculando os ICs dos estimadores pontuais com base na primeira amostra de Y:

alfa = 0.05
quantil_t = qt(alfa, n-2, lower.tail=F) 
S1 = sqrt(s2) / sqrt( Sxx )
S0 = sqrt(s2) * sqrt(1/n + mean(x)^2/Sxx )
c(beta0_chapeu - quantil_t*S0, beta0_chapeu + quantil_t*S0)
## [1]  9.963696 10.009735
c(beta1_chapeu - quantil_t*S1, beta1_chapeu + quantil_t*S1)
## [1] 1.973941 2.018778

1.3 Estimação do valor médio de \(Y\) em \(x_h\)

Como o valor ajustado também é combinação linear dos erros, segue que \[ \widehat{Y}(x_h) \sim N\left(\mu(x_h), \sigma \sqrt{\dfrac{1}{n} + \dfrac{(x_h-\bar{x})^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}} \right) \]

O intervalo de \(100(1-\alpha)\%\) de confiança para \(\mu(x_h)\) é dado por \[ \widehat{\mu}(x_h) \pm t_{\alpha/n,n-2}S \sqrt{\dfrac{1}{n} + \dfrac{(x_h-\bar{x})^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}} \] Mostrando esse resultado empiricamente:

y_chapeu_amostra = beta0_chapeu_amostra = beta1_chapeu_amostra = sigma2_chapeu_amostra = s2_amostra = NULL

for(j in 1:1500){
y3       = beta0+beta1*x + rnorm(n, 0, sqrt(sigma2))
beta1_chapeu_amostra[j]    = (sum(x*y3) - n*mean(x)*mean(y3))/(sum(x^2)-n*mean(x)^2)
beta0_chapeu_amostra[j]    = mean(y3) - beta1_chapeu_amostra[j]*mean(x)
s2_amostra[j]              = sum((y3-beta0_chapeu_amostra[j]-beta1_chapeu_amostra[j]*x)^2)/(n-2)
y_chapeu_amostra[j] = beta0_chapeu_amostra[j] + beta1_chapeu_amostra[j] * x[10]
}

par(mfrow=c(1,2))

hist(y_chapeu_amostra, freq=F, main=paste("Histograma para x=",round(x[10],2),sep=""), xlab=expression(mu), ylab="densidade")
mu = beta0 + beta1 * x[10]
u = (x[10]-mean(x))^2
curve(dnorm(x, mu, sqrt(sigma2* (1/n+u/Sxx))), add=T, lwd=2)

Smu = sqrt(s2_amostra* (1/n+u/Sxx))
hist((y_chapeu_amostra-mu)/Smu, freq=F, main="Histograma", xlab="mu transformado", ylab="densidade")
curve(dt(x, n-2), add=T, lwd=2)

1.4 Previsão de \(Y\) em \(x_h\)

Agora \[ Y(x_h) - \widehat{Y}(x_h) \sim N\left(0, \sigma \sqrt{1+\dfrac{1}{n} + \dfrac{(x_h-\bar{x})^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}} \right) \]

O intervalo predição de \(100(1-\alpha)\%\) para uma nova observação \(Y(x_h)\) é dado por \[ \widehat{y}(x_h) \pm t_{\alpha/n,n-2}S \sqrt{1+\dfrac{1}{n} + \dfrac{(x_h-\bar{x})^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}} \]

Mostrando esse resultado empiricamente:

xprev = 0.8

yprev = yprev_amostra = beta0_chapeu_amostra = beta1_chapeu_amostra = s2_amostra = NULL

for(j in 1:1500){
y3        = beta0+beta1*x + rnorm(n, 0, sqrt(sigma2))
yprev[j]     = beta0+beta1*xprev + rnorm(1, 0, sqrt(sigma2))
  
beta1_chapeu_amostra[j]    = (sum(x*y3) - n*mean(x)*mean(y3))/(sum(x^2)-n*mean(x)^2)
beta0_chapeu_amostra[j]    = mean(y3) - beta1_chapeu_amostra[j]*mean(x)
s2_amostra[j]              = sum((y3-beta0_chapeu_amostra[j]-beta1_chapeu_amostra[j]*x)^2)/(n-2)
yprev_amostra[j]           = beta0_chapeu_amostra[j] + beta1_chapeu_amostra[j] * xprev
}

par(mfrow=c(1,2))

hist(yprev-yprev_amostra, freq=F, main=paste("Histograma para x=",round(xprev,2),sep=""), xlab="y previsto transf 1", ylab="densidade")
mu = beta0 + beta1 * xprev
u = (xprev-mean(x))^2
aux = 1 + 1/n + u/Sxx
curve(dnorm(x, 0, sqrt(sigma2* aux)), add=T, lwd=2)

Smu = sqrt(s2_amostra* aux)
hist((yprev-yprev_amostra)/Smu, freq=F, main="Histograma", xlab="y previsto transf 2", ylab="densidade")
curve(dt(x, n-2), add=T, lwd=2)

1.5 Teste de hipótese sobre o coeficiente angular

Neste caso, o interesse é testar a seguintes hipóteses \[\begin{align*} & H_0: \beta_1 = 0 \\ & H_1: \beta_1 \neq 0 \end{align*}\]

A estatística \(T= \dfrac{\widehat{\beta}_1-\beta_1}{S_1} \sim t_{n-2}\) e, sob \(H_0\), \(T=\dfrac{\widehat{\beta}_1}{S_1} \sim t_{n-2}\).

Voltando a primeira amostra de Y:

ajuste_resu
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.58841 -0.29138 -0.01045  0.29827  1.38085 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.98672    0.01398   714.2   <2e-16 ***
## x            1.99636    0.01362   146.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4422 on 998 degrees of freedom
## Multiple R-squared:  0.9556, Adjusted R-squared:  0.9556 
## F-statistic: 2.149e+04 on 1 and 998 DF,  p-value: < 2.2e-16
ajuste2 = lm(y~rexp(n,3))
ajuste2_resu = summary(ajuste2)
ajuste2_resu
## 
## Call:
## lm(formula = y ~ rexp(n, 3))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9833 -1.4388 -0.0086  1.3974  5.8113 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.89980    0.09312 106.318   <2e-16 ***
## rexp(n, 3)   0.26513    0.18478   1.435    0.152    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.097 on 998 degrees of freedom
## Multiple R-squared:  0.002059,   Adjusted R-squared:  0.001059 
## F-statistic: 2.059 on 1 and 998 DF,  p-value: 0.1516

1.6 Coeficiente de determinação \(R^2\)

Suponha que não tivéssemos covariáveis. Nesse caso, nossa estimativa seria baseada apenas nos valores de \(Y\). E então o EMV de \(\mu\) seria \(\bar{y}\). Neste caso, os erros seriam dados por \(y_i-\bar{y}\) e a variação em torno da estimativa, medida pela variância amostral, seria proporcional a \[SQT = \sum_{i=1}^{n} (y_i - \bar{y})^2,\] denotada por { Soma de Quadrados Total}.

Mas se utilizarmos a informação sobre as variáveis explicativas, \(e_i = y_i - \widehat{y}_i\), a variação associada a estes erros pode ser medida por \[SQE = \sum_{i=1}^{n} (y_i - \widehat{y}_i)^2,\] que é a { Soma dos Quadrados dos Resíduos}.

A melhora das estimativas, como resultado do modelo linear, pode ser medida através da diferença \[ (y_i - \bar{y}) - (y_i - \widehat{y}_i) = (\widehat{y}_i - \bar{y}) \] e medimos a variação nesses erros a partir de \[ SQR = \sum_{i=1}^{n} (\widehat{y}_i-\bar{y})^2 \] conhecido como { Soma dos Quadrados da Regressão}

Pode-se mostrar que \[\begin{align*} \sum_{i=1}^{n} (y_i - \bar{y})^2 & = \sum_{i=1}^{n}(\widehat{y}_i - \bar{y}) ^2 + \sum_{i=1}^{n} (y_i - \widehat{y}_i)^2 \\ SQT & = SQR + SQE \end{align*}\]

Portanto, a variabilidade total é decomposta na variabilidade explicada pela regressão e a variabilidade residual (não ``explicável”)

Um índice bastante utilizado para resumir quão bom é o ajuste de um modelo é chamado de { coeficiente de determinação} e é dado por \[ R^2 = 1- \dfrac{SQE}{SQT} \]

y_chapeu = beta0_chapeu + beta1_chapeu*x
SQE = sum((y-y_chapeu)^2)
SQR = sum((y_chapeu-mean(y))^2)
SQT = SQR + SQE
R2 = 1-SQE/SQT
c(R2, ajuste_resu$r.squared)
## [1] 0.9556291 0.9556291
ajuste2_resu$r.squared
## [1] 0.002058739

1.7 Verificando bondade do ajuste

Equivale a realizar análise dos resíduos.

1º) Verifico Normalidade

  • qqplot para verificar a normalidade dos resíduos

  • Teste de Shapiro-Wilk

  • Teste de Kolmogorov-Smirnov

residuos = y-(beta0_chapeu + beta1_chapeu*x)
residuos2 = ajuste_resu$residuals
plot(residuos, residuos2)

ks.test(residuos, pnorm, mean(residuos), sd(residuos))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.019654, p-value = 0.8346
## alternative hypothesis: two-sided
shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.99896, p-value = 0.8524
par(mfrow=c(1,2))
hist(residuos, freq=F, main="Histograma", xlab="resíduos", ylab="densidade")
curve(dnorm(x, mean(residuos), sd(residuos)), add=T, lwd=2)

qqnorm(residuos,xlab="quantis teóricos", ylab="quantis amostrais", bty="n")
qqline(residuos,col="red", lwd=2)

2º) Verifico homocedasticidade

Sabe-se que \(V(e_i)=\sigma^2(1-h_{ii})\), \(i=1,\ldots,n\). Consequentemente, a variância dos resíduos não é constante. Portanto, é sugerido padronizar os resíduos, usando a estimativa do erro padrão de \(e_i\). Assim, o resíduo padrão como \[ r_i = \dfrac{e_i}{S\sqrt{(1-h_{ii})}}, i=1,\ldots,n \]

Quando avaliamos um scatter plot de \(\sqrt{|r_i|}\) \(\widehat{y}_i\), espera-se ver uma faixa constante.

#install.packages("lmtest")
library(lmtest)
bptest(ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 0.072171, df = 1, p-value = 0.7882
par(mfrow=c(1,2))
plot(residuos, ajuste$fitted.values)
plot(residuos, x)

3º) Verifico independência

Teste de Durbin-Watson, baseado na estatística \[D=\dfrac{\sum_{i=1}^{n}(e_i-e_{i-1})^2}{\sum_{j=1}^{n}e_j^2}\]

library(car)

# Realize o teste de Durbin-Watson
durbinWatsonTest(ajuste)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.03773921      1.924359   0.256
##  Alternative hypothesis: rho != 0

1.8 Alternativas quando o modelo não está adequado

  • A resposta média é não linear
  • A variância do erro NÃO é constante. Podemos utilizar alguma transformação para resolver o problema, do tipo \(Y^* = Y^{\lambda}\) (transformação de Box-Cox)
  • Se a variável resposta NÃO é normal, as mesmas transformações discutidas para estabilização da variância podem ser aplicadas

1.9 Outros diagnósticos

  • Observações influentes

Influenciam as estimativas, previsões e inferência. Uma pequena mudança em uma observação influente resulta numa grande mudança nas estimativas dos parâmetros ou inferência

  • Outliers

São aquelas observações que estão bem distantes dos dados. Podem indicar falta de ajuste do modelo de regressão

  • Leverages (Pontos de alavanca)

Os pontos de alavanca identificam observações que têm valores de \(x\) que estão distantes do restante dos dados. No modelo de regressão simples, a alavancagem de \(x_i\) é denotada por \(h_{ii}\) e definida por \[ h_{ii} = \dfrac{1}{n} + \dfrac{(x_i-\bar{x})^2}{\sum_{k=1}^{n}(x_i-\bar{x})^2} \] Se a distância de \(x_i\) e \(\bar{x}\) é grande, então \(h_{ii} \approx 1\). Uma sugestão é considerar valores de alavancagem altos que são duas vezes maiores do que o valor médio \(2/n\)

2 Regressão linear múltipla

O modelo de regressão linear múltipla usual é descrito por \[ Y_i = \sum_{j=0}^{p} x_{ij}\beta_j + \epsilon_i, \ \ i=1,\ldots,n; \ \ j=0,\ldots,p \]

  • \(y_i\) é a variável resposta registrada para a \(i\)-ésima observação
  • \(x_{ij}\) é a \(j\)-ésima variável explicativa, registrada para a \(i\)-ésima observação
  • \(\beta_j\) são os coeficientes de regressão
  • \(\epsilon_i\) é o termo de erro, independente e identicamente distribuído \(\forall i\). Em geral, \(\epsilon_i \sim N(0,\sigma^2)\)

Em forma matricial \[ \textbf{Y} = \textbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, \] com \(\boldsymbol{\epsilon} \sim NMV(\textbf{0},\sigma^2 \textbf{I}_n)\),

onde

\(\boldsymbol{Y} = \left[\begin{array}{c} y_1 \\ y_2 \\ \vdots\\ y_n \end{array} \right]\)
\(\boldsymbol{X} = \left[\begin{array}{c x x x x} 1 & x_{11} & x_{12} & ... & x_{1p} \\ 1 & x_{21} & x_{22} & ... & x_{2p} \\ \vdots & & & \ddots & \\ 1 & x_{n1} & x_{n2} & ... & x_{np} \\ \end{array} \right]\) \(\boldsymbol{\beta} = \left[\begin{array}{c} \beta_0 \\ \beta_1 \\ \vdots\\ \beta_p \end{array} \right]\) \(\boldsymbol{\epsilon} = \left[\begin{array}{c} \epsilon_1 \\ \epsilon_2 \\ \vdots\\ \epsilon_n \end{array} \right]\),

#gerando dados
n         = 1000
beta0     = 10
beta1     = 2
beta2     = -1
beta      = c(beta0, beta1, beta2)
sigma2    = 0.1
x0        = rep(1,n)
x1        = rnorm(n)
x2        = rgamma(n,1,1)
x         = cbind(x0,x1,x2)
e         = rnorm(n, 0, sqrt(sigma2))

y         = x%*%beta + e

#Analisando os dados gerados graficamente
par(mfrow=c(2,2))
plot(x1,y,lwd=2, bty="n")
plot(x2,y,lwd=2, bty="n")
hist(y, freq=F, main="Histograma", xlab="Y", ylab="densidade")
curve(dnorm(x, mean(y), sd(y)), add=T, lwd=2)
hist(e, freq=F, main="Histograma", xlab="e", ylab="densidade")
curve(dnorm(x, mean(e), sd(e)), add=T, lwd=2)

#Analisando os dados gerados numericamente
ks.test(y, pnorm, mean(y), sd(y))$p.value
## [1] 0.1305274
shapiro.test(y)$p.value
## [1] 0.0002856787
ks.test(e, pnorm, mean(e), sd(e))$p.value
## [1] 0.8013049
shapiro.test(e)$p.value
## [1] 0.9741092

2.1 Estimação por máxima verossimilhança

Assuma que foram observados \(\textbf{y} = (y_1, \ldots, y_n)\), seguindo o modelo, a função de verossimilhança é dada por

\[ L(\beta_0,\beta_1,\sigma^2;\textbf{y}) = (2\pi\sigma^2)^{-n/2}\exp\left\{-\dfrac{1}{2\sigma^2} (\textbf{y}-\textbf{X}\boldsymbol{\beta})'(\textbf{y}-\textbf{X}\boldsymbol{\beta}) \right \} \nonumber \] O EMV para \(\boldsymbol{\beta}\) é dado por \[ \widehat{\boldsymbol{\beta}}=(\textbf{X}'\textbf{X})^{-1}\textbf{X}'\textbf{y} \]

2.2 Estimativas pontuais da superfície de regressão

\[ \widehat{\textbf{y}} = \textbf{X}\widehat{\boldsymbol{\beta}} = \textbf{X}\left[(\textbf{X}'\textbf{X})^{-1}\textbf{X}'\textbf{y}\right] = \textbf{Hy} \] A matriz \(\textbf{H}\) é conhecida como matriz chapéu (hat matrix)

  • \(\textbf{H}\) é uma matriz quadrada e simétrica, de dimensão \(n\)

  • Os elementos da diagonal, \(h_{ii}\), satisfazem a \(0 \leq h_{ii} \leq 1\)

  • \(tr(\textbf{H})=p+1\)

  • \(\textbf{H}\) é uma matriz de projeção, de modo que \(\textbf{H}^2=\textbf{H}\), e o mesmo vale para \((\textbf{I}-\textbf{H})\)

2.3 Erro quadrático médio

\[ \textbf{e} = \textbf{y}-\widehat{\textbf{y}} = \textbf{y} - \textbf{Hy} = (\textbf{I}-\textbf{H})\textbf{y} \] - Pode ser mostrado que \(\textbf{e} \sim NMV(\textbf{0},\sigma^2(\textbf{I}-\textbf{H}))\)

  • \(SQE=\textbf{e}'\textbf{e} = \textbf{y}'(\textbf{I}-\textbf{H})(\textbf{I}-\textbf{H})\textbf{y}= \textbf{y}'(\textbf{I}-\textbf{H})\textbf{y}\)

  • \(\widehat{\sigma}^2 = \dfrac{SQE}{n-(p+1)}\)

  • O erro padrão residual é estimado por \(S=\sqrt{S^2}\)

2.4 Coeficiente de determinação múltiplo \(R^2\)

Pode ser mostrado que \(SQE=\textbf{y}'(\textbf{I}-\textbf{H})\textbf{y}\) e \(SQT=\textbf{y}'\left(\textbf{I}-\dfrac{1}{n}\textbf{J}\right)\textbf{y}\). Logo, \[ SQT = SQE + SQR \] e \[ R^2 = 1-\dfrac{SQE}{SQT} \]

  • O valor de \(R^2\) pode ser artificialmente inflado através da inclusão de outras covariáveis.

  • Uma alternativa é penalizar o valor de \(R^2\) pelo número de covariáveis, denominado \(R^2\) ajustado e calculado a partir de \[ R_a^2 = \left(R^2 - \dfrac{p}{n-1} \right) \left(\dfrac{n-1}{n-p-1} \right) \]

2.5 Teste F e ANOVA

Outra forma de verificar a bondade de ajuste é através de um teste F global, que testa as seguintes hipóteses \[\begin{align*} & H_0: \beta_1 = \beta_2 = \ldots = \beta_p = 0 \\ & H_1: \mbox{pelo menos um } \beta_j \neq 0 \end{align*}\] A estatística de teste é dada por \[ F=\dfrac{SQR/p}{SQE/(n-p-1)} \] Sib \(H_0\), \(F \sim F_{p,n-p-1}\). Rejeitamos \(H_0\) se F é muito grande. 
ajuste = lm(y~x-1)
ajuste_resu = summary(ajuste)
ajuste_resu
## 
## Call:
## lm(formula = y ~ x - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.00503 -0.20719 -0.00043  0.20355  1.07959 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## xx0  9.994919   0.014084   709.6   <2e-16 ***
## xx1  2.020491   0.010155   199.0   <2e-16 ***
## xx2 -0.996428   0.009824  -101.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3127 on 997 degrees of freedom
## Multiple R-squared:  0.9988, Adjusted R-squared:  0.9988 
## F-statistic: 2.864e+05 on 3 and 997 DF,  p-value: < 2.2e-16
residuos = ajuste_resu$residuals
par(mfrow=c(2,2))
plot(residuos)
plot(residuos,y)
plot(residuos,x1)
plot(residuos,x2)

ks.test(residuos, pnorm, mean(residuos), sd(residuos))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.019637, p-value = 0.8353
## alternative hypothesis: two-sided
shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.99922, p-value = 0.9597
par(mfrow=c(1,2))
hist(residuos, freq=F, main="Histograma", xlab="resíduos", ylab="densidade")
curve(dnorm(x, mean(residuos), sd(residuos)), add=T, lwd=2)

qqnorm(residuos,xlab="quantis teóricos", ylab="quantis amostrais", bty="n")
qqline(residuos,col="red", lwd=2)

bptest(ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste
## BP = 2.073, df = 2, p-value = 0.3547
durbinWatsonTest(ajuste)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.04997154       1.89887   0.132
##  Alternative hypothesis: rho != 0
#Usando uma única covariável para o ajuste:
ajuste2 = lm(y~x2)
ajuste2_resu = summary(ajuste2)
ajuste2_resu
## 
## Call:
## lm(formula = y ~ x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9786 -1.1453  0.0291  1.3568  6.5217 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.86257    0.08972  109.93   <2e-16 ***
## x2          -0.94179    0.06263  -15.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.994 on 998 degrees of freedom
## Multiple R-squared:  0.1847, Adjusted R-squared:  0.1839 
## F-statistic: 226.2 on 1 and 998 DF,  p-value: < 2.2e-16
residuos2 = ajuste2_resu$residuals
par(mfrow=c(2,2))
plot(residuos2)
plot(residuos2,y)
plot(residuos2,x1)
plot(residuos2,x2)

ks.test(residuos2, pnorm, mean(residuos2), sd(residuos2))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos2
## D = 0.034593, p-value = 0.1825
## alternative hypothesis: two-sided
shapiro.test(residuos2)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos2
## W = 0.99653, p-value = 0.02641
par(mfrow=c(1,2))
hist(residuos2, freq=F, main="Histograma", xlab="resíduos", ylab="densidade")
curve(dnorm(x, mean(residuos2), sd(residuos2)), add=T, lwd=2)

qqnorm(residuos2,xlab="quantis teóricos", ylab="quantis amostrais", bty="n")
qqline(residuos2,col="red", lwd=2)

bptest(ajuste2)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste2
## BP = 0.12432, df = 1, p-value = 0.7244
durbinWatsonTest(ajuste2)
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.03885523      2.074172   0.196
##  Alternative hypothesis: rho != 0
#usando uma covariável errada
x3 = rexp(n,9)
ajuste3 = lm(y~x3)
ajuste3_resu = summary(ajuste3)
ajuste3_resu
## 
## Call:
## lm(formula = y ~ x3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.1867  -1.3426   0.0718   1.5458   7.3692 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.8868     0.1022  86.917   <2e-16 ***
## x3            0.1453     0.6751   0.215     0.83    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.209 on 998 degrees of freedom
## Multiple R-squared:  4.644e-05,  Adjusted R-squared:  -0.0009555 
## F-statistic: 0.04635 on 1 and 998 DF,  p-value: 0.8296
residuos3 = ajuste3_resu$residuals
par(mfrow=c(3,3))
plot(residuos3)
plot(residuos3,y)
plot(residuos3,x3)

ks.test(residuos3, pnorm, mean(residuos3), sd(residuos3))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuos3
## D = 0.037784, p-value = 0.1151
## alternative hypothesis: two-sided
shapiro.test(residuos3)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos3
## W = 0.99359, p-value = 0.0002706
par(mfrow=c(1,3))

hist(residuos3, freq=F, main="Histograma", xlab="resíduos", ylab="densidade")
curve(dnorm(x, mean(residuos3), sd(residuos3)), add=T, lwd=2)

qqnorm(residuos3,xlab="quantis teóricos", ylab="quantis amostrais", bty="n")
qqline(residuos3,col="red", lwd=2)

bptest(ajuste3)
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste3
## BP = 1.1487, df = 1, p-value = 0.2838
durbinWatsonTest(ajuste3)
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.01379149      2.024095   0.672
##  Alternative hypothesis: rho != 0

2.6 Comparação de modelos

  • Critério de Informação de Akaike

    \[AIC = -2\log \widehat{L} + 2(p+1) \]

    Escolhe-se o modelo com o menor valor de \(AIC\)

  • Critério de Informação Bayesiano (de Schwarz) \[ BIC = -2\log \widehat{L} + (p+1) \log n \] O termo de penalização \(BIC\) é mais rigoroso do que o termos de penalização \(AIC\). Por exemplo, para \(n \geq 8\), \((p+1)\log n\) excede \(2(p+1)\). Aqui, também escolhemos o modelo com o menor valor de \(BIC\)

c(AIC(ajuste), AIC(ajuste2), AIC(ajuste3))
## [1]  517.9007 4222.2978 4426.5023
c(BIC(ajuste), BIC(ajuste2), BIC(ajuste3))
## [1]  537.5317 4237.0211 4441.2255