tv.dat<-scan("c:/temp/tv.dat",what=list(sub=0,homes=0,pcincm=0,instfee=0,svc=0,ncable=0,ntv=0))
attach(tv.dat)
sub<-tv.dat$sub
homes<-tv.dat$homes
pcincm<-tv.dat$pcincm
instfee<-tv.dat$instfee
svc<-tv.dat$svc
ncable<-tv.dat$ncable
ntv<-tv.dat$ntv
sub<-sub/1000
sub<-sqrt(sub)
homes<-homes/1000
tv<-data.frame(sub,homes,pcincm,instfee,svc,ncable,ntv)
Y<-cbind(sub,homes,pcincm,instfee,svc,ncable,ntv)
X<-cbind(homes,pcincm,instfee,svc,ncable,ntv)
cor(X)
pairs(Y)
###selecao do modelo baseado na soma extra de quadrados
tv.model<-lm(sub~homes+pcincm+instfee+svc+ncable+ntv,data=tv)
tv.model
lms<-summary(tv.model)
lms
anova(tv.model)
tv0.model<-lm(sub~homes+ncable+ntv,data=tv)
anova(tv0.model,tv.model)
summary(tv0.model)
##
   ###selecao do modelo pelo metodo do AIC
library(MASS)
##exemplo de sintaxe
stepAIC(tv.model)
stepAIC(tv.model,direction = c("backward"))
stepAIC(tv.model,scope= list(upper = ~homes+pcincm+instfee+svc+ncable+ntv, lower = ~1))
###selacao usando stepwise pelo AIC
saic<-stepAIC(tv.model,direction = c("both"))
model.aic<-lm(formula = sub ~ homes + ncable + ntv,  data = tv)
summary(model.aic)
####

fit.model<-model.aic
lms<-summary(fit.model)
n<-length(sub)
X<-model.matrix(fit.model)
p<-ncol(X)
s<-lms$sigma
r<-resid(lms)
h<-lm.influence(fit.model)$hat
ti<-r/(s*sqrt(1-h))
tsi<-ti*sqrt((n-p-1)/(n-p-ti^2))
yh<-fitted(fit.model)
###envelope####
form=sub~homes+ncable+ntv
envelope.normal(form,k=20,alfa=0.05)
win.graph()
par(mfrow=c(2,2))
#####grafico do residuos
plot(tsi,xlab="indice",ylab="Residuo Studentizado")
abline(-2,0,lty=2)
abline(2,0,lty=2)
#identify(tsi)
plot(yh,tsi,xlab="valores ajustados",ylab="Residuo Studentizado")
abline(-2,0,lty=2)
abline(2,0,lty=2)
#identify(yh,tsi)

#####graficos de diagnosticos
 win.graph()
par(mfrow=c(2,2))
plot(h,xlab="indice",ylab="Leverage")
abline(2*p/n,0,lty=2)
abline(3*p/n,0,lty=3)
identify(h,n=2)
di<-(ti^2)*h/(p*(1-h))
plot(di,xlab="indice",ylab="DCook")
identify(di)
dfit<-abs(tsi)*sqrt(h/(1-h))
plot(dfit,xlab="indice",ylab="DFFITS")
abline(2*sqrt(p/(n-p)),0,lty=2)
identify(dfit)
r<-resid(fit.model)
####grafico de influencia Local
R<-diag(r)
X<-model.matrix(fit.model)
H<-X%*%solve(t(X)%*%X)%*%t(X)
A<-R%*%H%*%R
Lmax<-eigen(A)$val[1]
dmax<-eigen(A)$vec[,1]
dmax<-dmax/sqrt(Lmax)
dmax<-abs(dmax)
plot(dmax,xlab="indice",ylab="dmax")
identify(dmax)
###################
###Calculando a variacao (%) quando retiramos o ponto
fit<-lm(sub~homes+ncable+ntv)
fit1<-lm(sub~homes+ncable+ntv,subset=-c(21))
sa<-100*((coef(fit1)-coef(fit)) /coef(fit)  )
sfit<-summary(fit)
sfit1<-summary(fit1)
sa1<-100*(sfit1$sigma*sqrt(diag(sfit1$cov.unscaled))-sfit$sigma*sqrt(diag(sfit$cov.unscaled)))/(sfit$sigma*sqrt(diag(sfit$cov.unscaled)))
round(sa,digits=2)
round(sa1,digits=2)