setwd("/Users/hlopes/Desktop/TEACHING/2016-2to6-Econometria-graduacaoeconomia/aulas/aula0330")

temco <-read.table("http://hedibert.org/wp-content/uploads/2014/04/temco.txt",header=TRUE)

attach(temco)

n = nrow(temco)

# creating dummies variables for department
dept4 = rep(0,n)
dept4[DEPT==4]=1


# p regressors: 2^p-1 linear models
# p=5: 31 linear models

reg1=lm(SALARIO~EDUC)
reg2=lm(SALARIO~ANOSEMP)
reg3=lm(SALARIO~EXPPREV)
reg4=lm(SALARIO~SEXO)
reg5=lm(SALARIO~dept4)
reg6=lm(SALARIO~EDUC+ANOSEMP)
reg7=lm(SALARIO~EDUC+EXPPREV)
reg8=lm(SALARIO~EDUC+SEXO)
reg9=lm(SALARIO~EDUC+dept4)
reg10=lm(SALARIO~ANOSEMP+EXPPREV)
reg11=lm(SALARIO~ANOSEMP+SEXO)
reg12=lm(SALARIO~ANOSEMP+dept4)
reg13=lm(SALARIO~EXPPREV+SEXO)
reg14=lm(SALARIO~EXPPREV+dept4)
reg15=lm(SALARIO~SEXO+dept4)
reg16=lm(SALARIO~EDUC+ANOSEMP+EXPPREV)
reg17=lm(SALARIO~EDUC+ANOSEMP+SEXO)
reg18=lm(SALARIO~EDUC+ANOSEMP+dept4)
reg19=lm(SALARIO~EDUC+EXPPREV+SEXO)
reg20=lm(SALARIO~EDUC+EXPPREV+dept4)
reg21=lm(SALARIO~EDUC+SEXO+dept4)
reg22=lm(SALARIO~ANOSEMP+EXPPREV+SEXO)
reg23=lm(SALARIO~ANOSEMP+EXPPREV+dept4)
reg24=lm(SALARIO~ANOSEMP+SEXO+dept4)
reg25=lm(SALARIO~EXPPREV+SEXO+dept4)
reg26=lm(SALARIO~EDUC+ANOSEMP+EXPPREV+SEXO)
reg27=lm(SALARIO~EDUC+ANOSEMP+EXPPREV+dept4)
reg28=lm(SALARIO~EDUC+ANOSEMP+SEXO+dept4)
reg29=lm(SALARIO~EDUC+EXPPREV+SEXO+dept4)
reg30=lm(SALARIO~ANOSEMP+EXPPREV+SEXO+dept4)
reg31=lm(SALARIO~EDUC+ANOSEMP+EXPPREV+SEXO+dept4)


r2=rbind(c(summary(reg1)$r.sq,summary(reg1)$adj.r),
c(summary(reg2)$r.sq,summary(reg2)$adj.r),
c(summary(reg3)$r.sq,summary(reg3)$adj.r),
c(summary(reg4)$r.sq,summary(reg4)$adj.r),
c(summary(reg5)$r.sq,summary(reg5)$adj.r),
c(summary(reg6)$r.sq,summary(reg6)$adj.r),
c(summary(reg7)$r.sq,summary(reg7)$adj.r),
c(summary(reg8)$r.sq,summary(reg8)$adj.r),
c(summary(reg9)$r.sq,summary(reg9)$adj.r),
c(summary(reg10)$r.sq,summary(reg10)$adj.r),
c(summary(reg11)$r.sq,summary(reg11)$adj.r),
c(summary(reg12)$r.sq,summary(reg12)$adj.r),
c(summary(reg13)$r.sq,summary(reg13)$adj.r),
c(summary(reg14)$r.sq,summary(reg14)$adj.r),
c(summary(reg15)$r.sq,summary(reg15)$adj.r),
c(summary(reg16)$r.sq,summary(reg16)$adj.r),
c(summary(reg17)$r.sq,summary(reg17)$adj.r),
c(summary(reg18)$r.sq,summary(reg18)$adj.r),
c(summary(reg19)$r.sq,summary(reg19)$adj.r),
c(summary(reg20)$r.sq,summary(reg20)$adj.r),
c(summary(reg21)$r.sq,summary(reg21)$adj.r),
c(summary(reg22)$r.sq,summary(reg22)$adj.r),
c(summary(reg23)$r.sq,summary(reg23)$adj.r),
c(summary(reg24)$r.sq,summary(reg24)$adj.r),
c(summary(reg25)$r.sq,summary(reg25)$adj.r),
c(summary(reg26)$r.sq,summary(reg26)$adj.r),
c(summary(reg27)$r.sq,summary(reg27)$adj.r),
c(summary(reg28)$r.sq,summary(reg28)$adj.r),
c(summary(reg29)$r.sq,summary(reg29)$adj.r),
c(summary(reg20)$r.sq,summary(reg30)$adj.r),
c(summary(reg31)$r.sq,summary(reg31)$adj.r))


par(mfrow=c(1,1))
plot(r2[,1],ylim=range(r2),xlab="Model",ylab="R2",pch=16,axes=FALSE)
axis(2);axis(1,at=1:31);box()
points(r2[,2],pch=16,col=2)

plot(sort(r2[,2]),axes=FALSE,xlab="Model",ylab="Adjusted R2",pch=16)
axis(2);axis(1,at=1:31,lab=order(r2[,2]));box()

          (Intercept)         EDUC      ANOSEMP      EXPPREV         SEXO        dept4         R2   Adjusted R2
Model 27     24620.07      1849.26       700.81        65.29                  -7806.79  0.8257088     0.8087048
Model 31     26053.12      1821.36       683.00        16.91     -1802.55     -7911.12  0.8257088     0.8087048
Model 18     24867.80      1858.30       699.60                               -7714.30  0.8252351     0.8127519
Model 28     26134.90      1823.20       682.40                  -1827.90     -7889.60  0.8318788     0.8154767


par(mfrow=c(1,2))
plot(reg1$fit,SALARIO,xlab="Fitted",ylab="Salary",main="Best model with one regressor\n EDUC",xlim=c(20000,80000),ylim=c(20000,80000))
abline(0,1,col=2)
abline(-2*summary(reg1)$sig,1,col=2,lty=2)
abline(2*summary(reg1)$sig,1,col=2,lty=2)
plot(reg18$fit,SALARIO,xlab="Fitted",ylab="Salary",main="Best overall model\n EDUC,ANOSEMP,DEPT4",xlim=c(20000,80000),ylim=c(20000,80000))
abline(0,1,col=2)
abline(-2*summary(reg18)$sig,1,col=2,lty=2)
abline(2*summary(reg18)$sig,1,col=2,lty=2)