#-------------------------------------------------------------------------------------------------------- # EXERCÍCIO 2 - ANÁLISE DE REGRESSÃO LINEAR # BANCO PAS #-------------------------------------------------------------------------------------------------------- # Altere o diretório de trabalho: setwd("") # Importe o banco de dados "PAS" e dê o nome de "dados" # Utilize o Menu # File -> Import Dataset -> From Excel... dados names(dados) # nomes das variáveis dim(dados) # dimensão da matriz/tabela/banco head(dados) # metade da matriz/tabela/banco #-------------------------------------------------------------------------------------------------------- # ANÁLISE DESCRITIVA #-------------------------------------------------------------------------------------------------------- # Produzindo um resumo rápido do conteúdo do banco de dados... summary(dados) # valor mínimo, primeiro quartil, mediana, média, terceiro quartil e valor máximo das variáveis # ou de uma variável... summary(dados$PAS) # valor mínimo, primeiro quartil, mediana, média, terceiro quartil e valor máximo de uma variável (pressão arterial sistólica) summary(dados$QUET) # valor mínimo, primeiro quartil, mediana, média, terceiro quartil e valor máximo de uma variável (índice de Quetelet) summary(dados$AGE) # valor mínimo, primeiro quartil, mediana, média, terceiro quartil e valor máximo de uma variável (idade) summary(dados$SMK) # valor mínimo, primeiro quartil, mediana, média, terceiro quartil e valor máximo de uma variável (fumo sim ou não) # DESCREVENDO UMA VARIÁVEL DE CADA VEZ #-------------------------------------------------------------------------------------------------------- # Medidas descritivas #-------------------------------------------------------------------------------------------------------- mean(dados$PAS) # média pressão arterial sistólica mean(dados$QUET) # média índice de quetelet mean(dados$AGE) # média idade mean(dados$SMK) # média fumo sim ou não var(dados$PAS) # variância pressão arterial sistólica var(dados$QUET) # variância índice de quetelet var(dados$AGE) # variância idade var(dados$SMK) # variância fumo sim ou não sd(dados$PAS) # desvio padrão pressão arterial sistólica sd(dados$QUET) # desvio padrão índice de quetelet sd(dados$AGE) # desvio padrão idade sd(dados$SMK) # desvio padrão fumo sim ou não quantile(dados$PAS, probs = c(0, 0.25, 0.5, 0.75, 1)) #quartis pressão arterial sistólica quantile(dados$QUET, probs = c(0, 0.25, 0.5, 0.75, 1)) #quartis índice de quetelet quantile(dados$AGE, probs = c(0, 0.25, 0.5, 0.75, 1)) #quartis idade quantile(dados$SMK, probs = c(0, 0.25, 0.5, 0.75, 1)) #quartis fumo sim ou não # Intervalo de confiança #-------------------------------------------------------------------------------------------------------- install.packages("epiDisplay") library(epiDisplay) ci.numeric(dados$PAS, alpha=0.05) # intervalo de confiança pressão arterial sistólica ci.numeric(dados$QUET, alpha=0.05) # intervalo de confiança índice de quetelet ci.numeric(dados$AGE, alpha=0.05) # intervalo de confiança idade ci.numeric(dados$SMK, alpha=0.05) # intervalo de confiança fumo sim ou não # Boxplot #-------------------------------------------------------------------------------------------------------- # Utilizando a sintaxe: boxplot(dados$PAS, ylab="PAS (mmHg)", col="orange") # boxplot pressão arterial sistólica boxplot(dados$QUET, ylab="QUET (kg/m³)", col="orange") # boxplot índice de quetelet boxplot(dados$AGE, ylab="AGE (anos)", col="orange") # boxplot idade boxplot(dados$SMK, ylab="SMK (sim ou não)", col="orange") # boxplot fumo sim ou não # Veja as cores disponíveis no R: colors() # Histograma #-------------------------------------------------------------------------------------------------------- # Utilizando a sintaxe: hist(dados$PAS, scale="density", col="lightblue", xlab="PAS (mmHg)", ylab="Densidade de frequência") # histograma pressão arterial sistólica hist(dados$QUET, scale="density", col="lightblue", xlab="QUET (kg/m³)", ylab="Densidade de frequência") # histograma índice de quetelet hist(dados$AGE, scale="density", col="lightblue", xlab="AGE (anos)", ylab="Densidade de frequência") # histograma anos hist(dados$SMK, scale="density", col="lightblue", xlab="SMK (sim ou não)", ylab="Densidade de frequência") # histograma fumo sim ou não # DESCREVENDO A RELAÇÃO/ASSOCIAÇÃO ENTRE PAS E QUET/AGE/SMK #-------------------------------------------------------------------------------------------------------- # Diagrama de dispersão #-------------------------------------------------------------------------------------------------------- # Utilizando a sintaxe: plot(PAS~QUET, xlab="QUET (kg/m³)", ylab="PAS (mmHg)", col="blue", data=dados) # diagrama dispersão PAS em função de QUET plot(PAS~AGE, xlab="AGE (anos)", ylab="PAS (mmHg)", col="blue", data=dados) # diagrama dispersão PAS em função de AGE plot(PAS~SMK, xlab="SMK (sim ou não)", ylab="PAS (mmHg)", col="blue", data=dados) # diagrama dispersão PAS em função de SMK # Coeficicente de correlação linear de Pearson #-------------------------------------------------------------------------------------------------------- # Utilizando a sintaxe: cor(dados$PAS, dados$QUET) # valor do teste coeficiente correlação linear de pearson cor.test(dados$PAS, dados$QUET, alternative="two.sided", method="pearson") # teste coeficiente correlação linear de Pearson cor(dados$PAS, dados$AGE) # valor do teste coeficiente correlação linear de pearson cor.test(dados$PAS, dados$AGE, alternative="two.sided", method="pearson") # teste coeficiente correlação linear de pearson # Análise de variáveis quantitativas e qualitativas (PAS + SMK) #-------------------------------------------------------------------------------------------------------- class(dados$SMK) dados$fsmk <- factor(dados$SMK, levels=c(0,1), labels=c("Não", "Sim")) # criar fator para smk (0 = não; 1 = sim) dados$fsmk table(dados$fsmk) tapply(dados$PAS, dados$fsmk, mean) # média tapply(dados$PAS, dados$fsmk, var) # variância bartlett.test(PAS~fsmk, data=dados) install.packages("car") library(car) leveneTest(PAS ~ fsmk, data = dados) t.test(PAS~fsmk, alternative="two.sided", conf.level=0.95, var.equal=TRUE, data=dados) #-------------------------------------------------------------------------------------------------------- # AJUSTE DO MODELO DE REGRESSÃO LINEAR SIMPLES #-------------------------------------------------------------------------------------------------------- # Checando as suposições do modelo #-------------------------------------------------------------------------------------------------------- hist(dados$PAS, scale="density", col="lightblue", xlab="PAS (mmHg)", ylab="Densidade de frequência") # Gráfico para avaliar normalidade qqPlot(dados$PAS, dist="norm") # Gráfico para avaliar normalidade shapiro.test(dados$PAS) # Teste de normalidadende Shapiro-Wilk shapiro.test(dados$QUET) # Teste de normalidadende Shapiro-Wilk shapiro.test(dados$AGE) # Teste de normalidadende Shapiro-Wilk ks.test(dados$PAS, "pnorm", mean = mean(dados$PAS), sd = sd(dados$PAS)) # Teste de normalidade de Kolmogorov-Smirnov ks.test(dados$QUET, "pnorm", mean = mean(dados$QUET), sd = sd(dados$QUET)) # Teste de normalidade de Kolmogorov-Smirnov ks.test(dados$AGE, "pnorm", mean = mean(dados$AGE), sd = sd(dados$AGE)) # Teste de normalidade de Kolmogorov-Smirnov #-------------------------------------------------------------------------------------------------------- # MODELO DE REGRESSÃO LINEAR MÚLTIPLA #-------------------------------------------------------------------------------------------------------- # MATRIZ DE CORRELAÇÃO ENTRE AS VARIÁVEIS #-------------------------------------------------------------------------------------------------------- install.packages("polycor") library(polycor) hetcor(dplyr::select(data.frame(dados), PAS, QUET, AGE, SMK)) # Matriz de correlação #-------------------------------------------------------------------------------------------------------- scatterplotMatrix(~PAS + QUET + AGE + SMK, regLine=TRUE, smooth=FALSE, diagonal=list(method="none"), data=dados) #-------------------------------------------------------------------------------------------------------- # Ajustando o modelo #-------------------------------------------------------------------------------------------------------- RegModel.1 <- lm(PAS~QUET, data=dados) summary(RegModel.1) RegModel.2 <- lm(PAS~AGE, data=dados) summary(RegModel.2) RegModel.3 <- lm(PAS~SMK, data=dados) summary(RegModel.3) RegModel.4 <- lm(PAS~AGE+QUET, data=dados) summary(RegModel.4) RegModel.5 <- lm(PAS~AGE+QUET+SMK, data=dados) summary(RegModel.5) # Extraindo objetos do modelo #-------------------------------------------------------------------------------------------------------- objects(RegModel.1) # PAS + QUET objects(RegModel.2) # PAS + AGE objects(RegModel.3) # PAS + SMK objects(RegModel.4) # PAS + AGE + QUET objects(RegModel.5) # PAS + AGE + QUET + SMK coefficients(RegModel.1) # Coeficientes da reta de regressão ajustada PAS + QUET coef(RegModel.1) # Coeficientes da reta de regressão ajustada PAS + QUET coefficients(RegModel.2) # Coeficientes da reta de regressão ajustada PAS + AGE coef(RegModel.2) # Coeficientes da reta de regressão ajustada PAS + AGE coefficients(RegModel.3) # Coeficientes da reta de regressão ajustada PAS + SMK coef(RegModel.3) # Coeficientes da reta de regressão ajustada PAS + SMK coefficients(RegModel.4) # Coeficientes da reta de regressão ajustada PAS + AGE + QUET coef(RegModel.4) # Coeficientes da reta de regressão ajustada PAS + AGE + QUET coefficients(RegModel.5) # Coeficientes da reta de regressão ajustada PAS + AGE + QUET + SMK coef(RegModel.5) # Coeficientes da reta de regressão ajustada PAS + AGE + QUET + SMK anova(RegModel.1) # Tabela de Análise de Variância PAS + QUET anova(RegModel.2) # Tabela de Análise de Variância PAS + AGE anova(RegModel.3) # Tabela de Análise de Variância PAS + SMK anova(RegModel.4) # Tabela de Análise de Variância PAS + AGE + QUET anova(RegModel.5) # Tabela de Análise de Variância PAS + AGE + QUET + SMK fitted(RegModel.1) # Valores ajustados PAS + QUET fitted(RegModel.2) # Valores ajustados PAS + AGE fitted(RegModel.3) # Valores ajustados PAS + SMK fitted(RegModel.4) # Valores ajustados PAS + AGE + QUET fitted(RegModel.5) # Valores ajustados PAS + AGE + QUET + SMK residuals(RegModel.1) # Resíduos PAS + QUET residuals(RegModel.2) # Resíduos PAS + AGE residuals(RegModel.3) # Resíduos PAS + SMK residuals(RegModel.4) # Resíduos PAS + AGE + QUET residuals(RegModel.5) # Resíduos PAS + AGE + QUET + SMK rstandard(RegModel.1) # Resíduos studentizados PAS + QUET rstandard(RegModel.2) # Resíduos studentizados PAS + AGE rstandard(RegModel.3) # Resíduos studentizados PAS + SMK rstandard(RegModel.4) # Resíduos studentizados PAS + AGE + QUET rstandard(RegModel.5) # Resíduos studentizados PAS + AGE + QUET + SMK # Obtendo o coeficiente de determinação (R2) #-------------------------------------------------------------------------------------------------------- anova(RegModel.1) # PAS + QUET sqreg <- anova(RegModel.1)[1,2] sqr <- anova(RegModel.1)[2,2] sqt <- sqreg + sqr sqreg; sqr; sqt R2 <- sqreg / sqt R2 sqrt(R2) anova(RegModel.2) # PAS + AGE sqreg <- anova(RegModel.2)[1,2] sqr <- anova(RegModel.2)[2,2] sqt <- sqreg + sqr sqreg; sqr; sqt R2 <- sqreg / sqt R2 sqrt(R2) anova(RegModel.3) # PAS + SMK sqreg <- anova(RegModel.3)[1,2] sqr <- anova(RegModel.3)[2,2] sqt <- sqreg + sqr sqreg; sqr; sqt R2 <- sqreg / sqt R2 sqrt(R2) anova(RegModel.4) # PAS + AGE + QUET sqreg <- anova(RegModel.4)[1,2] + anova(RegModel.4)[2,2] sqr <- anova(RegModel.4)[3,2] sqt <- sqreg + sqr sqreg; sqr; sqt R2 <- sqreg / sqt R2 sqrt(R2) anova(RegModel.5) # PAS + AGE + QUET + SMK sqreg <- anova(RegModel.5)[1,2] + anova(RegModel.4)[2,2] + anova(RegModel.5)[3,2] sqr <- anova(RegModel.5)[4,2] sqt <- sqreg + sqr sqreg; sqr; sqt R2 <- sqreg / sqt R2 sqrt(R2) # Inserindo a reta ajustada no diagrama de dispersão #-------------------------------------------------------------------------------------------------------- plot(PAS~QUET, xlab="QUET (kg/m³)", ylab="PAS (mmHg)", col="blue", data=dados) abline(RegModel.1, col="red") #reta ajustada no diagrama de dispersão PAS e QUET plot(PAS~AGE, xlab="AGE (anos)", ylab="PAS (mmHg)", col="blue", data=dados) abline(RegModel.2, col="red") #reta ajustada no diagrama de dispersão PAS e AGE plot(PAS~SMK, xlab="SMK (sim ou não)", ylab="PAS (mmHg)", col="blue", data=dados) abline(RegModel.3, col="red") #reta ajustada no diagrama de dispersão PAS e SMK # Intervalo de confiança para os parâmetros do modelo #-------------------------------------------------------------------------------------------------------- confint(RegModel.1) # PAS + QUET confint(RegModel.2) # PAS + AGE confint(RegModel.3) # PAS + SMK confint(RegModel.4) # PAS + AGE + QUET confint(RegModel.5) # PAS + AGE + QUET + SMK #-------------------------------------------------------------------------------------------------------- # ANÁLISE DE RESÍDUOS #-------------------------------------------------------------------------------------------------------- # Obtendo o resíduo padronizado #-------------------------------------------------------------------------------------------------------- mean(residuals(RegModel.1)) mean(residuals(RegModel.2)) mean(residuals(RegModel.3)) mean(residuals(RegModel.4)) mean(residuals(RegModel.5)) sd(residuals(RegModel.1)) sd(residuals(RegModel.2)) sd(residuals(RegModel.3)) sd(residuals(RegModel.4)) sd(residuals(RegModel.5)) respad1 <- residuals(RegModel.1) / sd(residuals(RegModel.1)) respad1 respad2 <- residuals(RegModel.2) / sd(residuals(RegModel.2)) respad2 respad3 <- residuals(RegModel.3) / sd(residuals(RegModel.3)) respad3 respad4 <- residuals(RegModel.4) / sd(residuals(RegModel.4)) respad4 respad5 <- residuals(RegModel.5) / sd(residuals(RegModel.5)) respad5 # Diagrama de dispersão dos resíduos x variável explicativa plot(dados$QUET, residuals(RegModel.1), xlab="QUET (kg/m³)", ylab="Resíduo") abline(h=0) plot(dados$QUET, respad1, xlab="QUET (kg/m³)", ylab="Resíduo Padronizado") abline(h=0) plot(dados$QUET, rstandard(RegModel.1), xlab="QUET (kg/m³)", ylab="Resíduo Studentizado") abline(h=0) plot(dados$AGE, residuals(RegModel.2), xlab="AGE (anos)", ylab="Resíduo") abline(h=0) plot(dados$AGE, respad2, xlab="AGE (anos)", ylab="Resíduo Padronizado") abline(h=0) plot(dados$AGE, rstandard(RegModel.2), xlab="AGE (anos)", ylab="Resíduo Studentizado") abline(h=0) plot(dados$SMK, residuals(RegModel.3), xlab="SMK (sim ou não)", ylab="Resíduo") abline(h=0) plot(dados$SMK, respad3, xlab="SMK (sim ou não)", ylab="Resíduo Padronizado") abline(h=0) plot(dados$SMK, rstandard(RegModel.3), xlab="SMK (sim ou não)", ylab="Resíduo Studentizado") abline(h=0) # Box-plot, histograma e qq-plot dos resíduos #-------------------------------------------------------------------------------------------------------- boxplot(rstandard(RegModel.1), ylab="Resíduo Studentizado") hist(rstandard(RegModel.1),scale="density", xlab="Resíduo Studentizado", ylab="Densidade de frequência", main="") qqPlot(rstandard(RegModel.1), dist="norm", xlab="Percentis da N(0,1)", ylab="Resíduo Studentizado") boxplot(rstandard(RegModel.2), ylab="Resíduo Studentizado") hist(rstandard(RegModel.2),scale="density", xlab="Resíduo Studentizado", ylab="Densidade de frequência", main="") qqPlot(rstandard(RegModel.2), dist="norm", xlab="Percentis da N(0,1)", ylab="Resíduo Studentizado") boxplot(rstandard(RegModel.3), ylab="Resíduo Studentizado") hist(rstandard(RegModel.3),scale="density", xlab="Resíduo Studentizado", ylab="Densidade de frequência", main="") qqPlot(rstandard(RegModel.3), dist="norm", xlab="Percentis da N(0,1)", ylab="Resíduo Studentizado") boxplot(rstandard(RegModel.4), ylab="Resíduo Studentizado") hist(rstandard(RegModel.4),scale="density", xlab="Resíduo Studentizado", ylab="Densidade de frequência", main="") qqPlot(rstandard(RegModel.4), dist="norm", xlab="Percentis da N(0,1)", ylab="Resíduo Studentizado") boxplot(rstandard(RegModel.5), ylab="Resíduo Studentizado") hist(rstandard(RegModel.5),scale="density", xlab="Resíduo Studentizado", ylab="Densidade de frequência", main="") qqPlot(rstandard(RegModel.5), dist="norm", xlab="Percentis da N(0,1)", ylab="Resíduo Studentizado") # Default do R #-------------------------------------------------------------------------------------------------------- par(mfrow=c(2,2)) plot(RegModel.1) par(mfrow=c(2,2)) plot(RegModel.2) par(mfrow=c(2,2)) plot(RegModel.3) par(mfrow=c(2,2)) plot(RegModel.4) par(mfrow=c(2,2)) plot(RegModel.5)