### Regressao multipla ###
##Leitura dos dados, obtidos de http://www.bio.ic.ac.uk/research/mjcraw/statcomp/data/
## Arquivo texto separado por tabulacao (argumento sep="\t")
poluicao <- read.table("Pollute.txt", header=T, sep="\t")

################################################################################
## 1. Primeiro modelo: efeito do numero de industrias
################################################################################

## 1.1.Graficos Exploratorios

par(mfrow=c(2,2))
plot(Pollution~Industry, data=poluicao)
plot(Pollution~Population, data=poluicao)
plot(Pollution~Temp, data=poluicao)
plot(Pollution~Wind, data=poluicao)
par(mfrow=c(1,1))

##1.2. Ajuste e avaliacao do modelo
pol.m1 <- lm(Pollution~Industry, data=poluicao)
##acrescentando a linha no grafico
plot(Pollution~Industry, data=poluicao)
abline(pol.m1, col="blue")
## A cidade mais industrializada tem uma alavancagem forte?
abline(lm(Pollution~Industry, data=poluicao, subset=Industry<max(Industry)), col="red")
## Nao tem. Vamos seguir em frente.
## Teste deste modelo
anova(pol.m1)
## Ha uma comparacao embutida neste teste, com o modelo:
pol.m0 <- lm(Pollution~1, data=poluicao)
##compare:
anova(pol.m1,pol.m0)
anova(pol.m1)
##Que raios e isto? Vamos ver o resumo do modelo:
summary(pol.m0)
## Duas dicas:
mean(poluicao$Pollution)
sd(poluicao$Pollution)
## Grafico dos dois modelos concorrentes
plot(Pollution~Industry, data=poluicao)
abline(pol.m1, col="blue")
abline(h=mean(poluicao$Pollution),col="red")
## Entendeu?

################################################################################
## 2. Complicando o modelo: efeito de variaveis climáticas
################################################################################

## 2.1. Ajuste dos modelo e comparando-os
## Segundo modelo: acrescentamos efeitos de temperatura
pol.m2 <- lm(Pollution~Industry+Temp, data=poluicao)
## Comparacao dos dois modelos
anova(pol.m1,pol.m2)
## Ha um ganho significativo de variancao explicada com a inclusao da temperatura
##O ganho pela adicao de cada variavel, NA ORDEM DE INCLUSAO eh o resultado do comando:
anova(pol.m2)
## Prosseguindo, vamos adicionar o efeito do vento
pol.m3 <- update(pol.m2,.~.+Wind)
anova(pol.m2,pol.m3)
anova(pol.m3)
## Sem efeitos, ficamos com o modelo 2.

## 2.2.Graficos de diagnostico
par(mfrow=c(2,2))
plot(pol.m2)
par(mfrow=c(1,1))
## Quem eh o ponto 41?
poluicao[41,c(1:3,5)]
## Nos graficos
par(mfrow=c(1,2))
plot(Pollution~Industry, data=poluicao)
points(Pollution~Industry, data=poluicao[41,], col="red", pch=19)
plot(Pollution~Temp, data=poluicao)
points(Pollution~Temp, data=poluicao[41,], col="red", pch=19)
par(mfrow=c(1,1))

## Ok, resumo do modelo
summary(pol.m2)

## 2.3. Pausa importante
###IMPORTANTE###
## Importante compreender a diferenca dos testes no resumo e em
anova(pol.m2)
##O comando ''anova'' sempre compara modelos. Se o argumneto é um único modelo, ele retorna o resultado de um **teste sequencial**, isto é, o aumento de variação explicada por um fator em relação a um modelo com os fatores **acima dele na tabela**, que segue a ordem da fórmula do modelo. Por isso, os dois comandos abaixo não retornam exatamente os mesmos resultados:
anova(lm(Pollution~Temp+Wind+Industry, data=poluicao))
anova(lm(Pollution~Industry+Wind+Temp, data=poluicao))
##O teste para Industry no resumo equivale a
anova(lm(Pollution~Temp+Industry, data=poluicao),lm(Pollution~Temp, data=poluicao))
## E para temp eh:
anova(lm(Pollution~Temp+Industry, data=poluicao),lm(Pollution~Industry, data=poluicao))


################################################################################
## 3. Representacao grafica do modelo multiplo (uma coisa dificil!)
################################################################################

##3.1. Grafico do plano da regressao (cuidado, falsa perspectiva engana muito nossa visão)
## 3.1.1. Etapas preparatorias (nao vimso isto em aula, nao se assute com estes códigos
## Matriz com combinacoes de valores das duas preditoras a intervalos regulares
pred.m2 <- expand.grid(Industry=seq(min(poluicao$Industry),max(poluicao$Industry),length=10),
                       Temp=seq(min(poluicao$Temp),max(poluicao$Temp),length=10))
##Valores previstos
fit.m2 <- matrix(predict(pol.m2,pred.m2),10,10)
## Inicio e fim do arquivo
head(pred.m2)
tail(pred.m2)
## 3.1.2. Os Graficos
## Um grafico de superficie
persp(x = seq(min(poluicao$Industry),max(poluicao$Industry),length=10),
      y= seq(min(poluicao$Temp),max(poluicao$Temp),length=10),
      z=fit.m2, theta = -30, phi = 30, expand = 0.5, col = "lightblue",
      xlab="Industries", ylab="Temp", zlab="Esperado")
## Curvas de nivel (isolinhas)
contour(x = seq(min(poluicao$Industry),max(poluicao$Industry),length=10),
        y= seq(min(poluicao$Temp),max(poluicao$Temp),length=10),
        z=fit.m2,xlab="Industries", ylab="Temp", labcex=1.5)

## 3.2. Plotando previstos x observados
##Uma solucao eh dividir os  pontos em subconjuntos por uma das variaveis e plotar a outra
##Coeficientes do modelo
cf.m2 <- coef(pol.m2)
##Tercis da variavel Temp
Temp.tercis <- quantile(poluicao$Temp,c(0,1/3,2/3,1))
Temp.tercis
##Criando um fator com estes breakpoints, com a funcao cut
poluicao$Temp.c <- cut(poluicao$Temp, breaks=Temp.tercis)
## Verificando: tres fatores ordenados de faixas de temperatura
levels(poluicao$Temp.c)
## E agora vamos tres graficos com os valores selecionados po faixa de temperatura
par(mfrow=c(1,3))
##Primeiro Grafico: primeiro nivel do fator faixa de temperatura
plot(Pollution~Industry, data=poluicao, subset=Temp.c==levels(poluicao$Temp.c)[1],
     main=paste("Temp: ",Temp.tercis[1]," a ",Temp.tercis[2],sep=""),
     xlim=c(min(poluicao$Industry),max(poluicao$Industry)),
     ylim=c(min(poluicao$Pollution),max(poluicao$Pollution)) )
##Linhas do previsto para ponto medio da faixa de temperatura
abline(cf.m2[1]+cf.m2[3]*mean(Temp.tercis[1:2]),cf.m2[2], col="blue", lwd=2)
##Segundo Grafico
plot(Pollution~Industry, data=poluicao, subset=Temp.c==levels(poluicao$Temp.c)[2],
     main=paste("Temp: ",round(Temp.tercis[2],1)," a ",round(Temp.tercis[3],1),sep=""),
     xlim=c(min(poluicao$Industry),max(poluicao$Industry)),
     ylim=c(min(poluicao$Pollution),max(poluicao$Pollution)) )
abline(cf.m2[1]+cf.m2[3]*mean(Temp.tercis[2:3]),cf.m2[2], col="blue", lwd=2)
##Terceiro grafico
plot(Pollution~Industry, data=poluicao, subset=Temp.c==levels(poluicao$Temp.c)[3],
     main=paste("Temp: ",round(Temp.tercis[3],1)," a ",round(Temp.tercis[4],1),sep=""),
     xlim=c(min(poluicao$Industry),max(poluicao$Industry)),
     ylim=c(min(poluicao$Pollution),max(poluicao$Pollution)) )
abline(cf.m2[1]+cf.m2[3]*mean(Temp.tercis[3:4]),cf.m2[2], col="blue", lwd=2)
par(mfrow=c(1,1))



################################################################################
## 4 . TESTANDO INTERACOES NO MODELO DA POLUICAO
################################################################################
## Verificando graficamente a interacao entre industria e vento
## definimos 4 faixas de temperatura parcialmente sobrepostas
## com a funcao equal.count e condicionamos o grafico por esta variavel categorica
## o argumento panel e uma funcao a ser aplicada em cada painel, no caso plotar os pontos
## e fazer uma regressao apenas com estes pontos (detalhes na ajuda do xyplot)
temp.shingle <- equal.count(poluicao$Temp,number=4)
summary(temp.shingle)
xyplot(Pollution~Industry|temp.shingle,data=poluicao,
   panel=function(x,y,...){
     panel.xyplot(x,y,...)
     panel.abline(lm(y~x),...)
   }
       )

## O mesmo para Wind
xyplot(Pollution~Industry|equal.count(poluicao$Wind,number=4),data=poluicao,
   panel=function(x,y,...){
     panel.xyplot(x,y,...)
     panel.abline(lm(y~x),...)
   }
       )
## Os dois graficos sugerem que ha interacao entre as variaveis climaticas e o numero de industrias


##Vamos acrescentar estas interacoes nos modelos
pol.m4 <- update(pol.m2,.~.+Industry:Temp)## acrescenta apenas a interacao
pol.m5 <- update(pol.m2,.~.+Wind+Industry:Wind)##acrescenta a interacao e seus componentes
##Comparacao destes novos modelos com o anterior
anova(pol.m2,pol.m4)
anova(pol.m2,pol.m5)
## nenhuma das interacoes melhora o modelo.



################################################################################
## 5 . O FANTASMA DA COLINEARIDADE
################################################################################

## 5.1. Ainda o modelo com poluicao
## Inspecionado ralacoes entre todas as variaveis
pairs(poluicao[,1:7])
## Ha forte correlacao entre populacao e n de indústrias.
## Mas vamos acrescentar populacao ao modelo mesmo assim
pol.m6 <- update(pol.m2,.~.+Population)
anova(pol.m2,pol.m6)
##Resumo do modelo
summary(pol.m6)
## Populacao com efeito negativo ????
## Como ela se comporta sozinha?
pol.m7 <- lm(Pollution~Population, data=poluicao)
summary(pol.m7) ## xiii ...


## 5.2.  Um exemplo extremo de colinearidade
## gerando os dados: visitantes em flores
altura <- runif(50,20,100)
nflores <- rnorm(50,mean= 10 + 3*altura, sd=10)
med.vis <- 2 + 0.2*nflores + 0.5*altura
nvisitas <- rnorm(50,mean=med.vis, sd=15)
## relacoes entre as variaveis
pairs(data.frame(altura,nflores,nvisitas))
##Modelo com nflores
vis.m1 <- lm(nvisitas~nflores)
## Com altura
vis.m2 <- lm(nvisitas~altura)
# Com as duas
vis.m3 <- lm(nvisitas~altura+nflores)
## Comparacoes
anova(vis.m1,vis.m3)
anova(vis.m2,vis.m3)
summary(vis.m3)
anova(vis.m1)
summary(vis.m1)
anova(vis.m2)
summary(vis.m2)
anova(lm(nvisitas ~ altura + nflores))
anova(lm(nvisitas ~ nflores + altura ))

################################################################################
##INTERVALO
################################################################################

################################################################################
## 6. A MATRIZ DO MODELO
################################################################################
## Um exemplo de matriz de modelo de regressao multipla:
## Nosso modelo selecionado para o conjunto de dados poluicao:
head(model.matrix(pol.m2))
## Compare com os primeiro valores das preditoras, na matriz
head(poluicao[,2:3])
##Multiplicacao matricial da matriz do modelo pelos coeficientes
prod.m2 <- model.matrix(pol.m2)%*%coef(pol.m2)
##Primeiro valores deste produto
prod.m2[1:5]
## Valores ajustados
fitted(pol.m2)[1:5]

################################################################################
## 7. REGRESSOES COM VARIAVEIS CATEGORICAS
################################################################################
## Conjunto de dados "babies"
## Qual o efeito do tempo de gestacao (continua) e do fatod a mae fumar (fator) sobre o peso do bebe?
##7.1. Modelo sem interacao
babies.m2 <- lm(bwt~gestation+smoke,data=babies)
summary(babies.m2)
cf.bm2 <- coef(babies.m2)
cf.bm2
head(model.matrix(babies.m2))
head(babies[,c(2,7)])
## O previsto para o primeiro valor: primeiro fazemos o produto de cada coeficiente
## com a primeira linha da matriz do modelo
model.matrix(babies.m2)[1,]
cf.bm2
model.matrix(babies.m2)[1,]*cf.bm2
sum(model.matrix(babies.m2)[1,]*cf.bm2)
## Conferindo
fitted(babies.m2)[1] ## ok!
## O mesmo para a terceira linha
model.matrix(babies.m2)[3,]
cf.bm2
model.matrix(babies.m2)[3,]*cf.bm2
sum(model.matrix(babies.m2)[3,]*cf.bm2)
## Conferindo
fitted(babies.m2)[3] ## ok!
## Logo, temos duas retas que diferem quanto ao intercepto
## Um grafico
plot(bwt~gestation, data=babies, subset=smoke==F, xlim=range(babies$gestation),
     ylim=range(babies$bwt), col="blue", cex=0.5)
abline(cf.bm2[1],cf.bm2[2],lwd=2, col="blue")
points(bwt~gestation, data=babies, subset=smoke==T, col="red", cex=0.5)
abline(cf.bm2[1]+cf.bm2[3],cf.bm2[2],lwd=2, col="red")

## 7.2. Modelo com interacao
### E o que significa Interacao? ###
babies.m3 <- update(babies.m2,.~.+gestation:smoke)
head(model.matrix(babies.m3))
cf.bm3 <- coef(babies.m3)
cf.bm3
## Qual e o valor previsto para a primeira observacao agora?
model.matrix(babies.m3)[1,]
model.matrix(babies.m3)[1,]*cf.bm3
sum(model.matrix(babies.m3)[1,]*cf.bm3)
## Conferindo
fitted(babies.m3)[1] ## ok!
## E para a terceira?
cf.bm3
model.matrix(babies.m3)[3,]
model.matrix(babies.m3)[3,]*cf.bm3
sum(model.matrix(babies.m3)[3,]*cf.bm3)
## Conferindo
fitted(babies.m3)[3] ## ok!
## Ou seja, agora temos duas retas com inclinacoes diferentes!
plot(bwt~gestation, data=babies, subset=smoke==F, xlim=range(babies$gestation),
     ylim=range(babies$bwt), col="blue", cex=0.5)
abline(cf.bm3[1],cf.bm3[2],lwd=2, col="blue")
points(bwt~gestation, data=babies, subset=smoke==T, col="red", cex=0.5)
abline(cf.bm3[1]+cf.bm3[3],cf.bm3[2]+cf.bm3[4],lwd=2, col="red")
## O mesmo, em dois graficos separados
xyplot(bwt~gestation|smoke,data=babies,
       panel=function(x,y,...){
         panel.xyplot(x,y,...)
         panel.abline(lm(y~x),...)
                    })
## Falta verificar se a interacao e significativa
anova(babies.m2,babies.m3) ## Sim!
## Diagnosticos do modelo
par(mfrow=c(2,2))
plot(babies.m3)
par(mfrow=c(1,1))
## resumo do modelo
summary(babies.m3)
## Intervalos de confianca dos coeficientes
confint(babies.m3)

################################################################################
##8. ANOVA REVISITADA
################################################################################
## Altura das mudas de duas especies de arvores, em funcao de substrato e especie
## A pergunta eh: ha diferenças entre 
## Uma ANOVA de dois fatores
## 8.1. Leitura e conversao da variaveis de  bloco ( que não usamos),
## e substrato, que sao numerica, em fator:
mudas <- read.csv("altura-mudas.csv")
mudas$bloco <- as.factor(mudas$bloco)
mudas$substrato <- as.factor(mudas$substrato)

## 8.2. Primeiro modelo: altura depende da especie
mudas.m1 <- lm(altura~especie, data=mudas)
## Resumo do modelo
summary(mudas.m1)
## Coeficientes do modelo
cf.mud1 <- coef(mudas.m1)
cf.mud1
## o que siginificam estes coeficientes?
## A resposta está na matriz do modelo
head(model.matrix(mudas.m1))
##Vamos usa-la para o calculo dos esperados
head(model.matrix(mudas.m1)%*%cf.mud1)
tail(model.matrix(mudas.m1)%*%cf.mud1)
## Jah entendeu o que sao os coeficientes?
## Isto pdoe ajudar: veja a medias das especies
tapply(mudas$altura,mudas$especie,mean)
## De novo o resumo
summary(mudas.m1)
## Entendeu?

## 8.3. Anova de dois fatores
mudas.m2 <- update(mudas.m1,.~.+substrato)
##Os coeficientes
cf.mud2 <- coef(mudas.m2)
cf.mud2
## A matriz do modelo
head(model.matrix(mudas.m2))
##resumo do modelo
summary(mudas.m2)
anova(mudas.m2)

##MORAL DA HISTORIA: uma anova é o teste de siginificancia marginal para um modelo linear com variaveis categoricas!