############################### Aula Prática 2 ##########################################
############################### Modelo para dados de preço ##############################

 detach(dados)

 dados<-read.table("http://wiki.icmc.usp.br/images/1/12/Pre%C3%A7os.txt", header=TRUE)

 attach(dados)

 X <- cbind(1, X1, X2, X3, X4, X5, X6, X7, X8, X9)
 
 lm(Y ~ -1 +X) # por que colocar -1 ?

 fit.model <- lm(Y ~ -1 +X)

 library(MASS)

 stepAIC(lm(Y ~ X1+X2+X3+X4+X5+X6+X7+X8+X9), direction = "both")
 
 fit.model <- lm(Y ~ X1+X2+X3+X4+X5+X6+X7+X8+X9)

 source("http://www.ime.usp.br/~giapaula/diag_norm")

 source("http://www.ime.usp.br/~giapaula/envel_norm")



############################### Modelo para dados de censo ##############################


dados<-read.table("http://wiki.icmc.usp.br/images/c/c0/Censo.txt", header=TRUE)

attach(dados)

Y <- renda

X <- escolaridade

fit.model <- lm(renda~escolaridade)

 source("http://www.ime.usp.br/~giapaula/diag_norm")

 source("http://www.ime.usp.br/~giapaula/envel_norm")


# Transformação em X

 plot(escolaridade^2, renda)

 par(new=TRUE)

 plot(lowess(escolaridade^2,renda), type="l", ylim=range(renda), Xlim=range(escolaridade^2), ylab="", Xlab="",  lty=2)

 fit.model <- lm(renda~escolaridade^2)

 source("http://www.ime.usp.br/~giapaula/diag_norm")

 source("http://www.ime.usp.br/~giapaula/envel_norm")



# Transformação em Y

 plot(escolaridade, log(renda))

 par(new=TRUE)

 plot(lowess(escolaridade, log(renda)), type="l", ylim=range(log(renda)), Xlim=range(escolaridade), ylab="", Xlab="",  lty=2)
 
 fit.model <- lm(log(renda)~escolaridade)

 source("http://www.ime.usp.br/~giapaula/diag_norm")

 source("http://www.ime.usp.br/~giapaula/envel_norm")


# Modelo gama

 fit.model = glm(renda ~ escolaridade, family=Gamma(link=log))

 summary(fit.censo)

 source("http://www.ime.usp.br/~giapaula/diag_gama")

 source("http://www.ime.usp.br/~giapaula/envel_gama")



############################### Modelo para dados de poliamina ##############################

# Transformação na variável resposta
# Modelo original Y ~ X

rm(list=ls(all=TRUE))

dados<-read.table("http://wiki.icmc.usp.br/images/d/d0/Poliamina.txt", header=TRUE)

attach(dados)

par(mfrow=c(1,3))

plot(X,Y, pch=16)

abline(lm(Y~X), lwd=2)

par(new=TRUE)

plot(lowess(X,Y), type="l", ylim=range(Y), xlim=range(X), ylab="", xlab="",  lty=2)

legend(locator(1), c("modelo linear", "lowess"), lwd=c(2,1), lty=c(1,2))

fit.model <-lm(Y~X)

 source("http://www.ime.usp.br/~giapaula/diag_norm")

 source("http://www.ime.usp.br/~giapaula/envel_norm")




# Transformação na variável preditora
# Modelo original log(Y) ~ X

rm(list=ls(all=TRUE))

dados<-read.table("http://wiki.icmc.usp.br/images/d/d0/Poliamina.txt", header=TRUE)

attach(dados)

par(mfrow=c(1,3))

plot(X,log(Y), pch=16)

abline(lm(log(Y)~X), lwd=2)

par(new=TRUE)

plot(lowess(X,log(Y)), type="l", ylim=range(log(Y)), xlim=range(X), ylab="", xlab="",  lty=2)

legend(locator(1), c("modelo linear", "lowess"), lwd=c(2,1), lty=c(1,2))

fit.model <-lm(log(Y)~X)


 source("http://www.ime.usp.br/~giapaula/diag_norm")

 source("http://www.ime.usp.br/~giapaula/envel_norm")


# Transformação de Box-Cox

rm(list=ls(all=TRUE))

dados<-read.table("http://wiki.icmc.usp.br/images/d/d0/Poliamina.txt", header=TRUE)

attach(dados)

boxcox(Y ~ X, lambda = seq(-2, 2, len = 20))

rm(list=ls(all=TRUE))

par(mfrow=c(1,3))

plot(X,1/sqrt(Y), pch=16)

abline(lm(1/sqrt(Y) ~ X), lwd=2)

par(new=TRUE)

plot(lowess(X,1/sqrt(Y)), type="l", ylim=range(log(Y)), xlim=range(X), ylab="", xlab="",  lty=2)

legend(locator(1), c("modelo linear", "lowess"), lwd=c(2,1), lty=c(1,2))

fit.model <-lm(1/sqrt(Y)~X)


 source("http://www.ime.usp.br/~giapaula/diag_norm")

 source("http://www.ime.usp.br/~giapaula/envel_norm")



############################### Modelo para dados de programação ##############################


dados<-read.table("http://wiki.icmc.usp.br/images/d/d4/Programa%C3%A7%C3%A3o.txt", header=TRUE)

attach(dados)

plot(tempo, sucesso, pch=16, xlim=c(0,40), ylim=c(0,1))

fit.model <- glm(sucesso ~ tempo, family=binomial(link=logit))

par(new=T)

curve(exp(-3.0597+ 0.1615 *x)/(1+exp(-3.0597+ 0.1615 *x)), 0, 40, xlim=c(0,40), ylim=c(0,1), xlab="", ylab="")


 source("http://www.ime.usp.br/~giapaula/diag_bino")

 source("http://www.ime.usp.br/~giapaula/envel_bino")