######################################## # # aula de 18 de Maio 09 # ######################################## #################################### bp <- read.table("lowbwtdata.dat", header = T) dim(bp) names(bp) <- tolower(names(bp)) head(bp) #variáveis categóricas bp$race <- factor(bp$race) bp$smoke <- factor(bp$smoke) bp$ht <- factor(bp$ht) bp$ui <- factor(bp$ui) bp$low <- factor(bp$low) bp$maekg<-bp$lwt*0.453 bp$bwt<-bp$bwt/1000 summary(bp) #Análise exploratória #Sumário dos dados as.data.frame(table(bp$ptl)) as.data.frame(table(bp$ftv)) #Distribuição do peso e de algumas covariavéis hist(bp$bwt,main="peso à nascença",xlab="btw (Kg)",ylab="nº de casos") plot(density(bp$bwt),main="peso à nascença") barplot(table(bp$ptl),main="PTL - historial de partos prematuros", xlab="PTL",ylab="nº de casos") barplot(table(bp$ftv),main="FTV - número de visitas médicas", xlab="FTV",ylab="nº de casos") # Observe os gráficos de associação entre o peso ao nascer e as covariáveis plot(bp$age, bp$bwt, pch=19,main="peso à nascença vs idade da mãe", xlab="age",ylab="bwt") plot(bp$maekg, bp$bwt, pch=19,main="peso à nascença vs peso da mãe",xlab="maekg",ylab="bwt") plot(bp$ptl, bp$bwt, pch=19,main="peso à nascença vs nº partos prematuros",xlab="PTL",ylab="bwt") plot(bp$ftv, bp$bwt, pch=19,main="peso à nascença vs nº visitas médicas",xlab="FTV",ylab="bwt") boxplot(bwt~race,data=bp,main="distribuição do peso à nascença por raça",xlab="raça",ylab="bwt" ) boxplot(bwt~smoke,data=bp,main="distribuição do peso à nascença por hábito de fumo",xlab="smoke",ylab="bwt" ) boxplot(bwt~ht,data=bp,main="distribuição do peso à nascença por história de hipertensão",xlab="ht",ylab="bwt") boxplot(bwt~ui,data=bp,main="distribuição do peso à nascença por irritabilidade uterina",xlab="ui",ylab="bwt") # O que podemos dizer desta análise gráfica? #Observando os primeiros dois gráficos, não parece evidente a existência #de qualquer relação entre bwt e age ou maekg. #A representação de bwt vs ptl ou ftv também não revela a existência #de associação de tipo linear. Como estas duas variáveis são de tipo #discreto e assumem uma diversidade muito pequena de valores, poderemos #pensar em categorizar as variáveis fzendo alguns agrupamentos nos seus #valores. #Analisando os box-plot: #bwt por race: As medianas das raças negra e outra são aproximadas, #sendo a correspondente à raça branca a mais elevadas das três. #A variabilidade dentro dos grupos parece ser algo diferene, sendo maior #na raça branca. Note-se, no entanto, que os grupos têm dimesões bastante #distintas: table(bp$race) #Para saber quais as classes referência, temos contrasts(bp$race) contrasts(bp$smoke) contrasts(bp$ht) contrasts(bp$ui) contrasts(bp$low) #modelo sugerido: bp.lm1 <- lm(bwt ~age+maekg+race+ftv, data = bp) summary(bp.lm1) names(bp.lm1) # anova para comparCompare os dois modelos bp.lm1 e bp.lm2. #Qual a sua conclusão? anova(bp.lm2, bp.lm1) bp.nulo<-lm(bwt ~1,data = bp) summary(bp.nulo) anova(bp.nulo, bp.lm1) # teste de wald para verificar quais podem sair summary(bp.lm1) # Estime um modelo sem as variáveis age e ftv bp.lm2 <- lm(bwt ~maekg+race, data = bp) summary(bp.lm2) #Compare os dois modelos. Qual a sua conclusão? anova(bp.lm2, bp.lm1) #Estime um modelo que tenha apenas a variável smoke como #variável independente. bp.lm3 <- lm(bwt ~ smoke, data = bp) summary(bp.lm3) #Utilize o modelo para calcular o peso médio ao nascer das crianças #para as mães que fumaram durante a gravidez e para as que não fumaram #Compare o resultado que encontrou na alínea anterior com o resultado #calculado empiricamente tapply(bp$bwt,list(bp$smoke),mean,na.rm=T) #Qual a vantagem de fazer um modelo de regressão linear em relação ao #cálculo empírico? #Resp.: É a facilidade de se obter um intervalo de confiança, #especialmente quando temos muitas covariáveis. names(bp.lm3) #intervalo de confiança assintotico li<-bp.lm3$coefficients[2]-1.96*0.10697 ls<-bp.lm3$coefficients[2]+1.96*0.10697 li ls #intervalo de confiança exacto li<-bp.lm3$coefficients[2]-qt(0.975,187)*0.10697 ls<-bp.lm3$coefficients[2]+qt(0.975,187)*0.10697 li ls # Agora suponha que temos um modelo com uma varíavel categórica com #mais de duas categorias. Neste caso todas as categorias estarão a ser #comparadas com a categoria referência. Para exemplificar, vamos usar a #covariável raça em que a categoria referência é a raça branca. bp.lm4 <- lm(bwt ~ race, data = bp) summary(bp.lm4) #Conferindo qual é a categoria referência. contrasts(bp$race) #Compare com o resultado anterior com o cálculo empírico tapply(bp$bwt,list(bp$race),mean,na.rm=T) #Admita agora o caso em que a covariável é contínua e possui uma relação #linear com a variável resposta. Estime o modelo que tem como variável #regressora apenas o peso da mãe aquando da última menstruação: maekg bp.lm5<-lm(bwt~maekg,data=bp) summary(bp.lm5) #Escreva a equação do modelo estimado #Calcule o valor médio do peso ao nascer para uma criança cuja mãe #apresenta maekg igual ao peso médio #Diga se existe diferença significativa entre o peso esperado à nascença para duas crianças cujas mães tenham valores de \code{maekg} com uma diferença igual a 5 Kg. Formalize a questão. # Considere agora um modelo envolvendo uma variável contínua e uma #variável categórica maekg e race: bp.lm6<-lm(bwt~maekg+race,data=bp) summary(bp.lm6) anova(bp.lm5,bp.lm6) #Compare o coeficiente da variável maekg deste modelo com o do #modelo contendo apenas a variável maekg (modelo anterior). #Que significado tem esse coeficiente em cada um dos casos? #Escreva a equação do modelo para cada uma das raças. #Ajuste modelos simples (só com a variável maekg) considerando #os subconjuntos dos dados que se obtêm separando-os por race. #Compare estes modelos com o que escreveu na alínea anterior. #Qual a melhor abordagem? Porquê? bp.lm6b<-lm(bwt~maekg,data=bp,subset=(race=="1")) summary(bp.lm6b) bp.lm6n<-lm(bwt~maekg,data=bp,subset=(race=="2")) summary(bp.lm6n) bp.lm6o<-lm(bwt~maekg,data=bp,subset=(race=="3")) summary(bp.lm6o) plot(bp$lwt[which(bp$race=="1")],bp$bwt[which(bp$race=="1")],col=1,xlab="lwt",ylab="bwt") points(bp$lwt[which(bp$race=="2")],bp$bwt[which(bp$race=="2")],col=2) points(bp$lwt[which(bp$race=="3")],bp$bwt[which(bp$race=="3")],col=3) abline(bp.lm6b,col=1) abline(bp.lm6n,col=2) abline(bp.lm6o,col=3) #(Selecção de modelos) #Admita que não tem ideia, à partida, sobre quais as variáveis que deve #incluir no modelo. Comecemos por construir o modelo contendo todas as #variáveis disponíveis com excepção de \code{id}, \code{low} e \code{lwt} bp.modcomp<-lm(bwt~.-id-low-lwt,data=bp) summary(bp.modcomp) #Utilize a função \code{step} para seleccionar o melhor modelo. #Analise os resultados parciais. step2<-step(bp.modcomp,direction="both") summary(step2) anova(bp.modcomp,bp.nulo) bp.mod7<-lm(bwt~.-id-low-lwt-ui,data=bp) anova(bp.modcomp,bp.mod7) #Recorde-se que as variáveis \code{ftv} e \code{ptl} são variáveis #quantitativas. No entanto, são de natureza discreta apresentando uma #diversidade de valores muito pequena. A questão que deve colocar-se é #se não seriam mais úteis se fossem categorizadas. Categorize as #variáveis e repita o procedimento de selecção do melhor modelo #utilizando o stepwise. Compare com os resultados anteriores. table(bp$ftv) table(bp$ptl) #Podemos juntar criar três categorias para ftv: #``nenhuma consulta",``uma consulta", ``duas ou mais consultas". bp$ftvcat<-0 bp$ftvcat[which(bp$ftv==1)]<-1 bp$ftvcat[which(bp$ftv>1)]<-2 bp$ftvcat<-factor(bp$ftvcat) table(bp$ftvcat) # Para a variável ptl criamos apenas duas categorias: # ``não"=0 e ``sim"=1. bp$ptlcat<-0 bp$ptlcat[which(bp$ptl>0)]<-1 bp$ptlcat<-factor(bp$ptlcat) table(bp$ptlcat) bp$ptlcat <- factor(bp$ptlcat) bp$ftvcat <- factor(bp$ftvcat) #E vamos repetir o stepwise. bp.modcomp2<-lm(bwt~.-id-low-lwt-ptl-ftv,data=bp) summary(bp.modcomp2) step3<-step(bp.modcomp2,direction="both") summary(step3) #Concluímos que, mesmo categorizadas, estas variáveis não são #importantes para modelar bwt # diagnóstico do modelo final #Análise de Resíduos # residuos padronizados bp.final<-step3 summary(bp.final) names(bp.final) #histograma dos resíduos hist(bp.final$residuals) bp.final$stdres<-bp.final$residuals/sd(bp.final$residuals) hist(bp.final$stdres) # verificação do pressuposto de normalidade: <<>>= library(car) qq.plot(bp.final$stdres,pch=19, main="Normal Q-Q plot", ylab="quantis empíricos", xlab="quantis teóricos da normal(0,1)") abline(0,1,lwd=2) #Verificação do pressuposto de homocedasticidade plot(bp.final$fitted.values,bp.final$stdres, main="resíduos padronizados vs valores ajustados", ylab="resíduos padronizados", xlab="valores ajustados",pch=21,bg=2,col=2) abline(h=0,lty=2) abline(h=-2,lty=3) abline(h=2,lty=3) #Verificação do pressuposto de independência plot(bp.final$fitted.values, main="valores ajustados vs index", ylab="valores ajustados",pch=21,bg=2,col=2) #Adequabilidade das variáveis regressoras #peso da mãe plot(bp$maekg,bp.final$stdres, main="resíduos padronizados vs peso da mãe", xlab="peso da mãe (maekg)",ylab="resíduos padronizados", pch=21,bg=2,col=2) abline(h=0,lty=2) hist(bp$maekg) #race boxplot(bp.final$stdres[which(bp$race==1)], bp.final$stdres[which(bp$race==2)], bp.final$stdres[which(bp$race==3)], col=c("green","blue","red"),names=c("raça branca","raça negra","raça=outra"), ylab="resíduos padronizados", main="distribuição dos resíduos padronizados por raça") plot(bp.final$fitted.values[which(bp$race==1)], main="valores ajustados vs index por raça", ylab="valores ajustados",pch="b",bg=2,col=2) points(bp.final$fitted.values[which(bp$race==2)],pch="n",bg=4,col=4) points(bp.final$fitted.values[which(bp$race==3)],pch="o",bg=3,col=3) abline(h=0,lty=2) #hábito de fumo boxplot(bp.final$stdres[which(bp$smoke==0)], bp.final$stdres[which(bp$smoke==1)], col=c("green","red"),names=c("não fuma","fuma"), ylab="resíduos padronizados", main="distribuição dos resíduos padronizados por hábito de fumo") plot(bp.final$fitted.values[which(bp$smoke==0)], main="valores ajustados vs index por hábito de fumo", ylab="valores ajustados",pch="n",bg=3,col=3) points(bp.final$fitted.values[which(bp$smoke==1)],pch="s",bg=2,col=2) abline(h=0,lty=2) #presença de hipertensão ht boxplot(bp.final$stdres[which(bp$ht==0)], bp.final$stdres[which(bp$ht==1)], col=c("green","red"),names=c("não hipertensa","hipertensa"), ylab="resíduos padronizados", main="distribuição dos resíduos padronizados por hipertensão") plot(bp.final$fitted.values[which(bp$ht==0)], main="valores ajustados vs index por hipertensão", ylab="valores ajustados",pch="n",bg=3,col=3) points(bp.final$fitted.values[which(bp$ht==1)],pch="s",bg=2,col=2) abline(h=0,lty=2) #O número de mulheres hipertensas é muito pequeno: table(bp$ht) #presença de irritabilidade uterina ui boxplot(bp.final$stdres[which(bp$ui==0)], bp.final$stdres[which(bp$ui==1)], col=c("green","red"),names=c("sem irritabilidade uterina","com irritabilidade uterina"), ylab="resíduos padronizados", main="distribuição dos resíduos padronizados por irritabilidade uterina") plot(bp.final$fitted.values[which(bp$ui==0)], main="valores ajustados vs index por irritabilidade uterina", ylab="valores ajustados",pch="n",bg=3,col=3) points(bp.final$fitted.values[which(bp$ui==1)],pch="s",bg=2,col=2) abline(h=0,lty=2) #O número de mulheres com irritabilidade uterina também é muito pequeno: table(bp$ui) #(Previsão)Considere um caso em que # maekg}=76.5, race=3,smoke=1,ht=1,ui=0, #Qual o peso esperado à nascença para uma criança cuja mãe apresenta #as características acima? #Obtenha um intervalo de predição a um nível de 95% de confiança para #a previsão. #Obtenha um intervalo de confiança a um nível de 95% de confiança para #o peso médio à nascença das crianças cujas mães possuem aquelas #características. predict(bp.final) new<-c(76.5,3,1,1,0) new<-data.frame(maekg=76.5,race="3",smoke="1",ht="1",ui="0") predict(bp.final, new,se.fit = TRUE) pred.w.plim <- predict(bp.final,new, interval="prediction") pred.w.plim pred.w.clim <- predict(bp.final,new, interval="confidence") pred.w.clim