# Análise de Regressão
# Instale o R a partir de www.r-project.org
# Utilize um editor de R, preferencialmente o https://www.rstudio.com/
# Deseja-se explicar o consumo de combustível para 48 estados norte-americanos 
# em 1971, com as seguintes variáveis explicativas 
# X1 o imposto sobre combustível (centavos de dólar/galão), 
# X2 o percentual da população com carteira de habilitação, 
# X3 renda per capita (milhares de dólares) e 
# X4 a quantidade de estradas federais pavimentadas (em milhares de milhas).
# A variável resposta é o consumo de combustível Y (galões por pessoa).
# Um galão equivale a 3,785 litros. (Weisberg, S. 1985. Applied Linear Regression, Wiley).

library(MASS)

# Leitura dos dados
combustivel <- read.table("http://wiki.icmc.usp.br/images/b/bc/Combustivel.txt", header=T)
attach(combustivel)
names(combustivel)

# Definindo a matriz X 
X <- cbind(1,X1,X2,X3,X4)

# Ajustando o modelo linear com a suposição de erros independentes
beta.hat <- solve(t(X) %*% X) %*% t(X) %*% Y
beta.hat

# Forma Alternativa 1
fit.model <- lm(Y ~ X1+X2+X3+X4)
beta.hat <- fit.model$coefficients
beta.hat

# Forma Alternativa 2 - X já tem uma coluna de 1's
fit.modelA <- lm(Y ~ -1 + X)
beta.hatA <- fit.modelA$coefficients
beta.hatA

# Significância da Regressão
anova(fit.model)
anova(fit.modelA)

# Novo modelo - sem X4
fit.model2 <- lm(Y ~ X1+X2+X3)
beta.hat2 <- fit.model2$coefficients
beta.hat2

beta.hat

# Matriz do modelo
model.matrix(fit.model2)

# Significância da Regressão
anova(fit.model2)

# Teste de Hipóteses
anova(fit.model, fit.model2, test="F")

fit.model0<-lm(Y ~ -1 + X1 + X2 + X3 + X4)

anova(fit.model, fit.model0, test="F")

# Análise de resíduos
# Resíduos x valores ajustados
plot(fitted(fit.model), stdres(fit.model),pch=16, ylim=range(stdres(fit.model),-4,4))
abline(qnorm(0.975),0,lty=2,col=2)
abline(qnorm(0.025),0,lty=2,col=2)

plot(fitted(fit.model2),stdres(fit.model2),pch=16, ylim=range(stdres(fit.model2),-4,4))
abline(qnorm(0.975),0,lty=2,col=2)
abline(qnorm(0.025),0,lty=2,col=2)

# Resíduos x variáveis explicativas
plot(X1, stdres(fit.model2),pch=16, ylim=range(stdres(fit.model),-4,4))
abline(qnorm(0.975),0,lty=2,col=2)
abline(qnorm(0.025),0,lty=2,col=2)

plot(X2, stdres(fit.model2),pch=16, ylim=range(stdres(fit.model),-4,4))
abline(qnorm(0.975),0,lty=2,col=2)
abline(qnorm(0.025),0,lty=2,col=2)

plot(X3, stdres(fit.model2),pch=16, ylim=range(stdres(fit.model),-4,4))
abline(qnorm(0.975),0,lty=2,col=2)
abline(qnorm(0.025),0,lty=2,col=2)

# Pontos de alavanca
H = X %*% solve(t(X) %*% X) %*% t(X)
diag(H)

plot(diag(H),pch=16)
lim = 2*(4+1)/48 

abline(lim, 0, lty=2, col=2)

# Seleção de modelos
library(MASS)
stepAIC(fit.model)

# Mínimos quadrados generalizados
V = diag(X2)
X<-cbind(1,X1,X2,X3)
beta.hatG <- solve(t(X) %*% solve(V) %*% X) %*% t(X) %*% solve(V) %*% Y

# Mínimos quadrados generalizados - Forma alternativa
fit.modelG2 <- lm(Y ~ X1+X2+X3, weights = 1/X2)
beta.hatG2 <- fit.model$coefficients
beta.hatG2

# Gráfico de probabilidade normal - modelo no Caso 1
qqnorm(stdres(fit.model2),pch=16)    

# Envelope para o modelo normal - modelo no Caso 1
fit.model<-fit.model2
source("http://www.ime.usp.br/~giapaula/envel_norm") 

# Diagnóstico no caso normal
source("http://www.ime.usp.br/~giapaula/diag_norm")

#Identificando o maior resíduo*
which(stdres(fit.model) == max(stdres(fit.model)))