#Exercícios 7 - Regressões Lineares Simples e Múltiplas
##Gustavo Agudelo. No. USP: 8893871

setwd("E:/Science/USP/Mestrado Fisiologia Geral_IB-USP/Disciplinas/Primeiro semestre_2014-1/Uso da Linguagem R para Análise de Dados em Ecologia_2014-1/Exercícios R/Exercícios 7 - Regressões Lineares Simples e Múltiplas")

#Parte A: Regressões Lineares Simples

#7.1 Altura na Infância e Na Vida Adulta

#Parte 1: regressão
alturas <- data.frame("a.dois.anos"=c(39, 30, 32, 34, 35, 36, 36, 30), "a.adulto"=c(71, 63, 63, 67, 68, 68, 70, 64))
alturas

alturas.lm <- lm(a.adulto ~ a.dois.anos, data=alturas)
summary(alturas.lm)
plot(alturas.lm)
anova(alturas.lm)
#Há uma relação significativa entre a altura aos dois anos de idade e a altura na idade adulta

alturas.conf <- confint(alturas.lm)
alturas.conf

#Parte 2: visualização

alturas$a.adulto.esp = (alturas$a.dois.anos*2)
alturas.esp.lm <- lm(a.adulto.esp ~ a.dois.anos, data=alturas)

plot(a.adulto ~ a.dois.anos, data=alturas)
abline(alturas.lm, col="blue", lwd=2)
abline(alturas.esp.lm, col="red", lwd=1)

anova(alturas.lm, alturas.esp.lm)

#Não, não os corroboram, porque os dois modelos são bem diferentes (p<0.01).


#7.2 O modelo mais simples possível

library(MASS)
data(Animals)
str(Animals)
summary(Animals)
anim.m2 <- lm(log(brain)~log(body),data=Animals, subset=!(log(Animals$body)>8&log(Animals$brain)<6))
anim.m0 <- lm(log(brain)~1, data=Animals, subset=!(log(Animals$body)>8&log(Animals$brain)<6))

anova(anim.m0,anim.m2)

#Respostas.

anova(anim.m2)

#1. A relação do comando anterior com aquele que compara os modelos anim.m0 e anim.m2 é que quando a função anova() é usada com o modelo anim.m2 aparece a tabela de análise de variância para este modelo.
# Neste caso, os resultados são muito parecidos porque a comparação dos modelos está dizendo que efetivamente o modelo anim.m2 é muito melhor para explicar a variação dos dados do que o modelo nulo (anim.m0).

summary(anim.m0)
mean(log(Animals$brain[!(log(Animals$body)>8&log(Animals$brain)<6)]))
sd(log(Animals$brain[!(log(Animals$body)>8&log(Animals$brain)<6)]))

#2. O valor da média obtido a partir do comando acima corresponde ao valor do coeficiente do intercepto
# Igualmente, o valor do desvio padrão corresponde ao erro associado aos residuais.


#7.3 Resíduos de Iris

iris.setosa <- subset(iris, iris$Species=="setosa")
lm.iris <- lm(iris.setosa$Sepal.Width ~ iris.setosa$Sepal.Length, data=iris.setosa)

lm.iris.coef <- coef(lm.iris)
lm.iris.coef

lm2.iris <- lm(iris.setosa$Sepal.Length ~ iris.setosa$Petal.Length, data=iris.setosa)
lm3.iris <- lm(iris.setosa$Sepal.Width ~ iris.setosa$Petal.Length, data=iris.setosa)

lm2.residuals <- as.numeric(residuals(lm2.iris))
lm3.residuals <- as.numeric(residuals(lm3.iris))

lm.iris.nopetal <- lm(lm3.residuals ~ lm2.residuals)

lm.iris.nopetal.coef <- coef(lm.iris.nopetal)
lm.iris.nopetal.coef 

summary.lm.iris <- summary(lm.iris)
summary.lm.iris.nopetal <- summary(lm.iris.nopetal)

r2.lm.iris = summary.lm.iris$r.squared #55.14% da variância é explicada por este modelo
r2.lm.iris.nopetal = summary.lm.iris.nopetal$r.squared #53.72% da variância é explicada por este modelo


#7.4 Pressão e temperatura

pressure

p <- pressure$pressure
t <- pressure$temperature

#A relação entre as variáveis não parece ser linear, mas uma parábola.

reg1 <- lm(p~t)
summary(reg1)

reg2 <- update(reg1,.~. +I(t^2))
summary(reg2)

reg3 <- update(reg2,.~. +I(t^3))
summary(reg3)

summary.reg3 <- summary(reg3)

r2 = summary.reg3$r.squared
r2

#Dependendo da natureza da relação entre as variáveis, aumentar a ordem da regressão pode se ajustar melhor o não. Neste caso, o modelo de terceira ordem foi o que explicou a maioria da variância de p em t. 

#7.5 Seriemas e Carcarás

aves <- read.table(file="http://ecologia.ib.usp.br/bie5782/lib/exe/fetch.php?media=dados:aves_cerrado.csv", header=T, row.names=1, sep=";", as.is=T)
str(aves)
head(aves)
is.na(aves[,2:4]) -> avesNAs
colSums(avesNAs)
aves[is.na(aves$urubu) | is.na(aves$carcara)| is.na(aves$seriema), ]
aves$urubu[is.na(aves$urubu)] <- 0
aves$carcara[is.na(aves$carcara)] <- 0
aves$seriema[is.na(aves$seriema)] <- 0
aves[aves$urubu==0|aves$carcara==0|aves$seriema==0,]
table(aves$fisionomia)
aves$fisionomia[aves$fisionomia=="ce"] <- "Ce"
table(aves$fisionomia)
class(aves$fisionomia)
aves$fisionomia <- factor(aves$fisionomia, levels=c("CL","CC","Ce"))
str(aves)

aves.ce <- subset(aves, aves$fisionomia=="Ce")
aves.cc <- subset(aves, aves$fisionomia=="CC")
aves.cl <- subset(aves, aves$fisionomia=="CL")

mod.ce <- lm(aves.ce$seriema ~ aves.ce$carcara, data=aves.ce)
summary(mod.ce) 

p.ce <- summary(mod.ce)$coefficients[2,4] 
coef.ce <- coef(mod.ce)
coef.ce 

mod.cc <- lm(aves.cc$seriema ~ aves.cc$carcara, data=aves.cc)
summary(mod.cc) 

p.cc <- summary(mod.cc)$coefficients[2,4] 
coef.cc <- coef(mod.cc)
coef.cc

mod.cl <- lm(aves.cl$seriema ~ aves.cl$carcara, data=aves.cl)
summary(mod.cl)

p.cl <- summary(mod.cl)$coefficients[2,4] 
coef.cl <- coef(mod.cl)
coef.cl

#Parte B: Regressão Múltipla

#7.1 Galileu estava Certo?

init.h = c(600, 700, 800, 950, 1100, 1300, 1500)
h.d = c(253, 337, 395, 451, 495, 534, 573)

plot(h.d~init.h)
mod1 <- lm(h.d~init.h)
abline(mod1)
mod2 <- update(mod1,.~. +I(init.h^2))
cf.m2 <- coef(mod2)
curve(cf.m2[1]+cf.m2[2]*x+cf.m2[3]*x^2, add=T, lty=2)
mod3 <- update(mod2,.~. +I(init.h^3))
cf.m3 <- coef(mod3)
curve(cf.m3[1]+cf.m3[2]*x+cf.m3[3]*x^2+cf.m3[4]*x^3, add=T, lty=3)

summary.m1 <- summary(mod1)
summary.m2 <- summary(mod2)
summary.m3 <- summary(mod3)

r2.m1 <- summary.m1$r.squared
r2.m2 <- summary.m2$r.squared
r2.m3 <- summary.m3$r.squared

anova.models <- anova(mod1, mod2, mod3)
RSS.models <- anova.models$RSS

#O modelo do polinômio de terceiro grau é o melhor para descrever os dados do experimento de Galileu, porque ele diminui a somatória dos residuais quadráticos, ou seja que tem o melhor ajuste dos dados com o modelo.
#Além disso, o r2 do polinômio de terceiro grau é maior neste modelo do que nos outros, explicando uma maior percentagem da variância na variável de resposta

#7.2 Massa de Recém-Nascidos
babies <- read.table("babies.txt", header=T,na.strings=999)
str(babies)
summary(babies)
table(is.na(babies))
babies$age[babies$age==99] <- NA
babies$height[babies$height==99] <- NA
babies$smoke[babies$smoke==9] <- NA
babies$smoke <- as.logical(babies$smoke)
babies <- babies[apply(is.na(babies),1,sum)==0,]
summary(babies)
str(babies)

#O primeiro modelo será considerado como o mais complicado. Por efeitos de parcimônia na interpretação, só consideram-se interações de segunda ordem

babies.mod1 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:parity+gestation:age+gestation:height+gestation:weight+gestation:smoke+
parity:age+parity:height+parity:weight+parity:smoke+
age:height+age:weight+age:smoke+
height:weight+height:smoke+
weight:smoke, data=babies)

summary(babies.mod1)
anova(babies.mod1)

#Tirando weight:smoke

babies.mod2 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:parity+gestation:age+gestation:height+gestation:weight+gestation:smoke+
parity:age+parity:height+parity:weight+parity:smoke+
age:height+age:weight+age:smoke+
height:weight+height:smoke, data=babies)

summary(babies.mod2)
anova(babies.mod2)

#Tirando parity:smoke

babies.mod3 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:parity+gestation:age+gestation:height+gestation:weight+gestation:smoke+
parity:age+parity:height+parity:weight+
age:height+age:weight+age:smoke+
height:weight+height:smoke, data=babies)

summary(babies.mod3)
anova(babies.mod3)

#Tirando age:height

babies.mod4 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:parity+gestation:age+gestation:height+gestation:weight+gestation:smoke+
parity:age+parity:height+parity:weight+
age:weight+age:smoke+
height:weight+height:smoke, data=babies)

summary(babies.mod4)
anova(babies.mod4)

#Tirando parity:weight

babies.mod5 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:parity+gestation:age+gestation:height+gestation:weight+gestation:smoke+
parity:age+parity:height+
age:weight+age:smoke+
height:weight+height:smoke, data=babies)

summary(babies.mod5)
anova(babies.mod5)

#Tirando height:weight

babies.mod6 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:parity+gestation:age+gestation:height+gestation:weight+gestation:smoke+
parity:age+parity:height+
age:weight+age:smoke+
height:smoke, data=babies)

summary(babies.mod6)
anova(babies.mod6)

#Tirando age:weight

babies.mod7 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:parity+gestation:age+gestation:height+gestation:weight+gestation:smoke+
parity:age+parity:height+
age:smoke+
height:smoke, data=babies)

summary(babies.mod7)
anova(babies.mod7)

#Tirando gestation:parity

babies.mod8 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:age+gestation:height+gestation:weight+gestation:smoke+
parity:age+parity:height+
age:smoke+
height:smoke, data=babies)

summary(babies.mod8)
anova(babies.mod8)

#Tirando gestation:height

babies.mod9 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:age+gestation:weight+gestation:smoke+
parity:age+parity:height+
age:smoke+
height:smoke, data=babies)

summary(babies.mod9)
anova(babies.mod9)

#Tirando age:smoke

babies.mod10 <- lm(bwt~gestation+parity+age+height+weight+smoke+
gestation:age+gestation:weight+gestation:smoke+
parity:age+parity:height+
height:smoke, data=babies)

summary(babies.mod10)
anova(babies.mod10)

#Tirando age + suas interações

babies.mod11 <- lm(bwt~gestation+parity+height+weight+smoke+
gestation:weight+gestation:smoke+
parity:height+
height:smoke, data=babies)

summary(babies.mod11)
anova(babies.mod11)

#Tirando height:smoke

babies.mod12 <- lm(bwt~gestation+parity+height+weight+smoke+
gestation:weight+gestation:smoke+
parity:height, data=babies)

summary(babies.mod12)
anova(babies.mod12)

#Tirando gestation:weight

babies.mod13 <- lm(bwt~gestation+parity+height+weight+smoke+
gestation:smoke+
parity:height, data=babies)

summary(babies.mod13)
anova(babies.mod13)

#Tirando parity:height

babies.mod14 <- lm(bwt~gestation+parity+height+weight+smoke+
gestation:smoke, data=babies)

summary(babies.mod14)
anova(babies.mod14)

##Eu comecei pelo modelo mais complicado, extraindo aquelas variáveis ou interações que tiveram os valores de p mais grandes e não significativos. 
##Finalmente, o modelo que faz a melhor previsão da massa de bebês ao nascer foi o modelo 14 (babies.mod14).