Introdução à Análise de Dados e Modelação Estatística
Lab#3 – Modelos de regressão linear
1 Mais diabetes . . .
Voltemos aos dados dos pacientes de diabetes.
url <- "https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/diabetes.txt"
db <- read.delim(url)
db$SEX <- as.factor(db$SEX)
str(db)'data.frame': 442 obs. of 11 variables:
$ AGE: int 59 48 72 24 50 23 36 66 60 29 ...
$ SEX: Factor w/ 2 levels "1","2": 2 1 2 1 1 1 2 2 2 1 ...
$ BMI: num 32.1 21.6 30.5 25.3 23 22.6 22 26.2 32.1 30 ...
$ BP : num 101 87 93 84 101 89 90 114 83 85 ...
$ S1 : int 157 183 156 198 192 139 160 255 179 180 ...
$ S2 : num 93.2 103.2 93.6 131.4 125.4 ...
$ S3 : num 38 70 41 40 52 61 50 56 42 43 ...
$ S4 : num 4 3 4 5 4 2 3 4.55 4 4 ...
$ S5 : num 4.86 3.89 4.67 4.89 4.29 ...
$ S6 : int 87 69 85 89 80 68 82 92 94 88 ...
$ Y : int 151 75 141 206 135 97 138 63 110 310 ...As covariáveis serão referidas por nome e, quando conveniente, por número de acordo com a tabela seguinte.
1 2 3 4 5 6 7 8 9 10
AGE SEX BMI BP S1 S2 S3 S4 S5 S6 1.1 Testes de hipóteses
Considere o modelo de primeira ordem com todas as covariáveis e teste a sua significância explicando todos os resultados relevantes.
modelo <- lm(Y ~ ., data = db) summ <- summary(modelo) summCall: lm(formula = Y ~ ., data = db) Residuals: Min 1Q Median 3Q Max -155.83 -38.54 -0.23 37.81 151.35 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -357.4268 67.0581 -5.33 1.6e-07 *** AGE -0.0364 0.2170 -0.17 0.8670 SEX2 -22.8596 5.8358 -3.92 0.0001 *** BMI 5.6030 0.7171 7.81 4.3e-14 *** BP 1.1168 0.2252 4.96 1.0e-06 *** S1 -1.0900 0.5733 -1.90 0.0579 . S2 0.7465 0.5308 1.41 0.1604 S3 0.3720 0.7825 0.48 0.6347 S4 6.5338 5.9586 1.10 0.2735 S5 68.4831 15.6697 4.37 1.6e-05 *** S6 0.2801 0.2733 1.02 0.3060 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 54.2 on 431 degrees of freedom Multiple R-squared: 0.518, Adjusted R-squared: 0.507 F-statistic: 46.3 on 10 and 431 DF, p-value: <2e-16Estamos a testar a hipótese \(H_0: \forall i >0: \beta_i=0\) contra \(H_1: \exists i >0: \beta_i\neq 0\) usando um teste-F, ou seja, uma estatística de teste com uma distribuição \(\mathcal{F}\) com \(p-1=10\) e \(n-p=431\) graus de liberdade. O valor observado da estatística de teste é \(F_o=46.27\) que conduz ao valor-p \(P(F>F_o\mid H_0)\approx 0\). Este resultado permite rejeitar \(H_0\) e concluir que as covariáveis têm alguma influência na variável resposta. Contudo, também temos \(R^2=0.5177\) o que mostra que, da forma que foram usadas, as covariáveis apenas contribuem para explicar 52% da variação observada na variável resposta. O valor observado do \(R^2\) ajustado mostra uma pequena penalização pelo uso de todas as 10 covariáveis.
Calculemos o valor de SST (soma dos quadrados total) como uma medida da variação observada na variável resposta:
SST <- var(db$Y) * (nrow(db) - 1)Juntamente com o sumário do output da função
lm, um outro importante instrumento de análise é a tabela ANOVA:anv <- anova(modelo) anvAnalysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) AGE 1 92527 92527 31.55 3.5e-08 *** SEX 1 293 293 0.10 0.75 BMI 1 826955 826955 281.98 < 2e-16 *** BP 1 129312 129312 44.09 9.4e-11 *** S1 1 1791 1791 0.61 0.43 S2 1 5058 5058 1.72 0.19 S3 1 237329 237329 80.93 < 2e-16 *** S4 1 1821 1821 0.62 0.43 S5 1 58856 58856 20.07 9.6e-06 *** S6 1 3080 3080 1.05 0.31 Residuals 431 1263986 2933 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nota: \(SST = SSE + SSR\)
SST[1] 2621009SSE <- tail(anv$`Sum Sq`, 1) SSE[1] 1263986SSR <- sum(head(anv$`Sum Sq`, -1)) SSR[1] 1357023SSR + SSE[1] 2621009Alternativa
modelo_nulo <- lm(Y ~ 1, data = db) anv_nulo <- anova(modelo_nulo) anv_nuloAnalysis of Variance Table Response: Y Df Sum Sq Mean Sq F value Pr(>F) Residuals 441 2621009 5943anova(modelo_nulo, modelo)Analysis of Variance Table Model 1: Y ~ 1 Model 2: Y ~ AGE + SEX + BMI + BP + S1 + S2 + S3 + S4 + S5 + S6 Res.Df RSS Df Sum of Sq F Pr(>F) 1 441 2621009 2 431 1263986 10 1357023 46.3 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Teste a remoção da variável
S3.modelo_red <- update(modelo, . ~ . - S3) anova(modelo_red, modelo)Analysis of Variance Table Model 1: Y ~ AGE + SEX + BMI + BP + S1 + S2 + S4 + S5 + S6 Model 2: Y ~ AGE + SEX + BMI + BP + S1 + S2 + S3 + S4 + S5 + S6 Res.Df RSS Df Sum of Sq F Pr(>F) 1 432 1264649 2 431 1263986 1 663 0.23 0.63Note-se que \(t_o^2 = 0.4754^2 = 0.226\), do sumário em a., é igual a \(F_o\), com ambos os testes a conduzirem ao mesmo valor-p.
Teste a hipótese \(H_0:\beta_3=5\) contra \(H_1:\beta_3\neq 5\).
modelo_red <- update(modelo, . ~ . - BMI + offset(5 * BMI)) anova(modelo_red, modelo)Analysis of Variance Table Model 1: Y ~ AGE + SEX + BP + S1 + S2 + S3 + S4 + S5 + S6 + offset(5 * BMI) Model 2: Y ~ AGE + SEX + BMI + BP + S1 + S2 + S3 + S4 + S5 + S6 Res.Df RSS Df Sum of Sq F Pr(>F) 1 432 1266059 2 431 1263986 1 2073 0.71 0.4Teste a hipótese de todos os coeficientes de regressão associados às covariáveis
S1atéS6serem iguais.modelo_red <- lm(Y ~ AGE + SEX + BMI + BP + I(S1 + S2 + S3 + S4 + S5 + S6), data = db) anova(modelo_red, modelo)Analysis of Variance Table Model 1: Y ~ AGE + SEX + BMI + BP + I(S1 + S2 + S3 + S4 + S5 + S6) Model 2: Y ~ AGE + SEX + BMI + BP + S1 + S2 + S3 + S4 + S5 + S6 Res.Df RSS Df Sum of Sq F Pr(>F) 1 436 1571921 2 431 1263986 5 307935 21 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1.2 Seleção de modelos e diagnóstico
Use os resultados dos testes-t do ajustamento do modelo de primeira ordem com todas as covariáveis para selecionar um modelo reduzido. Considere a inclusão de interações e aplique as técnicas mais comuns de diagnóstico ao último modelo selecionado.
summary(modelo)
Call:
lm(formula = Y ~ ., data = db)
Residuals:
Min 1Q Median 3Q Max
-155.83 -38.54 -0.23 37.81 151.35
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -357.4268 67.0581 -5.33 1.6e-07 ***
AGE -0.0364 0.2170 -0.17 0.8670
SEX2 -22.8596 5.8358 -3.92 0.0001 ***
BMI 5.6030 0.7171 7.81 4.3e-14 ***
BP 1.1168 0.2252 4.96 1.0e-06 ***
S1 -1.0900 0.5733 -1.90 0.0579 .
S2 0.7465 0.5308 1.41 0.1604
S3 0.3720 0.7825 0.48 0.6347
S4 6.5338 5.9586 1.10 0.2735
S5 68.4831 15.6697 4.37 1.6e-05 ***
S6 0.2801 0.2733 1.02 0.3060
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 54.2 on 431 degrees of freedom
Multiple R-squared: 0.518, Adjusted R-squared: 0.507
F-statistic: 46.3 on 10 and 431 DF, p-value: <2e-16mod2 <- lm(Y ~ .^2, data = db)
summary(mod2)
Call:
lm(formula = Y ~ .^2, data = db)
Residuals:
Min 1Q Median 3Q Max
-159.98 -32.05 -2.35 33.30 147.02
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.31e+03 8.31e+02 1.57 0.117
AGE -9.98e+00 6.83e+00 -1.46 0.145
SEX2 -1.10e+02 1.80e+02 -0.61 0.540
BMI -2.61e+01 2.15e+01 -1.21 0.226
BP 3.39e+00 7.32e+00 0.46 0.643
S1 4.71e+01 2.52e+01 1.87 0.062 .
S2 -4.66e+01 2.49e+01 -1.87 0.062 .
S3 -5.28e+01 2.60e+01 -2.03 0.043 *
S4 9.39e+00 1.31e+02 0.07 0.943
S5 -2.18e+02 2.08e+02 -1.05 0.295
S6 -7.36e+00 8.42e+00 -0.87 0.383
AGE:SEX2 1.21e+00 5.27e-01 2.30 0.022 *
AGE:BMI -3.78e-03 6.38e-02 -0.06 0.953
AGE:BP 8.56e-03 2.07e-02 0.41 0.679
AGE:S1 -2.65e-02 6.18e-02 -0.43 0.668
AGE:S2 9.17e-04 5.79e-02 0.02 0.987
AGE:S3 8.25e-02 8.44e-02 0.98 0.329
AGE:S4 5.87e-01 5.96e-01 0.99 0.325
AGE:S5 9.32e-01 1.52e+00 0.61 0.540
AGE:S6 3.17e-02 2.46e-02 1.29 0.199
SEX2:BMI 1.51e+00 1.68e+00 0.90 0.371
SEX2:BP 6.18e-01 5.30e-01 1.17 0.244
SEX2:S1 1.31e+00 1.59e+00 0.82 0.413
SEX2:S2 -1.25e+00 1.45e+00 -0.86 0.392
SEX2:S3 -8.73e-01 2.15e+00 -0.41 0.684
SEX2:S4 -7.71e+00 1.50e+01 -0.51 0.608
SEX2:S5 -3.03e+01 4.07e+01 -0.74 0.457
SEX2:S6 4.57e-01 6.20e-01 0.74 0.462
BMI:BP 1.26e-01 5.71e-02 2.21 0.028 *
BMI:S1 -2.40e-01 1.85e-01 -1.30 0.194
BMI:S2 2.27e-01 1.82e-01 1.25 0.212
BMI:S3 2.32e-01 2.55e-01 0.91 0.365
BMI:S4 -1.81e-01 1.88e+00 -0.10 0.923
BMI:S5 4.74e+00 5.07e+00 0.94 0.350
BMI:S6 5.54e-02 8.09e-02 0.69 0.494
BP:S1 4.29e-02 6.85e-02 0.63 0.532
BP:S2 -3.05e-02 6.63e-02 -0.46 0.646
BP:S3 -4.65e-02 7.95e-02 -0.58 0.559
BP:S4 -2.53e-01 5.38e-01 -0.47 0.639
BP:S5 -8.89e-01 1.86e+00 -0.48 0.634
BP:S6 -3.72e-02 2.48e-02 -1.50 0.134
S1:S2 4.97e-03 5.13e-03 0.97 0.333
S1:S3 -1.22e-02 2.94e-02 -0.41 0.679
S1:S4 -1.71e-01 7.54e-01 -0.23 0.820
S1:S5 -6.17e+00 4.10e+00 -1.51 0.133
S1:S6 -1.32e-02 7.34e-02 -0.18 0.858
S2:S3 4.15e-03 4.18e-02 0.10 0.921
S2:S4 -9.13e-02 7.28e-01 -0.13 0.900
S2:S5 6.40e+00 4.07e+00 1.57 0.116
S2:S6 -2.57e-03 6.85e-02 -0.04 0.970
S3:S4 -4.12e-01 8.08e-01 -0.51 0.610
S3:S5 6.34e+00 4.32e+00 1.47 0.144
S3:S6 7.82e-02 9.95e-02 0.79 0.432
S4:S5 -7.24e+00 2.99e+01 -0.24 0.809
S4:S6 1.05e+00 6.74e-01 1.57 0.118
S5:S6 5.53e-01 1.94e+00 0.29 0.775
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 53.1 on 386 degrees of freedom
Multiple R-squared: 0.585, Adjusted R-squared: 0.525
F-statistic: 9.87 on 55 and 386 DF, p-value: <2e-16mod3 <- lm(Y ~ . + AGE:SEX + BMI:BP, data = db)
summary(mod3)
Call:
lm(formula = Y ~ . + AGE:SEX + BMI:BP, data = db)
Residuals:
Min 1Q Median 3Q Max
-153.91 -39.88 -0.52 31.92 140.12
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -37.8238 119.3984 -0.32 0.75156
AGE -0.6352 0.2863 -2.22 0.02703 *
SEX2 -91.5355 20.1828 -4.54 7.5e-06 ***
BMI -6.3967 3.9315 -1.63 0.10447
BP -2.2674 1.1169 -2.03 0.04296 *
S1 -1.1764 0.5604 -2.10 0.03640 *
S2 0.8919 0.5199 1.72 0.08697 .
S3 0.5080 0.7653 0.66 0.50714
S4 6.7161 5.8244 1.15 0.24951
S5 72.0231 15.3260 4.70 3.5e-06 ***
S6 0.3276 0.2672 1.23 0.22086
AGE:SEX2 1.3977 0.3968 3.52 0.00047 ***
BMI:BP 0.1255 0.0403 3.11 0.00198 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 52.9 on 429 degrees of freedom
Multiple R-squared: 0.542, Adjusted R-squared: 0.529
F-statistic: 42.3 on 12 and 429 DF, p-value: <2e-16mod4 <- lm(Y ~ . - S3 - S4 - S6 + AGE:SEX + BMI:BP, data = db)
summary(mod4)
Call:
lm(formula = Y ~ . - S3 - S4 - S6 + AGE:SEX + BMI:BP, data = db)
Residuals:
Min 1Q Median 3Q Max
-156.27 -39.52 -0.26 31.82 138.07
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.7770 103.5687 -0.02 0.9863
AGE -0.5905 0.2836 -2.08 0.0379 *
SEX2 -88.8838 20.1273 -4.42 1.3e-05 ***
BMI -6.3048 3.9269 -1.61 0.1091
BP -2.2568 1.1147 -2.02 0.0435 *
S1 -1.0241 0.2176 -4.71 3.4e-06 ***
S2 0.8932 0.2269 3.94 9.6e-05 ***
S5 74.3075 7.1652 10.37 < 2e-16 ***
AGE:SEX2 1.3688 0.3961 3.46 0.0006 ***
BMI:BP 0.1258 0.0403 3.12 0.0019 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 52.9 on 432 degrees of freedom
Multiple R-squared: 0.539, Adjusted R-squared: 0.529
F-statistic: 56 on 9 and 432 DF, p-value: <2e-16plot(mod4, which = c(1, 2, 5))
par(mfrow = c(1, 2))
res <- rstudent(mod4)
hist(res)
boxplot(res)
par(mfrow = c(1, 1))2 Maquinaria pesada
Um estudo experimental foi realizado numa fábrica para investigar como o desgaste de uma peça de uma máquina (y, quantitativa, medida em mm) varia com a velocidade de operação (x1, quantitativa, medida em rot/seg) e o tipo de máquina (x2, qualitativa com 3 níveis: M1, M2 e M3). Os dados estão disponíveis no ficheiro:
https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/machines.txt (255B)
Leia os dados e represente-os graficamente.
url <- "https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/machines.txt" mp <- read.delim(url, sep = " ") summary(mp)i y x1 x2 Min. : 1.00 Min. : 4.30 Min. :200 Length:18 1st Qu.: 5.25 1st Qu.: 5.97 1st Qu.:200 Class :character Median : 9.50 Median : 7.75 Median :300 Mode :character Mean : 9.50 Mean : 7.77 Mean :300 3rd Qu.:13.75 3rd Qu.: 9.30 3rd Qu.:400 Max. :18.00 Max. :12.70 Max. :400mp$i <- NULLlibrary(ggplot2) g <- ggplot(mp) + geom_point(aes(x1, y, color = x2), size = 3, alpha = 0.5) + theme_light() g
Ajuste um modelo de regressão apenas com a covariável quantitativa e sobreponha a reta de regressão ajustada aos dados.
mod0 <- lm(y ~ x1, data = mp) summary(mod0)Call: lm(formula = y ~ x1, data = mp) Residuals: Min 1Q Median 3Q Max -1.972 -1.143 -0.114 0.899 2.736 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.19722 1.25644 0.95 0.35 x1 0.02192 0.00404 5.42 5.6e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.4 on 16 degrees of freedom Multiple R-squared: 0.648, Adjusted R-squared: 0.626 F-statistic: 29.4 on 1 and 16 DF, p-value: 5.63e-05coefs <- coef(mod0) g + geom_abline(intercept = coefs[1], slope = coefs[2])
Ajuste um modelo com ambas as covariáveis usando os valores (0, 1, 2) para codificar os níveis da variável categórica. Produza um gráfico do modelo ajustado e comente.
mp$x3 <- ifelse(mp$x2 == "M1", 0, ifelse(mp$x2 == "M2", 1, 2)) mod1 <- lm(y ~ x1 + x3, data = mp) summary(mod1)Call: lm(formula = y ~ x1 + x3, data = mp) Residuals: Min 1Q Median 3Q Max -1.972 -0.993 0.257 0.817 1.961 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.42222 1.19624 0.35 0.729 x1 0.02192 0.00366 5.98 2.5e-05 *** x3 0.77500 0.36627 2.12 0.051 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.27 on 15 degrees of freedom Multiple R-squared: 0.729, Adjusted R-squared: 0.692 F-statistic: 20.1 on 2 and 15 DF, p-value: 5.64e-05coef <- coef(mod1) coef(Intercept) x1 x3 0.42222 0.02192 0.77500# primeiras 3 cores do ggplot cores <- scales::hue_pal()(3) g + geom_abline(intercept = coef[1], slope = coef[2], color = cores[1]) + geom_abline(intercept = coef[1] + coef[3], slope = coef[2], color = cores[2]) + geom_abline(intercept = coef[1] + 2 * coef[3], slope = coef[2], color = cores[3])
Ajuste um modelo com ambas as covariáveis usando a codificação de nível de referência para a variável categórica. Produza um gráfico do modelo ajustado.
mp$x4 <- as.factor(mp$x2) contrasts(mp$x4)M2 M3 M1 0 0 M2 1 0 M3 0 1mod2 <- lm(y ~ x1 + x4, data = mp) summary(mod2)Call: lm(formula = y ~ x1 + x4, data = mp) Residuals: Min 1Q Median 3Q Max -1.142 -0.352 -0.025 0.215 1.275 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.10833 0.68665 1.61 0.1288 x1 0.02192 0.00207 10.59 4.6e-08 *** x4M2 -1.28333 0.41406 -3.10 0.0078 ** x4M3 1.55000 0.41406 3.74 0.0022 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.717 on 14 degrees of freedom Multiple R-squared: 0.919, Adjusted R-squared: 0.902 F-statistic: 53 on 3 and 14 DF, p-value: 6.88e-08coef <- coef(mod2) coef(Intercept) x1 x4M2 x4M3 1.10833 0.02192 -1.28333 1.55000g + geom_abline(intercept = coef[1], slope = coef[2], color = cores[1]) + geom_abline(intercept = coef[1] + coef[3], slope = coef[2], color = cores[2]) + geom_abline(intercept = coef[1] + coef[4], slope = coef[2], color = cores[3])
Ajuste um modelo com ambas as covariáveis usando a codificação de nível de referência para a variável categórica e incluindo interações. Produza um gráfico do modelo ajustado e compare com o modelo anterior.
mod3 <- lm(y ~ x1 * x4, data = mp) summary(mod3)Call: lm(formula = y ~ x1 * x4, data = mp) Residuals: Min 1Q Median 3Q Max -1.0083 -0.2250 -0.0792 0.4292 0.9917 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.15833 1.14229 1.01 0.33 x1 0.02175 0.00367 5.92 7e-05 *** x4M2 -0.45833 1.61545 -0.28 0.78 x4M3 0.57500 1.61545 0.36 0.73 x1:x4M2 -0.00275 0.00520 -0.53 0.61 x1:x4M3 0.00325 0.00520 0.63 0.54 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.735 on 12 degrees of freedom Multiple R-squared: 0.927, Adjusted R-squared: 0.897 F-statistic: 30.6 on 5 and 12 DF, p-value: 1.98e-06coef <- coef(mod3) coef(Intercept) x1 x4M2 x4M3 x1:x4M2 x1:x4M3 1.15833 0.02175 -0.45833 0.57500 -0.00275 0.00325g + geom_abline(intercept = coef[1], slope = coef[2], color = cores[1]) + geom_abline(intercept = coef[1] + coef[3], slope = coef[2] + coef[5], color = cores[2]) + geom_abline(intercept = coef[1] + coef[4], slope = coef[2] + coef[6], color = cores[3])
anova(mod2, mod3)Analysis of Variance Table Model 1: y ~ x1 + x4 Model 2: y ~ x1 * x4 Res.Df RSS Df Sum of Sq F Pr(>F) 1 14 7.20 2 12 6.48 2 0.722 0.67 0.53
