Introdução à Análise de Dados e Modelação Estatística
Lab#4 – Modelos de regressão logística
1 Poços no Bangladesh
Muitos dos poços de água potável usados no Bangladesh e noutros países do sul da Ásia estão contaminados com arsénico natural (\(\ce{As2O3}\)), o que se estima que afete cerca de 100 milhões de pessoas. O arsénico é um veneno cumulativo, o que faz com que o risco de cancro e de outras doenças aumente proporcionalmente com a exposição.
Qualquer povoação pode ter poços com níveis de arsénico variados. Infelizmente, mesmo que um poço vizinho seja seguro, isso não quer dizer que o seu poço também o seja. Em contrapartida, se o seu poço tiver um elevado nível de arsénico, é possível que possa encontrar um poço seguro nas proximidades – desde que esteja disposto a percorrer a distância necessária e o seu vizinho aceite partilhar o poço dele.
Uma equipa de investigação mediu o nível de arsénico em todos os poços numa pequena região do Bangladesh e classificou cada um deles como “seguro” (menos de 50 μg/L) ou “inseguro” (pelo menos 50 μg/L). Os agregados familiares que se abasteciam de poços inseguros foram encorajados a mudar para poços seguros, de vizinhos ou comunitários, ou a construirem novos poços. Alguns anos mais tarde, os investigadores voltaram à região para averiguarem quem tinha mudado de poço.
O ficheiro de dados encontra-se disponível em:
https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/wells.csv (82.3KB).
Nas cinco colunas desse ficheiro encontram-se os valores das seguintes variáveis:
switch: indica se um agregado familiar mudou ou não de poço (variável resposta)
arsenic: nível de arsénico do poço do agregado (100 μg/L)
dist: distância até ao poço seguro mais próximo (m)
assoc: indica se algum membro do agregado participa ativamente em associações comunitárias
educ: nível de educação do representante do agregado (anos de ensino)
Neste trabalho, pretendem-se ajustar possíveis modelos de regressão logística para relacionar a resposta binária com as restantes informações recolhidas sobre cada agregado familiar.
1.1 Análise exploratória
Leia, inspecione o ficheiro de dados e converta a variável dist para hectómetros. Analise as associações entre as covariáveis e a variável resposta.
library(tidyverse)
theme_set(theme_light())
url <- "https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/wells.csv"
poços <- read.csv(url)
str(poços)'data.frame': 3020 obs. of 5 variables:
$ switch : int 1 1 0 1 1 1 1 1 1 1 ...
$ arsenic: num 2.36 0.71 2.07 1.15 1.1 3.9 2.97 3.24 3.28 2.52 ...
$ dist : num 16.8 47.3 21 21.5 40.9 ...
$ assoc : int 0 0 0 0 1 1 1 0 1 1 ...
$ educ : int 0 0 10 12 14 9 4 10 0 0 ...poços$assoc <- as.factor(poços$assoc)
poços$switch <- as.factor(poços$switch)
poços$dist <- poços$dist / 100
summary(poços) switch arsenic dist assoc educ
0:1283 Min. :0.51 Min. :0.00387 0:1743 Min. : 0.00
1:1737 1st Qu.:0.82 1st Qu.:0.21117 1:1277 1st Qu.: 0.00
Median :1.30 Median :0.36762 Median : 5.00
Mean :1.66 Mean :0.48332 Mean : 4.83
3rd Qu.:2.20 3rd Qu.:0.64041 3rd Qu.: 8.00
Max. :9.65 Max. :3.39531 Max. :17.00 ggplot(poços) +
geom_boxplot(aes(arsenic, color = switch)) +
theme(legend.position = "top")
ggplot(poços) +
geom_boxplot(aes(dist, color = switch)) +
theme(legend.position = "top")
tab <- with(poços, table(switch, assoc))
prop.table(tab, margin = 2) assoc
switch 0 1
0 0.4096 0.4456
1 0.5904 0.5544poços |>
summarise(n=n(), .by = c(switch, educ)) |>
ggplot() +
geom_col(aes(educ, n, fill = switch))
poços |>
summarise(n=n(), .by = c(switch, educ)) |>
filter(n > 1) |>
mutate(prop = n / sum(n), .by = educ) |>
ggplot() +
geom_col(aes(educ, prop, fill = switch))
1.2 Regressão logística
Comece por ajustar o modelo de regressão logística de primeira ordem com todas as covariáveis e, de seguida, considere um modelo com todas as interações de segunda ordem.
modelo1 <- glm(switch ~ ., family = "binomial", data = poços) summary(modelo1)Call: glm(formula = switch ~ ., family = "binomial", data = poços) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.15671 0.09960 -1.57 0.12 arsenic 0.46702 0.04160 11.23 < 2e-16 *** dist -0.89611 0.10458 -8.57 < 2e-16 *** assoc1 -0.12430 0.07697 -1.61 0.11 educ 0.04245 0.00959 4.43 9.5e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 4118.1 on 3019 degrees of freedom Residual deviance: 3907.8 on 3015 degrees of freedom AIC: 3918 Number of Fisher Scoring iterations: 4modelo2 <- glm(switch ~ .^2, family = "binomial", data = poços) summary(modelo2)Call: glm(formula = switch ~ .^2, family = "binomial", data = poços) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.0577 0.1765 -0.33 0.7435 arsenic 0.5337 0.0945 5.65 1.6e-08 *** dist -1.2161 0.2793 -4.35 1.3e-05 *** assoc1 0.1389 0.1892 0.73 0.4629 educ -0.0132 0.0214 -0.62 0.5365 arsenic:dist -0.1100 0.1032 -1.07 0.2862 arsenic:assoc1 -0.1626 0.0842 -1.93 0.0536 . arsenic:educ 0.0174 0.0110 1.58 0.1132 dist:assoc1 0.2188 0.2148 1.02 0.3083 dist:educ 0.0838 0.0268 3.13 0.0018 ** assoc1:educ -0.0276 0.0198 -1.39 0.1641 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 4118.1 on 3019 degrees of freedom Residual deviance: 3883.4 on 3009 degrees of freedom AIC: 3905 Number of Fisher Scoring iterations: 4Com base nos dois modelos ajustados procure um modelo reduzido que pareça adequado. Analise a variação de medidas da qualidade de ajustamento dos sucessivos modelos e compare os modelos com o maior e o menor números de graus de liberdade usando um teste de hipóteses apropriado.
modelo3 <- glm(switch ~ . - assoc + dist:educ, family = "binomial", data = poços) summary(modelo3)Call: glm(formula = switch ~ . - assoc + dist:educ, family = "binomial", data = poços) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.000496 0.109615 0.00 0.9964 arsenic 0.480599 0.041987 11.45 < 2e-16 *** dist -1.389852 0.171884 -8.09 6.2e-16 *** educ -0.002077 0.015255 -0.14 0.8917 dist:educ 0.095636 0.025680 3.72 0.0002 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 4118.1 on 3019 degrees of freedom Residual deviance: 3896.2 on 3015 degrees of freedom AIC: 3906 Number of Fisher Scoring iterations: 4modelo4 <- glm(switch ~ arsenic + dist + dist:educ, family = "binomial", data = poços) summ4 <- summary(modelo4) summ4Call: glm(formula = switch ~ arsenic + dist + dist:educ, family = "binomial", data = poços) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.00974 0.07978 -0.12 0.9 arsenic 0.48041 0.04196 11.45 < 2e-16 *** dist -1.37566 0.13647 -10.08 < 2e-16 *** dist:educ 0.09292 0.01612 5.77 8.2e-09 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 4118.1 on 3019 degrees of freedom Residual deviance: 3896.2 on 3016 degrees of freedom AIC: 3904 Number of Fisher Scoring iterations: 4medidas <- data.frame(deviance = modelo1$deviance, aic = modelo1$aic, gl = modelo1$df.residual) medidas |> bind_rows(data.frame(deviance = modelo2$deviance, aic = modelo2$aic, gl = modelo2$df.residual)) -> medidas medidas |> bind_rows(data.frame(deviance = modelo3$deviance, aic = modelo3$aic, gl = modelo3$df.residual)) -> medidas medidas |> bind_rows(data.frame(deviance = modelo4$deviance, aic = modelo4$aic, gl = modelo4$df.residual)) -> medidas medidasdeviance aic gl 1 3908 3918 3015 2 3883 3905 3009 3 3896 3906 3015 4 3896 3904 3016anova(modelo4, modelo2)Analysis of Deviance Table Model 1: switch ~ arsenic + dist + dist:educ Model 2: switch ~ (arsenic + dist + assoc + educ)^2 Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 3016 3896 2 3009 3883 7 12.8 0.078 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1.3 Predição
Separe o conjunto de dados em dois subconjuntos: treino e teste de uma forma repetível. Este último conjunto deverá conter 10% das observações.
Recorra ao modelo reduzido anterior para predizer os valores da variável switch das observações no conjunto de teste. Use o valor de corte 0.5 para classificar cada agregado familiar de acordo com a probabilidade de mudança estimada e calcule as proporções de classificações corretas.
set.seed(1234)
n <- nrow(poços)
indices <- sample(n, size = 0.1 * n)
treino <- poços[-indices,]
teste <- poços[indices,]
modelo <- glm(switch ~ arsenic + dist + dist:educ, family = "binomial", data = treino)
pred <- predict(modelo, teste, type = "response")
prediçoes <- data.frame(obs = teste$switch, pred = (pred > 0.5))
prediçoes$pred <- as.factor(prediçoes$pred)
levels(prediçoes$pred) <- levels(prediçoes$obs)
tab <- with(prediçoes, table(obs, pred))
prop.table(tab, margin = 1) pred
obs 0 1
0 0.2932 0.7068
1 0.1716 0.8284