############################## ####7b Regressão Multipla ############################## ##1. Uma estimativa da incerteza na previsao do modelo bebes <- read.table("babies.csv", header= TRUE, as.is = TRUE, sep= "\t") #carregando os dados str(bebes) #checando a estrutura dos dados lm.beb <- lm(bwt ~ gestation, data=bebes) #construindo o modelo ##Construindo o modelo #Definindo os parametros gerais x11() par(mar=c(6,6,2,2)) par(pch= 16, tcl= 0.3) par(las= 1, bty= "n", bg= "gray70") #Executando a plotagem vazia plot(NA, cex= 1.4, axes = FALSE, ann= FALSE, xlim= range(bebes$gestation), ylim= range(bebes$bwt)) rect(par()$usr[1],par()$usr[3], par()$usr[2],par()$usr[4], col="gray90") mtext("Peso dos bebês (oz)", side=2, cex= 1.7, line= 4, family="serif", las= 0) mtext("Tempo de gestação (dias)", side=1, cex= 1.7, line= 4, family="serif") #Adicionando as observacoes e eixos points(bwt ~ gestation, data = bebes, pch = 19, cex= 1.5, col= rgb(0,1,0, 0.4)) axis(1) axis(2) #Calculando CI 95% x.nov <- seq(min(bebes$gestation)-10, max(bebes$gestation)+10, length.out=100) ss <- sum((bebes$gestation - mean(bebes$gestation))^2) s2 <- var(bebes$bwt) n <- length(bebes$gestation) s.nov <- (x.nov - mean(bebes$gestation))^2 se.y <- sqrt(s2*((1/n) + (s.nov)/ss)) t <- qt(0.025, length(bebes$gestation)-2) se.y.t <- abs(t)*se.y fit <- coef(lm.beb)[[1]] + coef(lm.beb)[[2]]*x.nov predito <- data.frame(fit= fit, lwr= fit - se.y.t, upr= fit + se.y.t) #Adicionando o CI 95% ao grafico polygon(c(rev(x.nov), x.nov), c(rev(predito[ ,3]), predito[ ,2]), col = rgb(0,0,0,0.2), border = NA) abline(lm.beb, lwd = 2.5, col = "black") ############################## ##2. Galileu estava Certo? ############################## #Partindo do tutorial Ajuste de Polinomios, avalie se um polinômio de terceiro grau é um melhor modelo para descrever os dados do experimento de Galileu. A equação para este modelo é: #As alturas de lançamentos e as distâncias percorridas obtidas por Galileu sao: init.h = c(600, 700, 800, 950, 1100, 1300, 1500) h.d = c(253, 337, 395, 451, 495, 534, 573) mod1 <- lm(h.d ~ 1) #Melhor modelo do tutorial mod2 <- lm(h.d ~ init.h + I(init.h^2)) #Modelo do exercicio (polinomio de terceiro grau) mod3 <- lm(h.d ~ init.h + I(init.h^2) + I(init.h^3)) #Comparando os dois modelos anova(mod2, mod3) #RESP: O polinomio de terceiro grau (mod3) e de fato o melhor modelo. A adicao do novo parametro (init.h^3) #explicou muita variacao da variavel resposta, ja que a quantidade de variacao nao explicada pelo #mod3 (Residual Sum of Squares, RSS= 48.25) e muito menor que a variacao nao explicada pelo mod2 (polinomio #de segundo grau; RSS= 744.08). Logo, o teste estatistico F sobre essa comparacao de variancia particionada e #significativo ao nivel de 5%, ja que a probabilidade de estarmos errados com relacao a conclusao aqui explosta #e de 0.715%. ############################## ##3. Massa de Recem-Nascidos ############################# #Utilize todas as variáveis disponíveis no conjunto de dados Massa ao Nascer e Dados da mãe para chegar ao melhor modelo de previsão da massa de bebês ao nascer. #Parta do arquivo acima e não do que foi usado no tutorial. Antes de começar, elimine as observações com dados faltantes, veja descrição do arquivo. #Carregando os dados babies <- read.table("babies.txt", header= TRUE, as.is = TRUE, sep= " ") head(babies) #checando se a leitura esta correta dim(babies) #checando o numero de linhas e colunas do arquivo de entrada #Checando a presenca de observacoes com dados faltantes babies[which(babies$bwt==999),] #nenhum dado faltante babies[which(babies$gestation==999),] babies[which(babies$parity==9),] #nenhum dado faltante babies[which(babies$age==99),] babies[which(babies$height==99),] babies[which(babies$weight==999),] babies[which(babies$smoke==9),] #Removendo observacoes com dados faltantes babies <- babies[-which(babies$gestation==999|babies$age==99|babies$height==99|babies$weight==999| babies$smoke==9),] dim(babies) #62 linhas (observacoes) removidas #Analisando a correlacao entre as preditoras round(cor(babies[,c(2,4,5,6)], method="pearson"), 2) #RESULT: para evitar qualquer problema com colinearidade, a variavel altura das gestante (height) #foi removida antes de criarmos os modelos ja que sua correlacao com o peso (weight) e razoavel (r= 0.44). #Ainda que a correlacao possa ser considerada moderada, optei por evitar essa associacao entre preditoras. ##Criando o modelo full e analisando as estimativas #Assim como a variavel height, a variavel parity foi removida antes de criar o modelo full pois nao acredito, #a priori, que a primeira gestacao afete o peso do recem nascido diferente das demais gestacoes m.full <- lm(bwt ~ gestation * age * weight * smoke, data= babies) summary(m.full) #RESULT: como nao obtivemos nenhuma estimativa significativa, vamos iniciar a reducao do modelo #removendo as interacoes de maior grau. ##Removendo a interacao de quarta ordem e avaliando sua importancia m.01 <- update(m.full, .~. - gestation:age:weight:smoke, data= babies) anova(m.01, m.full) summary(m.01) #RESULT: esta interacao nao explicou uma quantidade de variacao suficiente para ser mantida (Pr(>F)= 0.655). ##Analisando e avaliando as interacoes de terceira ordem m.02 <- update(m.01, .~. - gestation:age:weight, data= babies) anova(m.01, m.02) m.03 <- update(m.01, .~. - gestation:age:smoke, data= babies) anova(m.01, m.03) m.04 <- update(m.01, .~. - gestation:weight:smoke, data= babies) anova(m.01, m.04) m.05 <- update(m.01, .~. - age:weight:smoke, data= babies) anova(m.01, m.05) #RESULT: novamente, nenhuma das interacoes explicou uma quantidade de variacao suficiente para serem mantidas. #Modelo atual sem as interacoes de terceira ordem m.06 <- update(m.01, .~. - (gestation:age:weight + gestation:age:smoke + gestation:weight:smoke + age:weight:smoke), data= babies) summary(m.06) #Como ainda temos muitos parametros estimados nao significativos, vamos seguir reduzindo o modelo. ##Analisando e avaliando as interacoes de segunda ordem m.07 <- update(m.06, .~. - gestation:age, data= babies) #interacao deve ser mantida anova(m.06, m.07) m.08 <- update(m.06, .~. - gestation:weight, data= babies) anova(m.06, m.08) m.09 <- update(m.06, .~. - age:weight, data= babies) anova(m.06, m.09) m.10 <- update(m.06, .~. - gestation:smoke, data= babies) #interacao deve ser mantida anova(m.06, m.10) m.11 <- update(m.06, .~. - age:smoke, data= babies) anova(m.06, m.11) m.12 <- update(m.06, .~. - weight:smoke, data= babies) anova(m.06, m.12) #RESULT: apenas as interacoes gestation:age e gestation:smoke explicaram uma razoavel quantidade de variacao, #sendo apropriado, portanto, mante-las. ##Logo, o modelo que provavelmente melhor preve a massa de bebes ao nascer e o modelo: #bwt ~ gestation + age + weight + smoke + gestation:age + gestation:smoke. m.13 <- update(m.06, .~. - (gestation:weight + age:weight + age:smoke + weight:smoke), data= babies) summary(m.13)