###Lendo arquivo babies
## O arquivo deve estar no seu diretorio de trabalho
babies <- read.table("babies.txt", header=T,na.strings=999)
## Corrigindo os outros dados com código de NA (9 ou 99, dependendo da variavel)
babies$age[babies$age==99] <- NA
babies$height[babies$height==99] <- NA
babies$smoke[babies$smoke==9] <- NA
babies$smoke <- as.logical(babies$smoke)
## Eliminando linhas com algum dado faltante
babies <- babies[apply(is.na(babies),1,sum)==0,]
summary(babies)


##Figura do Slide 1 ###
plot(bwt~gestation, data=babies, cex=0.5)
##Regressao em funcao do tempo de gestacao
babies.m1 <- lm(bwt~gestation, data=babies)
##Adicionando a linha da regressao ao grafico
abline(babies.m1, col="blue", lwd=2)

##Tabela de ANOVA da regressao
anova(babies.m1)

##### Figura do slide 2 #####
## Grafico demonstrativo dos residuos
x <- runif(8,0,10)
y <- rnorm(8,mean=2+3*x,sd=2.5)
m1 <- lm(y~x)
my <- mean(y)
##A figura
par(mfrow=c(1,2))
##Grafico com desvios em relacao aa media
plot(y~x, pch=19, col="blue", ylim=c(0,30), main="SS Total")
abline(h=mean(y), lty=2)
segments(x,rep(my,8),x,y,col="blue")
## Grafico com desvios em relacao a reta da regressao
plot(y~x, pch=19, col="blue", ylim=c(0,30), main="SS Erro")
abline(m1, lty=2)
segments(x,fitted(m1),x,y,col="blue")
par(mfrow=c(1,1))

################################################################################
### Calculos da ANOVA da Regressao
################################################################################
##Conteudo do objeto de regressao
names(babies.m1)
coef(babies.m1)#coeficientes
## Calculo dos previstos e dos residuos
## Guardando os coeficientes em um objeto numerico
cf.m1 <- as.numeric(coef(babies.m1))
## Valores previstos
babies$pred.m1 <- cf.m1[1]+cf.m1[2]*babies$gestation
## Residuos
babies$res.m1 <- babies$bwt - babies$pred
## Mostrando todos juntos
print(head(babies[,c(2,1,8,9)]),digits=3)
## A funcao residuals ja calcula os residuos
babies$res.m1[1:4]
residuals(babies.m1)[1:4]
## E tambem a funcao fitted para os previstos
babies$pred.m1[1:4]
fitted(babies.m1)[1:4]
## Soma dos quadrados total
babies.SST <- sum((babies$bwt-mean(babies$bwt))^2)
##Soma dos quadrados dos resi­duos (= soma dos quadrados do erro)
babies.m1.SSE <- sum(residuals(babies.m1)^2)
## Soma dos quadrados da regressao
babies.m1.SSR <- babies.SST - babies.m1.SSE
## Calculo do MSE
babies.m1.MSE <- babies.m1.SSE/(length(babies$bwt)-2)
## Valor de F
babies.m1.F <- babies.m1.SSR/babies.m1.MSE
## Conferindo com o resultado da anova
## Nossos calculos
##MSR (= SSR /1)
babies.m1.SSR
##MSE
babies.m1.MSE
##Valor de F
babies.m1.F
## Significancia do F
pf(babies.m1.F,df1=1, df2=(length(babies$bwt)-2), lower.tail=F)
## E agora a anova
anova(babies.m1)
## Coeficiente de determinacaoo: explicado / total
babies.m1.SSR/babies.SST
################################################################################

###Slide 6####
##Sintese: Objeto lm e seus extratores mais comuns:
## O objeto
names(babies.m1)
babies.m1$coefficients
babies.m1$residuals[1:4]
babies.m1$fitted.values[1:4]
babies.m1$call

###Slide 6####
## Funcoes extratoras
names(babies.m1)
coef(babies.m1)
confint(babies.m1)
residuals(babies.m1)[1:4]
fitted(babies.m1)[1:4]
logLik(babies.m1) ## pacote MASS
AIC(babies.m1)
summary(babies.m1)

######################################################################
####Demonstracao da linha de minimos quadrados
######################################################################

## Gerando os pontos: valor esperado de y = 2 + 3*x, com desvio=padrao de 2.5
## Oito valores de x tomados de uma distribuicao uniforme com minimo=0 e maximo=10
x1 <- runif(8,0,10)
## Valores de y tem media = 2+3*x e distribuicao normal
## em torno desta media, com desvio-padrao = 2.5
y1 <- rnorm(8,mean=2+3*x1,sd=2.5)
## Medias de x e y
my1 <- mean(y1)
mx1 <- mean(x1)

### Aqui feche todas as janelas graficas ###

##Primeiro Grafico: regressao com ponto na intersecao das medias
## "fulcro" da alavanca
## Abrindo a primeira janela para colocar este grafico
x11(height=4.5, width=5)
plot(y1~x1, col="blue", pch=19)
points(mx1,my1,cex=2, col="red"); points(mx1,my1,pch=19,cex=1/3,col="red")
## reta da regressao
abline(lm(y1~x1), col="red", lwd=2)

##Peparando o grafico com os pontos e "apoio da alavanca"
plot(y1~x1, col="blue", pch=19)
points(mx1,my1,cex=2, col="red"); points(mx1,my1,pch=19,cex=1/3,col="red")
## Preparando o grafico do SSE por valor da inclinacao
## Abrindo a segunda janela grafica
x11(height=3.5, width=4)
## Sequencia de valores de inclinacao a avaliar
intervalo <- seq(-6,15,by=0.5) 
##Valores maximos e minimos de SSE, para o intervalo de valores de inclinacoes
## Estes valores sao necessarios para definir a escala do eixo y
##Valor maximo de sse: o máximo dos extremos da sequencia
b.max <- max(intervalo)
a.max <- my1-b.max*mx1
b.min <- min(intervalo)
a.min <- my1-b.min*mx1
pred.max <- a.max+b.max*x1
pred.min <- a.min+b.min*x1
sse.max <- max(c(sum((y1-pred.max)^2),sum((y1-pred.min)^2)))
## Valor minimo
sse.min <- sum((y1-fitted(lmy1~x1))^2)
## O grafico
plot(0,0,xlim=range(intervalo),ylim=c(sse.min,sse.max),xlab="Inlinacao", ylab="SSE", type="n")

###Aqui comeca a animacao! ###
## A interacao: os SSEs sao calculados para valores de b definidos no intervalo
## Dentro de um loop. Para cada valor, os dois graficos sao atualizados
## Todos os valores
sse <- c()
i <- 0
devAskNewPage(TRUE)
for(b in intervalo){
  i <- i+1
  pred.y1 <- (my1-b*mx1)+b*x1
  sse[i] <- sum((y1-pred.y1)^2)
  dev.set(2)
  plot(y1~x1, col="blue", pch=19, ylim=c(-10,50))
  points(mx1,my1,cex=2, col="red"); points(mx1,my1,pch=19,cex=1/3,col="red")
  abline(my1-b*mx1,b, col="red", lty=2)
  segments(x1,pred.y1,x1,y1,col="blue")
  dev.set(3)
  points(b,sse[i],pch=19, cex=0.5)
}
devAskNewPage(F)
###Fim da animacao ###

##Qual o valor minimo?
## Pelo locator
locator(1)
b.est.g <- locator(1)$x
##Pelo mi­nimo dos valores
b.est <- intervalo[sse==min(sse)]
## Comparando
b.est.g
b.est
## Plotando no grafico
dev.set(2)
## Grafico com o "fulcro"
plot(y1~x1, col="blue", pch=19)
points(mx1,my1,cex=2, col="red"); points(mx1,my1,pch=19,cex=1/3,col="red")
##Linha da regressao ajustada pelo lm
abline(lm(y1~x1), lty=2, col="red")
## Linha com o valor de inclinação obtido do grafico
abline(my1-b.est*mx1,b.est, lty=2)



################################################################################
##Diagnosticos###
################################################################################
### Slide 11 ###
## Figura com os quatro graficos basicos
par(mfrow=c(2,2))
plot(babies.m1)
par(mfrow=c(1,1))

### Slide 12 ###
##Primeiro grafico: residuos x estimados
##Exemplo de relacao que nao segue uma reta
x <- runif(30,5,10)
y <- rnorm(30,mean=2*exp(x/1.8),sd=25)
##Grafico de residuos
par(mfrow=c(1,2))
plot(lm(y~x),which=1)
plot(y~x)
abline(lm(y~x),col="blue")
par(mfrow=c(2,2))
###Slide 13 ###
##Exemplo de residuos nao normais
x <- runif(100,1,10)
y <- rexp(100,rate=2/x)
###Slide 14####
## Grafico qqnorm dos residuos
par(mfrow=c(1,2))
plot(lm(y~x),which=2)
plot(y~x)
abline(lm(y~x),col="blue")
par(mfrow=c(1,1))
###Slide 15###
##Exemplo de variancia que aumenta com a media
par(mfrow=c(1,3))
plot(lm(y~x),which=1, add.smooth=F, pch=19)
plot(lm(y~x),which=3, pch=19)
plot(y~x, pch=19)
abline(lm(y~x),col="blue")
par(mfrow=c(1,1))
###Slide 16 ###
##Exemplo de ponto influente
x <- c(rnorm(20),10,15,19)
y <- c(rnorm(20,mean=2),rnorm(3,mean=c(20,30,38),sd=15))
par(mfrow=c(1,3))
plot(lm(y~x),which=c(4,5),add.smooth=F, pch=19)
plot(y~x, pch=19)
abline(lm(y~x),col="blue")
par(mfrow=c(1,1))


################################################################################
## Exemplo de analise: massa cerebral e corporal de 28 animais
################################################################################
## o objeto Animals e do pacote MASS
library(MASS)
##Uma copia no workspace
data(Animals)
## O objeto 
Animals
##Grafico
plot(brain~body,data=Animals)
## Claramente nao linear, tentando em escala log
plot(brain~body,data=Animals, log="xy")
plot(log(brain)~log(body), data=Animals)
##Ok, mas ha 3 pontos fora do padrao. Quem sao?
Animals[log(Animals$body)>8&log(Animals$brain)<6,] # sao os dinossauros!
## Ok, vamos ajustar uma regressao linear com os dados transformados para log
anim.m1 <- lm(log(brain)~log(body),data=Animals)
## E vamos adicionar a linha ao grafico
abline(anim.m1, col="blue")
## Graficos de diagnostico (um a um)
plot(anim.m1)
## Varios problemas. Vamos verificar como fica sem os dinossauros
anim.m2 <- lm(log(brain)~log(body),data=Animals, subset=!(log(Animals$body)>8&log(Animals$brain)<6))
##Grafico com as retas dos dois modelos
plot(log(brain)~log(body), data=Animals)
abline(anim.m1, col="blue", lty=2)
abline(anim.m2, col="red")
##Novos Diagnosticos (4 graficos de uma vez)
par(mfrow=c(2,2))
plot(anim.m2)
par(mfrow=c(1,1))
## Agora humanos e macacos parecem ser influentes, mas vamos em frente
## ANOVA do modelo
anova(anim.m2)
## Resumo do modelo
summary(anim.m2)
## Coeficientes e intervalos de confianca
coef(anim.m2)
confint(anim.m2)
## Agora vamos acrescentar e planiha transformar os valores previstos
## de volta na escala original
cf <- as.numeric(coef(anim.m2))
Animals$prev <- round(exp(cf[1])*Animals$body^cf[2],1)
Animals

################################################################################
### Tutorial 2: Uma simulacao com dados fortemente nao normais
################################################################################
## Vamos criar uma semente para os numeros aletorios padronizada,
## para que todos tenham os mesmos resultados
set.seed(22)
##Valores da preditora: 30 valores tomados de uma distribuicao uniforme
x <- runif(30,0,10)
##Variavel resposta: tem distribuicao exponencial, com media = 2+5*preditora
## Nota: o parametro rate da exponencial e o inverso da media
media.y = 2+5*x
y <- rexp(30,rate=1/media.y)
## Um grafico da distribuicao exponencial com media =10
curve(dexp(x,rate=1/10),0,50, ylab="Densidade", lwd=2, col="blue")
## Acrescentando a media
abline(v=10, lty=2)
abline
## Um grafico
plot(y~x)
##Ajuste do modelo
m1 <- lm(y~x)
## Diagnostico
par(mfrow=c(2,2))
plot(m1)
par(mfrow=c(1,1))
## Todos os problemas!
## Quais os coeficientes estimados?
coef(m1)## Bem longe do valor real
## Agora vamos simular 2000 repeticoes desta amostra de y,
## ajustando o modelo para cada uma delas
## e guardando os valores dos coeficientes em cada uma
## e tambem a se os intervalos de confianca dos coeficientes incluem os valores verdadeiros
## Uma matriz para guardar os resultados
result <- matrix(ncol=6,nrow=2000) #cada linha guarda o resultado de uma randomizacao

##As simulacoes
for(i in 1:2000){
  yr <- rexp(30,rate=1/(2+5*x))
  m <- lm(yr~x)
  CIs <- confint(m)
  result[i,1:2] <- as.numeric(coef(m))# guarda os coeficientes nas duas primeiras colunas
  result[i,3:4] <- CIs[1,]#guarda limites dos ICs do intercepto
  result[i,5:6] <- CIs[2,]#guarda limites dos ICs da inclinacao
}

## Agora vamos ver as distribuicao dos valores dos interceptos estimados
## que estao guardadas na primeira coluna da matriz de resultados
hist(result[,1])
##Acrescentamos a media destes valores
abline(v=mean(result[,1]), lty=2, col="blue")
## E o valor real de a = 2
abline(v=2, col="red", lty=2)
## Agora vamos ver as distribuicao dos valores das inclinacoes estimadas
## que estao guardadas na segunda coluna da matriz de resultados
hist(result[,2])
##Acrescentamos a media destes valores
abline(v=mean(result[,2]), lty=2, col="blue")
## E o valor real de b = 5
abline(v=5, col="red", lty=2)
## Nos dois casos o vies e muito pequeno.
## Vamos calcular este vies  em percentual em relacao ao valor real
##Do intercepto
100*(2-mean(result[,1]))/2 ## -3%
##Da inclinacao
100*(5-mean(result[,2]))/5 ## 1%


##Agora vamos verificar as estimativas pelo intervalo de confianca
## Quantas vezes os intervalos estimados incluíram o valor real?
##Para o intercepto:
sum(result[,3]<=2&result[,4]>=2)
## Em 1963 das 2000 simulacoes o intervalo incluiu a média, ou seja
## a significancia do intervalo e de 1/2000 = 0,0005
## menor do que a significância nominal de 0,05

##Para a inclinacao:
sum(result[,5]<=5&result[,6]>=5)
## Aqui as simulacoes mostram que o intervalo de confianca teve uma siginificancia
## um pouco maior que a nominal (alfa = 0,086)

###################################################################
## Exerci­cios
###################################################################
## Altura aos dois anos
## em polegadas
dois.anos <- c(39,30,32,34,35,36,36,30)
adulto <- c(71,63,63,67,68,68,70,64)
anova(lm(dois.anos~adulto))
confint(lm(dois.anos~adulto))
plot(dois.anos~adulto)

## Iris
head(iris)
rm(iris)
setosa <- iris[iris$Species=="setosa",]
sepalL <- residuals(lm(Sepal.Length~Petal.Length,data=setosa))
sepalW <- residuals(lm(Sepal.Width~Petal.Length,data=setosa))
anova(lm(sepalL~sepalW))