#Aula 8 - Modelo Linear multiplo



#Exercicio 7.b Estimativa da incerteza na previsao do modelo

#Lendo a tabela e vendo sua estrutura
bebes=read.table("babies.txt", header=T, sep=" ", as.is = T)
head(bebes)
tail(bebes)
str(bebes)
summary(bebes)

#Identificando os valores 999
range(bebes$gestation)
range(bebes$bwt)

#Remocao de 999 da gestation
bebes2=bebes[bebes$gestation<999, ]
bebes2

#O peso dos animais foram anotados em "oncas" e o exercicio pede que o modelo seja feito com o peso do bebe em Kg. 1 libra = 0,45392kg. Assim, vou criar uma coluna chamada bebe peso = bebe.p
bebes2$bebe.p=bebes2$bwt * 0.0283495
head(bebes2)  #conferindo

#Criando um intervalo de 100 pontos no eixo x (tempo de gestacao). Que sera os valores de x utilizados na formula do erro padrao
min(bebes2$gestation)  #valor minimo de gestacao=148 dias
max(bebes2$gestation)  #valor maximo de gestacao=353dias
vetor.prev=seq(from = 148, to=353, length.out=100)  #vetor previsto
vetor.prev

#Calculo do erro padrao (se):
(med.gest=mean(bebes2$gestation))  #Tempo media de gestacao=279.3385 dias
dif.x=((vetor.prev - med.gest)^2)  #Calculo da diferenca de x: (x-xmed)^2
ssx=sum((bebes2$gestation-med.gest)^2)  #Soma dos desvios quadraticos
(n=length(bebes2$bwt))  #numero de amostras
s2=var(bebes2$bebe.p)   #variancia do peso dos bebes (kg)
se=sqrt(s2*(1/n + dif.x/ssx))  #Erro padrao

#Para estimar o intervalo de confianca (ic), precisa calcular o valor de t de Student
?qt
t=qt(0.025, df= n - 2, lower.tail = F)  #teste t
ic=se*t
ic

#Modelo linear peso do bebe com o tempo de gestacao
ml=lm(bebes2$bebe.p ~ bebes2$gestation)
summary(ml)
str(ml)
plot(bwt~gestation, data=bebes2)
abline(bebes2, col="blue")

coef.lm=coef(ml)  #coeficientes do modelo: intercept - valor do y, onde x=0 e o outro coeficiente e a inclinacao da reta
coef.lm

ic.max=(coef.lm[1]+(coef.lm[2])*vetor.prev)+ic  #reta do intervalo de confianca maximo (95+2.5=97.5%), baseado no coeficiente do modelo
ic.min=(coef.lm[1]+(coef.lm[2])*vetor.prev) - ic  #intervalo de confianca minima (2.5%)
ic=c(ic.max, ic.min)

#grafico
par(mfrow=c(1,1))
plot(bebes2$bebe.p ~ bebes2$gestation, xlim = c(min(vetor.prev), max(vetor.prev)))
abline(ml, col="blue")
lines(x=vetor.prev, y=ic.max, col="red")
lines(x=vetor.prev, y=ic.min, col="red")


################################################################

#Exercicio 7.b Galileu estava certo?
#avaliando os modelos

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)

##modelo linear:y=a+bx
mod1 = lm(h.d~init.h)  
anova(mod1)

#modelo y=a+bx+cx^2
mod2 = update(mod1,.~. +I(init.h^2))  
?update
anova(mod1,mod2)  #comparando os modelos
summary(mod2)

#Grafico modelo2
abline(mod1)
cf.m2 = coef(mod2)
curve(cf.m2[1]+cf.m2[2]*x+cf.m2[3]*x^2, add=T, lty=2)

#modelo y=a+bx+cx^2+dx^3
mod3 = update(mod2,.~. +I(init.h^3))
summary(mod3)
anova(mod2,mod3)

#Grafico modelo3
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=2)

#A analise dos modelos mostraram que o modelo da equacao de 3 grau (mod3) teve maior ajuste aos dados comparado com os demais (mod1 e mod2). O valor de R^2 foi mais alto e os valores de p foram menores no mod3.


#################################################################

#Exercicio 7.b Massa de recem-nascidos

#Lendo a tabela e vendo sua estrutura
bebes3=read.table("babies.txt", header=T, sep=" ", as.is = T)

#Identificacao de dados faltantes (999, 99 e 9)
range(bebes3$gestation)  #999
range(bebes3$bwt)        #ok
range(bebes3$parity)     #ok
range(bebes3$age)        #99
range(bebes3$height)     #99
range(bebes3$weight)     #999
range(bebes3$smoke)      #9

#Transformando os dados faltanets em NA
bebes3$gestation[bebes3$gestation==999] = NA
bebes3$age[bebes3$age==99] = NA
bebes3$height[bebes3$height==99] = NA
bebes3$weight[bebes3$weight==999] = NA
bebes3$smoke[bebes3$smoke==9] = NA

is.na(bebes3)
?na.omit
na.omit(bebes3)

summary(bebes3)

#Pressupostos e hipoteses:
#Fator gestacao - H0: O tempo de gestacao influencia no peso da crianca
#Fator "parity" - foi removido do modelo. Considerei o pressuposto de que o numero de gestacoes anteriores nao influencia na gestacao atual.
#Fator idade da mae - foi removido do modelo. Pressuposto: a idade nao influencia no peso da criancas
#Fator altura da mae - tambem foi removido do modelo.
#Fator peso da mae antes da gestacao - H0: a qualidade nutricional da mae antes da gestacao ira se manter durante a gravidez
#Fator cigarro - H0: mulheres fumantes sentem menos fome, consequentemente comerao menos, e seus bebes terao menor peso

mod0=lm(bwt~gestation+weight+smoke+gestation:weight+weight:smoke, data=bebes3)
summary(mod0)
anova(mod0)   #comparacao com o modelo nulo

#remocao da interacao weight:smoke, pois teve o maior valor de p

mod1=lm(bwt~gestation+weight+smoke+gestation:weight, data=bebes3)
anova(mod0,mod1)  #comparando os modelos, o p foi >0.05. Portanto a interacao que foi removida realmente nao exercia efeito sobre a variavel resposta.
summary(mod1)     #o efeito da gestacao e do cigarro ficou mais robusta neste modelo, com menores valores de p

#remocao da interacao gestation:weight
mod2=lm(bwt~gestation+weight+smoke, data=bebes3)
anova(mod1,mod2)   #nao teve diferenca entre os modelos, entao fico com o mod2 que e mais simples
summary(mod2)
plot(mod2)

#Portanto a melhor previsao do modelo de massas de bebes e o mod2.

############################# THE END ##################################