##### 7b. Exercícios de Regressão Múltipla ##### ### 1. incerteza do modelo para a rela??o peso do beb? ### ## Construa esse gr?fico para avaliar a incerteza do modelo para a rela??o peso do # beb? (kg) ~ tempo de gesta??o (dias) dos dados Massa ao Nascer e Dados da m?e, usando a # f?rmula acima para estimar o intervalo de confian?a de 95% bebes <- read.table("babies.csv", header = TRUE, as.is = TRUE, sep = "\t") bebes <- bebes[bebes$bwt != "999" & bebes$gestation != "999" & bebes$parity != "9" & bebes$age != "99" & bebes$height != "99" & bebes$weight != "999" & bebes$smoke != "9", ] head(bebes) str(bebes) ## A. Modelo para peso do beb? (kg) em fun??o do tempo de gesta??o (dias) ml02 <- lm(bwt ~ gestation, data = bebes) summary(ml02) anova(ml02) coef(ml02) # Criando os valores preditos newx <- seq(min(bebes$gestation) - 10, max(bebes$gestation) + 10, length.out = 1174) ypredic <- (coef(ml02)[2] * newx) + coef(ml02)[1] ## B. Erro padrão da previsão do modelo # Com os valores observados da variavel x não deu certo #SSX <- sum((bebes$gestation - mean(bebes$gestation))^2) #Xx <- (bebes$gestation - mean(bebes$gestation))^2 #recip <- 1 / length(bebes$gestation) #s2 <- var(bebes$bwt) #n <- length(bebes$bwt) SSX <- sum((newx - mean(newx))^2) Xx <- (newx - mean(newx))^2 recip <- 1 / length(bebes$gestation) s2 <- var(bebes$bwt) n <- length(bebes$bwt) se <- sqrt(s2*(recip + (Xx / SSX))) # Multiplicar pelo valor t de Student, para n-2 graus de liberdade t.val <- qt(p = 0.975, df = (n - 2)) ## C. Gr?fico do intervalo de confian?a de 95% # Criando o intervalo de confiança # Acho que este intervao calculado a mão tem que converger com os valores obtidos pela função `predict` # porém eu não consegui. Procurei na web e achei muitas fórmulas diferentes para calcular o intervalo # de confiança, então não sei onde ficou o problema intervalo <- data.frame(ypredic, ylow = (ypredic - (t.val * se)), yup = (ypredic + (t.val * se))) head(intervalo) #preds <- predict(ml02, newdata = data.frame(gestation = newx), # interval = 'confidence') #head(preds) # Construindo o gr?fico com o intervalo de confian?a plot(bwt ~ gestation, data = bebes, xlab = "Tempo de gestação (dias)", ylab = "Peso do bebê (kg)", col = "darkblue") abline(ml02) polygon(c(rev(newx), newx), c(rev(intervalo[ ,3]), intervalo[ ,2]), col = rgb(0,0,0,0.2), border = NA) ### 2. Galileu estava Certo? ### ## Partindo do tutorial Ajuste de Polin?mios, 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 ? 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) # Ajustando modelo de primeiro grau mod1 <- lm(h.d~init.h) summary(mod1) ## Ajustando modelo de segundo grau mod2 <- lm(h.d ~ init.h + I(init.h^2)) summary(mod2) # O modelo de par?bola ? claramente superior, o que fica evidente se acrescentamos as # linhas dos previstos pelos dois modelos ao gr?fico plot(h.d ~ init.h) abline(mod1) cf.m2 <- coef(mod2) curve(cf.m2[1] + cf.m2[2] * x + cf.m2[3] * x^2, add = TRUE, lty = 2, col = "blue") ## Ajustando modelo do terceiro grau mod3 <- lm(h.d ~ init.h + I(init.h^2) + I(init.h^3)) summary(mod3) ## Comparando modelos # R/ O modelo de terceiro orden ? melhor porque explica uma porcao maior da variacao # (tem um RSS menor) dos daos. O anova resulta em differencas signficativas entre a # variacao explicada pelos dois modelos. Isso pode ser conferido no gr?fico anova(mod2, mod3) cf.m3 <- coef(mod3) curve(cf.m3[1] + cf.m3[2] * x + cf.m3[3] * x^2 + cf.m3[4] * x^3, add = TRUE, lty = 2, col = "red") ### 3. Massa de Rec?m-Nascidos ### # 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 bebes <- read.csv(file ="Planilha de dados Babies.csv", header = TRUE, as.is = TRUE, sep = " ") bebes <- bebes[bebes$bwt != "999" & bebes$gestation != "999" & bebes$parity != "9" & bebes$age != "99" & bebes$height != "99" & bebes$weight != "999" & bebes$smoke != "9", ] head(bebes) ## A. Modelo cheio # O resultado do summary indica que nenhuma das compara??es ? significativa, por?m o resultado do anova # indica que h? umas delas que podem explicar parte da varian?a, ent?o vou tentar incluir elas num modelo # mais simples do que o modelo cheio mfull <- lm(bwt ~ gestation * parity * age * height * weight * smoke, data = bebes) summary(mfull) a <- as.data.frame(summary(mfull)$coefficients) a[a$`Pr(>|t|)`<= 0.05, ] # Nenhum valor ? menor o igual a 0.05 anova(mfull) # Um outro jeito de comparar as variaveis str(anova(mfull)) anova(mfull)[anova(mfull)[5] <= 0.05, ] # Variaveis que explican parte da varia??o ## B. Modelo mais simples, removendo primeiro as interacoes m01 <- lm(bwt ~ gestation + parity + age + height + weight + smoke + gestation:age + gestation:weight + gestation:smoke, data = bebes) summary(m01) ## C. Comparando modelos anova(mfull, m01) # N?o h? diferen?a entre modelos, fico com o modelo mais simples ## D. Revomendo ainda mais variaveis que n?o explican a varia??o. Primeiro as intera?oes m02 <- lm(bwt ~ gestation + parity + age + height + weight + smoke + gestation:age + gestation:smoke, data = bebes) summary(m02) ## E. Comparando modelos anova(m01, m02) # N?o h? diferen?a entre modelos, fico com o modelo mais simples ## F. Revomendo ainda mais variaveis que n?o explican a varia??o. Agora a variavel gesta??o, no entanto # as intera??es com gesta??o s?o significativas, quais casos isso pode acontecer? m03 <- lm(bwt ~ parity + age + height + weight + smoke + gestation:age + gestation:smoke, data = bebes) summary(m03) ## G. Comparando modelos anova(m02, m03) # N?o h? diferen?a entre modelos, fico com o modelo mais simples ## H. Revomendo ainda mais variaveis que n?o explican a varia??o. Agora a variavel peso da m?e m04 <- lm(bwt ~ parity + age + height + smoke + gestation:age + gestation:smoke, data = bebes) summary(m04) ## I. Comparando modelos anova(m03, m04) # Parece que n?o h? diferen?a entre modelos, fico com o modelo mais simples ## J. Revomendo ainda mais variaveis que n?o explican a varia??o. Agora uma das intera??es por vez, # por?m elas s?o significativa m05 <- lm(bwt ~ parity + age + height + smoke + gestation:age, data = bebes) summary(m05) m06 <- lm(bwt ~ parity + age + height + smoke + gestation:smoke, data = bebes) summary(m06) ## K. Comparando modelos anova(m04, m05) # H? diferen?a entre modelos, ent?o o modelo com as intera??es ? melhor anova(m04, m06) # H? diferen?a entre modelos, ent?o o modelo com as intera??es ? melhor ## L. Interpretando o modelo m?nimo adecuado m04 <- lm(bwt ~ parity + age + height + smoke + gestation:age + gestation:smoke, data = bebes) summary(m04) anova(m04) coef(m04) ## M. Avalia??o do modelo # R/ Os residuais do modelo est?o bem ajustados par(mfrow = c(2,2)) plot(m04) par(mfrow = c(1,1))