Este trabalho tem como objetivo analisar os fatores que influenciam a mortalidade infantil (variável resposta) em diferentes países, utilizando um modelo de regressão múltipla, a partir das variáveis Expectativa de vida, médicos para cada mil pessoas e continente do país (variáveis explicativas).
O banco de dados que vamos analisar é
Global Country Information Dataset 2023
disponivel em https://www.kaggle.com/datasets/nelgiriyewithana/countries-of-the-world-2023/data
Foram removidos alguns países da análise, entre eles: Eswatini, Nauru, Sudão do Sul, Andorra, Mônaco, Vaticano, Palestina, Tuvalu, pois não haviam informações sobre as variáveis estudadas
library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)
##
## Attaching package: 'rnaturalearthdata'
## The following object is masked from 'package:rnaturalearth':
##
## countries110
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(summarytools)
## system might not have X11 capabilities; in case of errors when using dfSummary(), set st_options(use.x11 = FALSE)
library(kableExtra)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
ggplot(data = world_data) +
geom_sf(aes(fill = mort_inf), color = "white") +
scale_fill_gradient(
low = "green",
high = "red",
na.value = "lightgrey", # Cor para países sem dados
name = "Mortalidade infantil"
) +
labs(
title = "Mapa da taxa de mortalidade infantil",
) +
theme_light() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12),
legend.title = element_text(size = 10)
)
desc_mort <- summarytools::descr(wor$mort_inf)
desc_mort <- desc_mort[c(1,2,3,4,5,6,7)]
desc_mort <- data_frame("Estatística" = c("Média", "Desvio padrão", "Minimo", "Q1", "Mediana", "Q3", "Máximo"), "Valor" = round(desc_mort, 3))
kable(
desc_mort,
col.names = c(" Estatística ", " Valor "),
format = "html",
align = "c"
) %>%
kable_styling(
full_width = FALSE,
position = "center"
)
| Estatística | Valor |
|---|---|
| Média | 21.625 |
| Desvio padrão | 19.888 |
| Minimo | 1.400 |
| Q1 | 6.000 |
| Mediana | 13.900 |
| Q3 | 33.900 |
| Máximo | 84.500 |
ggplot(data = wor)+
geom_histogram(aes(x = mort_inf), fill = "lightblue", col = "black")+
labs(title = "Distribuição da taxa de mortalidade infantil")+
ylab("Frequência")+
xlab("Taxa de mortalidade infantil")+
theme_light()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data = world_data) +
geom_sf(aes(fill = expect_vida), color = "white") +
scale_fill_gradient(
low = "red",
high = "green",
na.value = "lightgrey", # Cor para países sem dados
name = "Expectativa de vida"
) +
labs(
title = "Mapa da expectativa de vida",
) +
theme_light() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12),
legend.title = element_text(size = 10)
)
desc_exp <- summarytools::descr(wor$expect_vida)
desc_exp <- desc_exp[c(1,2,3,4,5,6,7)]
desc_exp <- data_frame("Estatística" = c("Média", "Desvio padrão", "Minimo", "Q1", "Mediana", "Q3", "Máximo"), "Valor" = round(desc_exp, 3))
kable(
desc_exp,
col.names = c(" Estatística ", " Valor "),
format = "html",
align = "c"
) %>%
kable_styling(
full_width = FALSE,
position = "center"
)
| Estatística | Valor |
|---|---|
| Média | 72.208 |
| Desvio padrão | 7.511 |
| Minimo | 52.800 |
| Q1 | 66.700 |
| Mediana | 73.400 |
| Q3 | 77.600 |
| Máximo | 84.200 |
ggplot(data = wor)+
geom_histogram(aes(x = expect_vida), fill = "lightblue", col = "black")+
labs(title = "Distribuição da expectativa de vida")+
ylab("Frequência")+
xlab("Expectativa de vida")+
theme_light()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data = world_data) +
geom_sf(aes(fill = medic), color = "white") +
scale_fill_gradient(
low = "red",
high = "green",
na.value = "lightgrey", # Cor para países sem dados
name = "Médicos por 1000 hab."
) +
labs(
title = "Mapa da Quantidade de Médicos para cada habitante por País",
) +
theme_light() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12),
legend.title = element_text(size = 10)
)
desc_med <- summarytools::descr(wor$mort_inf)
desc_med <- desc_med[c(1,2,3,4,5,6,7)]
desc_med <- data_frame("Estatística" = c("Média", "Desvio padrão", "Minimo", "Q1", "Mediana", "Q3", "Máximo"), "Valor" = round(desc_med, 3))
kable(
desc_med,
col.names = c(" Estatística ", " Valor "),
format = "html",
align = "c"
) %>%
kable_styling(
full_width = FALSE,
position = "center"
)
| Estatística | Valor |
|---|---|
| Média | 21.625 |
| Desvio padrão | 19.888 |
| Minimo | 1.400 |
| Q1 | 6.000 |
| Mediana | 13.900 |
| Q3 | 33.900 |
| Máximo | 84.500 |
ggplot(data = wor)+
geom_histogram(aes(x = medic), fill = "lightblue", col = "black")+
labs(title = "Distribuição da quantidade de médicos para cada 1000 habitantes")+
ylab("Frequência")+
xlab("Médicos / 1000 hab.")+
theme_light()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
cor1 = cor(wor$mort_inf, wor$expect_vida, use = "complete.obs")
cor2 = cor(wor$mort_inf, wor$medic, use = "complete.obs")
cor3 = cor(wor$medic, wor$expect_vida, use = "complete.obs")
cov1 = cov(wor$mort_inf, wor$expect_vida, use = "complete.obs")
cov2 = cov(wor$mort_inf, wor$medic, use = "complete.obs")
cov3 = cov(wor$medic, wor$expect_vida, use = "complete.obs")
cor_matrix <- matrix(c(
1, cor1, cor2,
cor1, 1, cor3,
cor2, cor3, 1
), nrow = 3, byrow = TRUE, dimnames = list(
c("Mortalidade Infantil", "Expectativa de Vida", "Médicos"),
c("Mortalidade Infantil", "Expectativa de Vida", "Médicos")
))
cov_matrix <- matrix(c(
var(wor$mort_inf), cov1, cov2,
cov1, var(wor$expect_vida), cov3,
cov2, cov3, var(wor$medic, na.rm = TRUE)
), nrow = 3, byrow = TRUE, dimnames = list(
c("Mortalidade Infantil", "Expectativa de Vida", "Médicos"),
c("Mortalidade Infantil", "Expectativa de Vida", "Médicos")
))
library(knitr)
kable(cor_matrix, caption = "Tabela de Correlação", digits = 3)
| Mortalidade Infantil | Expectativa de Vida | Médicos | |
|---|---|---|---|
| Mortalidade Infantil | 1.000 | -0.933 | -0.696 |
| Expectativa de Vida | -0.933 | 1.000 | 0.700 |
| Médicos | -0.696 | 0.700 | 1.000 |
kable(cov_matrix, caption = "Tabela de Covariância", digits = 3)
| Mortalidade Infantil | Expectativa de Vida | Médicos | |
|---|---|---|---|
| Mortalidade Infantil | 395.528 | -139.306 | -22.516 |
| Expectativa de Vida | -139.306 | 56.409 | 8.570 |
| Médicos | -22.516 | 8.570 | 2.697 |
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
ggpairs(data = wor, columns = c(6,7,9))
Vamos analisar os diagramas de dispersão entre as variáveis qualitativas que estamos interessados em analisar
Mortalidade infantil vs expectativa de vida
Mortalidade infantil vs quantidade de médicos para cada 1000 pessoas
Além disso, vamos analisar os boxplots da mortalidade infantil para cada continente
ggplot(data = wor) +
geom_point(aes(x = expect_vida, y = mort_inf, color = cont)) +
scale_color_manual(
values = c(
"africa" = "orange",
"oceania" = "blue",
"america" = "green",
"asia" = "yellow",
"europa" = "lightblue"
)
) +
labs(
title = "Dispersão entre expectativa de vida e mortalidade infantil",
x = "Expectativa de vida",
y = "Mortalidade infantil",
color = "Continente"
) +
theme_light()
ggplot(data = wor) +
geom_point(aes(x = medic, y = mort_inf, color = cont)) +
scale_color_manual(
values = c(
"africa" = "orange",
"oceania" = "blue",
"america" = "green",
"asia" = "yellow",
"europa" = "lightblue"
)
) +
labs(
title = "Dispersão entre quantidade de médicos e mortalidade infantil",
x = "Médicos para cada 1000 habitantes",
y = "Mortalidade infantil",
color = "Continente"
) +
theme_light()
Há um problema:
uma relação não linear …
… para resolver, vamos fazer uma transformação
ggplot(data = wor) +
geom_point(aes(x = log(medic+0.1), y = mort_inf, color = cont)) +
scale_color_manual(
values = c(
"africa" = "orange",
"oceania" = "blue",
"america" = "green",
"asia" = "yellow",
"europa" = "lightblue"
)
) +
labs(
title = "Dispersão entre o log da quantidade de médicos e mortalidade infantil",
x = "log(Médicos para cada 1000 habitantes + 0.1 )",
y = "Mortalidade infantil",
color = "Continente"
) +
theme_light()
ggplot(wor)+
geom_boxplot(aes(cont, mort_inf), fill = cores)+
labs(title = "Quartís da mortalidade infantil por continente")+
xlab("Continentes") + ylab("Mortalidade infantil")+
theme_light()
\[ \begin{aligned} \text{Mortalidade infantil} = & \beta_0 \\ &\quad + \beta_1 \cdot \text{Expectativa de vida} \\ &\quad + \beta_2 \cdot \log(\text{Médicos por mil pessoas}+0.1) \\ &\quad + \beta_{3,j} I_{\text{Continente}_j} \hspace{4mm}, \text{para} \ j \in\{0,1,2,3,4\} \\ &\quad + \epsilon \end{aligned} \]
\[ \begin{aligned} Y &:= \text{Mortalidade infantil} \\ X_1 &:= \text{Expectativa de vida} \\ X_2 &:= \log(\text{Médicos para cada 1000 habitantes}+0.1) \\ X_3 &:= \text{Continente} \end{aligned} \]
\[ Y = \beta_0 + \beta_1 \cdot X_1 + \beta_2 \cdot X_2 + \beta_{3,j} I_{continente_j} \]
mod = lm(mort_inf ~ expect_vida + log(medic+0.1) + cont, data = wor)
mod
##
## Call:
## lm(formula = mort_inf ~ expect_vida + log(medic + 0.1) + cont,
## data = wor)
##
## Coefficients:
## (Intercept) expect_vida log(medic + 0.1) contamerica
## 165.6259 -1.9765 -3.4498 -1.2235
## contasia conteuropa contoceania
## -0.7514 -0.9235 -4.9725
\[ \begin{aligned} \beta_0 &= 165.6259 \\ \beta_1 &= -1.9765\\ \beta_2 &= -3.4497\\ X_{3,0} = \text{África} \Rightarrow \beta_{3,0} &= 0\\ X_{3,1} = \text{Ásia} \Rightarrow \beta_{3,1} &= -0.7514\\ X_{3,2} = \text{América} \Rightarrow \beta_{3,2} &= -1.2235\\ X_{3,3} = \text{Europa} \Rightarrow \beta_{3,3} &= -0.9235\\ X_{3,4} = \text{Oceania} \Rightarrow \beta_{3,4} &= -4.9725 \end{aligned} \]
\[ \begin{aligned} \hat{Y} = 165.6259 - 1.9765 \cdot X_1 - 3.4498 \cdot X_2 \\ + \quad \begin{cases} -1.2235 & \text{se } \text{Continente } = \text{América} \\ -0.7514 & \text{se } \text{Continente } = \text{Ásia} \\ -0.9235 & \text{se } \text{Continente } = \text{Europa} \\ -4.9725 & \text{se } \text{Continente } = \text{Oceania} \\ 0 & \text{se } \text{Continente } = \text{África} \end{cases} \end{aligned} \]
separada por continente e separada pelas variáveis quantitativas
af_ex = function(X1) { 165.6259 - 1.9765 * X1 }
as_ex = function(X1) { 164.8745 - 1.9765 * X1 }
am_ex = function(X1) { 164.4024 - 1.9765 * X1 }
eu_ex = function(X1) { 164.7024 - 1.9765 * X1 }
oc_ex = function(X1) { 160.6534 - 1.9765 * X1 }
af_me = function(X2) { 165.6259 - 3.4498 * X2 }
as_me = function(X2) { 164.8745 - 3.4498 * X2 }
am_me = function(X2) { 164.4024 - 3.4498 * X2 }
eu_me = function(X2) { 164.7024 - 3.4498 * X2 }
oc_me = function(X2) { 160.6534 - 3.4498 * X2 }
expectativa_vida <- seq(50, 90, by = 0.1)
data_ex <- data.frame(
Expectativa_Vida = expectativa_vida,
Af = af_ex(expectativa_vida),
As = as_ex(expectativa_vida),
Am = am_ex(expectativa_vida),
Eu = eu_ex(expectativa_vida),
Oc = oc_ex(expectativa_vida)
)
medicos_por_mil <- seq(-3, 3, by = 0.01)
data_me <- data.frame(
Medicos_por_Mil = medicos_por_mil,
Af = af_me(medicos_por_mil),
As = as_me(medicos_por_mil),
Am = am_me(medicos_por_mil),
Eu = eu_me(medicos_por_mil),
Oc = oc_me(medicos_por_mil)
)
library(patchwork)
p1 = ggplot(data_ex, aes(x = Expectativa_Vida)) +
geom_line(aes(y = Af, color = "África")) +
geom_line(aes(y = As, color = "Ásia")) +
geom_line(aes(y = Am, color = "América")) +
geom_line(aes(y = Eu, color = "Europa")) +
geom_line(aes(y = Oc, color = "Oceania")) +
labs(title = "Mortalidade infantil vs \nExpectativa de Vida",
x = "Expectativa de Vida",
y = "Mortalidade infantil",
color = "Continente") +
scale_color_manual(values = c("orange", "green", "yellow", "lightblue", "blue")) +
theme_light() +
theme(legend.position = "bottom")
p2 = ggplot(data_me, aes(x = Medicos_por_Mil)) +
geom_line(aes(y = Af, color = "África")) +
geom_line(aes(y = As, color = "Ásia")) +
geom_line(aes(y = Am, color = "América")) +
geom_line(aes(y = Eu, color = "Europa")) +
geom_line(aes(y = Oc, color = "Oceania")) +
labs(title = "Mortalidade infantil vs \nMédicos por Mil Pessoas",
x = "log(Médicos por Mil Pessoas+0.1)",
y = "Mortalidade infantil",
color = "Continente") +
scale_color_manual(values = c("orange", "green", "yellow", "lightblue", "blue")) +
theme_light() +
theme(legend.position = "bottom")
(p1 + p2) +
plot_layout(guides = "collect") &
theme(legend.position = "bottom")
summary(mod)
##
## Call:
## lm(formula = mort_inf ~ expect_vida + log(medic + 0.1) + cont,
## data = wor)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.8057 -5.3314 -0.3174 4.7588 25.2128
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 165.6259 9.3559 17.703 < 2e-16 ***
## expect_vida -1.9765 0.1340 -14.749 < 2e-16 ***
## log(medic + 0.1) -3.4498 0.9061 -3.807 0.000196 ***
## contamerica -1.2235 1.9222 -0.636 0.525303
## contasia -0.7514 1.7589 -0.427 0.669763
## conteuropa -0.9235 2.2451 -0.411 0.681328
## contoceania -4.9725 2.4990 -1.990 0.048204 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.919 on 171 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8807, Adjusted R-squared: 0.8765
## F-statistic: 210.4 on 6 and 171 DF, p-value: < 2.2e-16
library(plotly)
plot_ly(data = wor,
x = ~expect_vida,
y = ~log(medic + 0.1),
z = ~mort_inf,
color = ~cont,
colors = cores,
type = 'scatter3d',
mode = 'markers',
marker = list(size = 3)) %>%
layout(title = "Expectativa de Vida, Log(Médicos) e Mortalidade Infantil",
scene = list(
xaxis = list(title = "Expectativa de Vida"),
yaxis = list(title = "log(Médicos por 1000 habitantes)"),
zaxis = list(title = "Mortalidade Infantil")
),
legend = list(font = list(size = 14),
itemsizing = 'constant'))
grid <- expand.grid(
expect_vida = seq(min(wor$expect_vida, na.rm = TRUE), max(wor$expect_vida, na.rm = TRUE), length.out = 30),
medic = seq(min(wor$medic, na.rm = TRUE), max(wor$medic, na.rm = TRUE), length.out = 30),
cont = unique(wor$cont, na.rm = TRUE) # Cada continente será representado
)
wor <- wor[-179, ]
# Criando a grade de valores para o gráfico
grid <- expand.grid(
expect_vida = seq(min(wor$expect_vida, na.rm = TRUE), max(wor$expect_vida, na.rm = TRUE), length.out = 30),
medic = seq(min(wor$medic, na.rm = TRUE), max(wor$medic, na.rm = TRUE), length.out = 30),
cont = unique(wor$cont) # Cada continente será representado
)
# Calculando o log de 'medic' e fazendo as previsões
grid$log_medic <- log(grid$medic + 0.1)
grid$pred <- predict(mod, newdata = grid)
# Plot com o plano de regressão
plot_ly(data = wor,
x = ~expect_vida,
y = ~log(medic + 0.1),
z = ~mort_inf,
color = ~cont, # Aplicando a cor do continente
colors = cores,
type = 'scatter3d',
mode = 'markers', # Corrigido para 'mode = markers'
marker = list(size = 3) # Definindo o tamanho do marcador
) %>%
add_trace(data = grid,
x = ~expect_vida,
y = ~log_medic,
z = ~pred,
type = 'mesh3d', # Usando 'mesh3d' para representar o plano
opacity = 0.5,
colorscale = "Viridis",
color = ~cont, # Aplicando a cor do continente também no plano
name = "Plano de Regressão") %>%
layout(title = "Expectativa de Vida, Log(Médicos) e Mortalidade Infantil",
scene = list(
xaxis = list(title = "Expectativa de Vida"),
yaxis = list(title = "log(Médicos por 1000 habitantes)"),
zaxis = list(title = "Mortalidade Infantil")
),
legend = list(font = list(size = 14),
itemsizing = 'constant'))
shapiro.test(mod$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.9839, p-value = 0.03813
\[ H_0 : resíduos \sim normal \\ H_1 : resíduos \nsim normal \]
Em um nível de 1% de significância, não há evidencias para rejeitar a hipótese nula, isto é, podemos dizer que os resíduos seguem uma distribuição normal
Pelo gráfico abaixo, podemos ver que o resíduo se aproxima de uma normal
plot(mod, 2)
\[ H_0 : \rho = 0 \\ H_1 : \rho \neq 0 \]
dwtest(mod, alternative = "two.sided")
##
## Durbin-Watson test
##
## data: mod
## DW = 1.912, p-value = 0.5494
## alternative hypothesis: true autocorrelation is not 0
A autocorrelação dos resíduos não é estatisticamente significativa, isto é, a autocorrelação é proxima de 0, portanto, podemos assumir independencia dos resíduos
summary(rstandard(mod))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.163453 -0.794164 -0.047800 0.001015 0.699533 3.705340
plot(mod, 4)
Austrália, Nova Zelandia, Kuwait, Chade e o Paquistão são outliers
\[ H_0 : resíduos \rightarrow homocedástico \\ H_1 : resíduos \nrightarrow homocedástico \]
library(lmtest)
bptest(mod)
##
## studentized Breusch-Pagan test
##
## data: mod
## BP = 30.708, df = 6, p-value = 2.881e-05
Podemos perceber que o pressuposto da homocedasticidade não foi atendido, pois o p-valor é muito baixo O pressuposto foi quebrado, o modelo pode não ser bom/suficiente
Vamos verificar se a heterocedasticidade é dada pelos outliers
Vamos remover a Austrália, Nova Zelandia, Kuwait, Chade e o Paquistão, rodar o modelo novamente e fazer o teste de Breusch-Pagan
wor_so <- wor[-c(25, 37,95,148, 121), ]
mod_so = lm(mort_inf ~ expect_vida + log(medic+0.1) + cont, data = wor_so)
bptest(mod_so)
##
## studentized Breusch-Pagan test
##
## data: mod_so
## BP = 32.537, df = 6, p-value = 1.287e-05
O p-valor continua muito próximo de 0, portanto, não são os outliers que estão atrapalhando a homocedasticidade dos resíduos, mas o modelo é insuficiente para explicar a mortalidade infantil
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
vif(mod)
## GVIF Df GVIF^(1/(2*Df))
## expect_vida 3.685555 1 1.919780
## log(medic + 0.1) 4.030801 1 2.007685
## cont 2.563368 4 1.124868
Como os valores de \(VIF_j\) estão bem abaixo de 10, podemos assumir que não há multicolinearidade
A análise indicou a presença de não \(homocedasticidade\), evidenciada pela variação não constante dos resíduos ao longo dos valores preditos, o que compromete a precisão das estimativas. O \(R^2\) ajustado foi moderado, sugerindo que o modelo explica uma fração significativa da variação na taxa de mortalidade infantil, mas ainda há variação não explicada. Isso reforça a necessidade de considerar outras variáveis ou utilizar um modelo mais sofisticado para melhorar o ajuste e a precisão das previsões.