Introdução à Análise de Dados e Modelação Estatística
Lab#2 – Modelos de regressão linear
1 Diabetes
Foram registados os valores de 10 variáveis clínicas para cada um de 442 doentes de diabetes, assim como uma medida quantitativa da progressão da doença ao fim de um ano (Y).
Os dados encontram-se disponíveis em:
https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/diabetes.txt (18.1KB),
e a descrição das variáveis (em Inglês) pode ser lida aqui.
1.1 Leitura e análise exploratória
Leia o ficheiro de dados para o
Re inspecione a data frame criada. Melhore a estrutura dos dados, reparando possíveis incorreções.url <- "https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/diabetes.txt" db <- read.delim(url) str(db)'data.frame': 442 obs. of 11 variables: $ AGE: int 59 48 72 24 50 23 36 66 60 29 ... $ SEX: int 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 ...db$SEX <- as.factor(db$SEX) summary(db)AGE SEX BMI BP S1 Min. :19.0 1:235 Min. :18.0 Min. : 62.0 Min. : 97 1st Qu.:38.2 2:207 1st Qu.:23.2 1st Qu.: 84.0 1st Qu.:164 Median :50.0 Median :25.7 Median : 93.0 Median :186 Mean :48.5 Mean :26.4 Mean : 94.6 Mean :189 3rd Qu.:59.0 3rd Qu.:29.3 3rd Qu.:105.0 3rd Qu.:210 Max. :79.0 Max. :42.2 Max. :133.0 Max. :301 S2 S3 S4 S5 S6 Min. : 41.6 Min. :22.0 Min. :2.00 Min. :3.26 Min. : 58.0 1st Qu.: 96.0 1st Qu.:40.2 1st Qu.:3.00 1st Qu.:4.28 1st Qu.: 83.2 Median :113.0 Median :48.0 Median :4.00 Median :4.62 Median : 91.0 Mean :115.4 Mean :49.8 Mean :4.07 Mean :4.64 Mean : 91.3 3rd Qu.:134.5 3rd Qu.:57.8 3rd Qu.:5.00 3rd Qu.:5.00 3rd Qu.: 98.0 Max. :242.4 Max. :99.0 Max. :9.09 Max. :6.11 Max. :124.0 Y Min. : 25 1st Qu.: 87 Median :140 Mean :152 3rd Qu.:212 Max. :346Analise a associação entre a variável resposta,
Y, e cada uma das restantes variáveis, recorrendo a medidas numéricas e a representações gráficas.library(ggplot2) library(dplyr) db |> tidyr::pivot_longer(-c(Y, SEX), names_to = "Covariável", values_to = "Valor") |> ggplot() + geom_point(aes(Valor, Y, color = SEX), size = 1/3) + facet_wrap(~ Covariável, scales = "free_x") + theme_light() + theme(legend.position = "top")
db |> select(-SEX) |> cor()AGE BMI BP S1 S2 S3 S4 S5 S6 AGE 1.00000 0.1851 0.3354 0.26006 0.2192 -0.07518 0.2038 0.2708 0.3017 BMI 0.18508 1.0000 0.3954 0.24978 0.2612 -0.36681 0.4138 0.4462 0.3887 BP 0.33543 0.3954 1.0000 0.24246 0.1855 -0.17876 0.2577 0.3935 0.3904 S1 0.26006 0.2498 0.2425 1.00000 0.8967 0.05152 0.5422 0.5155 0.3257 S2 0.21924 0.2612 0.1855 0.89666 1.0000 -0.19646 0.6598 0.3184 0.2906 S3 -0.07518 -0.3668 -0.1788 0.05152 -0.1965 1.00000 -0.7385 -0.3986 -0.2737 S4 0.20384 0.4138 0.2577 0.54221 0.6598 -0.73849 1.0000 0.6179 0.4172 S5 0.27077 0.4462 0.3935 0.51550 0.3184 -0.39858 0.6179 1.0000 0.4647 S6 0.30173 0.3887 0.3904 0.32572 0.2906 -0.27370 0.4172 0.4647 1.0000 Y 0.18789 0.5865 0.4415 0.21202 0.1741 -0.39479 0.4305 0.5659 0.3825 Y AGE 0.1879 BMI 0.5865 BP 0.4415 S1 0.2120 S2 0.1741 S3 -0.3948 S4 0.4305 S5 0.5659 S6 0.3825 Y 1.0000GGally::ggcorr(db)
GGally::ggpairs(db) + theme_light()
1.2 Regressão linear no
A função
lmé usada noRpara ajustar modelos de regressão, incluindo os modelos multivariados. Use essa função para ajustar um MRLS com a covariávelBMIe explore o objeto devolvido pela função.modelo <- lm(Y ~ BMI, data = db) modeloCall: lm(formula = Y ~ BMI, data = db) Coefficients: (Intercept) BMI -117.8 10.2sm <- summary(modelo) smCall: lm(formula = Y ~ BMI, data = db) Residuals: Min 1Q Median 3Q Max -164.92 -43.57 -8.65 46.34 154.88 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -117.773 18.019 -6.54 1.8e-10 *** BMI 10.233 0.674 15.19 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 62.5 on 440 degrees of freedom Multiple R-squared: 0.344, Adjusted R-squared: 0.342 F-statistic: 231 on 1 and 440 DF, p-value: <2e-16sm$coefficientsEstimate Std. Error t value Pr(>|t|) (Intercept) -117.77 18.0189 -6.536 1.752e-10 BMI 10.23 0.6738 15.187 3.466e-42sm$sigma[1] 62.52sm$r.squared[1] 0.3439Calcule intervalos de confiança a 96% para os parâmetros de regressão e para o valor esperado da variável resposta para pessoas em que os valores de
BMIsão iguais a 25 e 40.confint(modelo, level = 0.96)2 % 98 % (Intercept) -154.890 -80.66 BMI 8.845 11.62novos_dados <- data.frame(BMI = c(25, 40)) predict(modelo, novos_dados, level = 0.96, interval = "confidence")fit lwr upr 1 138.1 131.6 144.5 2 291.6 271.7 311.4predict(modelo, novos_dados, level = 0.96, interval = "prediction")fit lwr upr 1 138.1 9.123 267.0 2 291.6 161.255 421.8
1.3 A especificação de modelos
No R, os modelos são especificados numa forma simbólica compacta. Da documentação:
The
~operator is basic in the formation of such models. An expression of the formy ~ modelis interpreted as a specification that the responseyis modelled by a linear predictor specified symbolically bymodel. Such a model consists of a series of terms separated by+operators. The terms themselves consist of variable and factor names separated by:operators. Such a term is interpreted as the interaction of all the variables and factors appearing in the term.
Defina os seguintes modelos:
um modelo de primeira ordem com todas as covariáveis exceto
SEX;f1 <- Y ~ . - SEX terms(f1, data = db)Y ~ (AGE + SEX + BMI + BP + S1 + S2 + S3 + S4 + S5 + S6) - SEX attr(,"variables") list(Y, AGE, SEX, BMI, BP, S1, S2, S3, S4, S5, S6) attr(,"factors") AGE BMI BP S1 S2 S3 S4 S5 S6 Y 0 0 0 0 0 0 0 0 0 AGE 1 0 0 0 0 0 0 0 0 SEX 0 0 0 0 0 0 0 0 0 BMI 0 1 0 0 0 0 0 0 0 BP 0 0 1 0 0 0 0 0 0 S1 0 0 0 1 0 0 0 0 0 S2 0 0 0 0 1 0 0 0 0 S3 0 0 0 0 0 1 0 0 0 S4 0 0 0 0 0 0 1 0 0 S5 0 0 0 0 0 0 0 1 0 S6 0 0 0 0 0 0 0 0 1 attr(,"term.labels") [1] "AGE" "BMI" "BP" "S1" "S2" "S3" "S4" "S5" "S6" attr(,"order") [1] 1 1 1 1 1 1 1 1 1 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") <environment: R_GlobalEnv>um modelo para o logaritmo da variável resposta que inclua as variáveis
BMIeBPe a sua interação;f2 <- log(Y) ~ BMI * BP terms(f2, data = db)log(Y) ~ BMI * BP attr(,"variables") list(log(Y), BMI, BP) attr(,"factors") BMI BP BMI:BP log(Y) 0 0 0 BMI 1 0 1 BP 0 1 1 attr(,"term.labels") [1] "BMI" "BP" "BMI:BP" attr(,"order") [1] 1 1 2 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") <environment: R_GlobalEnv>um modelo com as covariáveis
S1atéS6e todas as suas interações de segunda ordem;f3 <- Y ~ (S1 + S2 + S3 + S4 + S5 + S6)^2 terms(f3, data = db)Y ~ (S1 + S2 + S3 + S4 + S5 + S6)^2 attr(,"variables") list(Y, S1, S2, S3, S4, S5, S6) attr(,"factors") S1 S2 S3 S4 S5 S6 S1:S2 S1:S3 S1:S4 S1:S5 S1:S6 S2:S3 S2:S4 S2:S5 S2:S6 Y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 S1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 S2 0 1 0 0 0 0 1 0 0 0 0 1 1 1 1 S3 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 S4 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 S5 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 S6 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 S3:S4 S3:S5 S3:S6 S4:S5 S4:S6 S5:S6 Y 0 0 0 0 0 0 S1 0 0 0 0 0 0 S2 0 0 0 0 0 0 S3 1 1 1 0 0 0 S4 1 0 0 1 1 0 S5 0 1 0 1 0 1 S6 0 0 1 0 1 1 attr(,"term.labels") [1] "S1" "S2" "S3" "S4" "S5" "S6" "S1:S2" "S1:S3" "S1:S4" [10] "S1:S5" "S1:S6" "S2:S3" "S2:S4" "S2:S5" "S2:S6" "S3:S4" "S3:S5" "S3:S6" [19] "S4:S5" "S4:S6" "S5:S6" attr(,"order") [1] 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") <environment: R_GlobalEnv>um modelo com as covariáveis
BMI,BPe(S1+S2+S3).f4 <- Y ~ BMI + BP + I(S1 + S2 + S3) terms(f4, data = db)Y ~ BMI + BP + I(S1 + S2 + S3) attr(,"variables") list(Y, BMI, BP, I(S1 + S2 + S3)) attr(,"factors") BMI BP I(S1 + S2 + S3) Y 0 0 0 BMI 1 0 0 BP 0 1 0 I(S1 + S2 + S3) 0 0 1 attr(,"term.labels") [1] "BMI" "BP" "I(S1 + S2 + S3)" attr(,"order") [1] 1 1 1 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") <environment: R_GlobalEnv>
Notas
a matriz de delineamento pode ser obtida com a função
model.matrix.md1 <- model.matrix(f1, data = db) # Retirar o comentário da linha abaixo se gostar de olhar para grandes matrizes. # md1as fórmulas também podem ser escritas como texto.
f <- paste("Y ~", paste0("V", 1:20, collapse = " + ")) terms(formula(f))Y ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V11 + V12 + V13 + V14 + V15 + V16 + V17 + V18 + V19 + V20 attr(,"variables") list(Y, V1, V2, V3, V4, V5, V6, V7, V8, V9, V10, V11, V12, V13, V14, V15, V16, V17, V18, V19, V20) attr(,"factors") V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20 Y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 V1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 V2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 V3 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 V4 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 V5 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 V6 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 V7 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 V8 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 V9 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 V10 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 V11 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 V12 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 V13 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 V14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 V15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 V16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 V17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 V18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 V19 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 V20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 attr(,"term.labels") [1] "V1" "V2" "V3" "V4" "V5" "V6" "V7" "V8" "V9" "V10" "V11" "V12" [13] "V13" "V14" "V15" "V16" "V17" "V18" "V19" "V20" attr(,"order") [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") <environment: R_GlobalEnv>
1.4 Ajustamento de modelos
Ajuste o modelo definido em 1.3 a. e analise os resultados.
model1 <- lm(f1, data = db)
summary(model1)
Call:
lm(formula = f1, data = db)
Residuals:
Min 1Q Median 3Q Max
-147.04 -39.03 -4.17 38.20 152.40
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -363.899 68.142 -5.34 1.5e-07 ***
AGE -0.121 0.220 -0.55 0.583
BMI 6.004 0.721 8.32 1.1e-15 ***
BP 0.951 0.225 4.23 2.9e-05 ***
S1 -0.981 0.582 -1.68 0.093 .
S2 0.658 0.539 1.22 0.223
S3 0.514 0.794 0.65 0.518
S4 4.660 6.037 0.77 0.441
S5 68.947 15.927 4.33 1.9e-05 ***
S6 0.203 0.277 0.73 0.465
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 55 on 432 degrees of freedom
Multiple R-squared: 0.501, Adjusted R-squared: 0.49
F-statistic: 48.1 on 9 and 432 DF, p-value: <2e-16par(mfrow = c(2, 2))
plot(model1)
par(mfrow = c(1, 1))1.5 Treino e teste
Comece por separar o conjunto de dados em dois, o conjunto de treino e o conjunto de teste. Este último deve incluir 20% das observações, sorteadas ao acaso mas de uma forma reprodutível.
set.seed(2025) n <- nrow(db) indices <- sample(n, size = 0.2 * n) db.teste <- db[indices,] db.treino <- db[-indices,] # Alternativa com dplyr # db.teste <- slice_sample(db, prop = 0.2) # db.treino <- anti_join(db, db.teste)Nota: seria razoável investigar se há diferenças notórias entre os dois conjuntos.
summary(db.treino) summary(db.teste)Use o modelo anterior para predizer o valor da variável resposta de cada observação no conjunto de teste. Avalie o poder preditivo do modelo usando a raíz do desvio quadrático médio (root mean square error) entre os valores preditos e observados da variável resposta.
\[RMSE = \sqrt{\frac{\sum_{i=1}^n \left(y_i-\hat{y}_i\right)^2}{n}}\]
observado <- db.teste$Y treino <- lm(f1, data = db.treino) pred <- predict(treino, db.teste) desv <- pred - observado summary(desv)Min. 1st Qu. Median Mean 3rd Qu. Max. -126.24 -46.67 1.56 -4.16 38.72 95.14rmse <- sqrt(mean(desv^2)) rmse[1] 55.09res <- data.frame(Y = observado, Desvios = desv, Predições = pred) ggplot(res) + geom_abline(intercept = 0, slope = 1) + geom_point(aes(Y, Predições), color = "tomato") + coord_fixed() + theme_light()
ggplot(res) + geom_hline(yintercept = 0) + geom_point(aes(Y, Desvios), color = "tomato") + theme_light()
ComentárioO modelo tem um fraco poder preditivo. Temos de fazer melhor!
