O presente post refere-se a um modelo de regressão linear simples com dados em seção cruzada, cuja variável dependente é a expectativa de vida e a variável independente é a renda per capita. Em outras palavras, este modelo pretende averiguar o poder explicativo da renda per capita dos países sobre a expectativa de vida. Foram utilizados dados do Banco Mundial.
###Pcotes necessários
{
library(readxl)
library(ggplot2)
library(dplyr)
library(lmtest)
library(sandwich)
library(stargazer)
}
##
## Anexando pacote: 'dplyr'
## Os seguintes objetos são mascarados por 'package:stats':
##
## filter, lag
## Os seguintes objetos são mascarados por 'package:base':
##
## intersect, setdiff, setequal, union
## Carregando pacotes exigidos: zoo
##
## Anexando pacote: 'zoo'
## Os seguintes objetos são mascarados por 'package:base':
##
## as.Date, as.Date.numeric
##
## Please cite as:
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
Os dados foram baixados diretamente do site do Banco Mundial e organizados em uma planiilha excel, mantendo-se uma amostra de 169 países.
### Extraindo os dados
dados <- read_excel("DataSet_BancoMundial.xlsx")
attach(dados)
dput(dados)
## structure(list(`Country Name` = c("Afghanistan", "Albania", "Algeria",
## "Angola", "Argentina", "Armenia", "Aruba", "Australia", "Austria",
## "Azerbaijan", "Bahamas, The", "Bangladesh", "Barbados", "Belarus",
## "Belgium", "Belize", "Benin", "Bhutan", "Bolivia", "Bosnia and Herzegovina",
## "Botswana", "Brazil", "Brunei Darussalam", "Bulgaria", "Burkina Faso",
## "Burundi", "Cabo Verde", "Cambodia", "Cameroon", "Canada", "Central African Republic",
## "Chad", "Chile", "China", "Colombia", "Comoros", "Congo, Dem. Rep.",
## "Congo, Rep.", "Costa Rica", "Cote d'Ivoire", "Croatia", "Cyprus",
## "Czechia", "Denmark", "Djibouti", "Dominica", "Dominican Republic",
## "Ecuador", "Egypt, Arab Rep.", "El Salvador", "Equatorial Guinea",
## "Estonia", "Eswatini", "Ethiopia", "Fiji", "Finland", "France",
## "Gabon", "Gambia, The", "Georgia", "Germany", "Ghana", "Greece",
## "Guatemala", "Guinea", "Guinea-Bissau", "Guyana", "Haiti", "Honduras",
## "Hungary", "Iceland", "India", "Indonesia", "Iran, Islamic Rep.",
## "Iraq", "Ireland", "Israel", "Italy", "Jamaica", "Japan", "Jordan",
## "Kazakhstan", "Kenya", "Kiribati", "Korea, Rep.", "Kyrgyz Republic",
## "Lao PDR", "Latvia", "Lebanon", "Lesotho", "Liberia", "Libya",
## "Lithuania", "Luxembourg", "Madagascar", "Malawi", "Malaysia",
## "Maldives", "Mali", "Mauritania", "Mauritius", "Mexico", "Micronesia, Fed. Sts.",
## "Moldova", "Mongolia", "Montenegro", "Morocco", "Mozambique",
## "Myanmar", "Namibia", "Nepal", "Netherlands", "New Zealand",
## "Nicaragua", "Niger", "Nigeria", "North Macedonia", "Norway",
## "Oman", "Pakistan", "Panama", "Papua New Guinea", "Paraguay",
## "Peru", "Philippines", "Poland", "Portugal", "Qatar", "Romania",
## "Russian Federation", "Rwanda", "Samoa", "Sao Tome and Principe",
## "Senegal", "Serbia", "Seychelles", "Sierra Leone", "Singapore",
## "Slovak Republic", "Slovenia", "Solomon Islands", "Somalia",
## "South Africa", "Spain", "Sri Lanka", "St. Lucia", "St. Vincent and the Grenadines",
## "Sudan", "Suriname", "Sweden", "Switzerland", "Tajikistan", "Tanzania",
## "Thailand", "Timor-Leste", "Togo", "Tonga", "Trinidad and Tobago",
## "Tunisia", "Turks and Caicos Islands", "Uganda", "Ukraine", "United States",
## "Uruguay", "Uzbekistan", "Vanuatu", "Viet Nam", "Zambia", "Zimbabwe"
## ), `Adjusted net national income per capita` = c(339.682114402906,
## 5262.91359620389, 2846.64050751421, 1232.94823320045, 8993.54741048778,
## 4023.25142894731, 24994.772987977, 45710.8904505974, 43281.7334797461,
## 3965.42439576196, 24192.5278094747, 2541.71174641412, 14810.0109491981,
## 5812.02174051461, 41743.0142804966, 4823.4817941724, 1210.8868870626,
## 2937.61750425413, 2625.31992475194, 5971.24253225669, 5061.46917681828,
## 5647.13391471309, 23060.0681318064, 10114.5235213919, 666.858377543748,
## 155.914279774751, 2820.86578831135, 1374.93488560507, 1362.77025686338,
## 41986.2328595884, 391.039615426431, 486.337982345425, 11095.8802840128,
## 9015.34800413676, 5094.62555383117, 1446.85922913064, 338.936292799504,
## 1079.5371809839, 10968.7313742907, 2334.27058955257, 14553.4842135772,
## 24687.482421875, 19843.782905675, 58796.3599323691, 2786.2924734597,
## 6519.00038493433, 7367.81018035442, 4392.18984789583, 3234.21182214225,
## 3676.34957226722, 2915.37922276238, 22462.3647535878, 2921.08259985864,
## 798.719118620326, 3749.90030646452, 44416.2001989398, 36489.8040864125,
## 4505.83502872903, 608.964552530311, 4084.97057411529, 42982.481891127,
## 1823.90292195275, 17262.586536424, 4308.91309993237, 865.899329488847,
## 677.418415520699, 6359.17320489518, 1724.94720964418, 2407.05060702826,
## 14799.3838039183, 54204.4243535276, 1916.12574705723, 3224.24989077319,
## 2999.38949002625, 3692.20017995738, 49607.1754153272, 43494.5020003156,
## 29681.5571261882, 4704.2840342503, 30520.0704825716, 3905.90675280648,
## 6545.72172734553, 1796.80669264333, 2791.01701400099, 28058.8301659895,
## 880.467336081222, 1844.84098010349, 16238.8897578519, 3094.62387529833,
## 1028.33206084876, 343.424530178291, 4913.90556321537, 19957.147930238,
## 77781.175704672, 432.303587120567, 567.503935132422, 7867.55481345638,
## 8303.11656234299, 662.140185239304, 1818.26062178628, 8167.69674265462,
## 7498.93242803608, 3691.02659650773, 4753.25524379629, 2792.16696296598,
## 8483.01598085187, 3361.27465820313, 366.255665060131, 1041.6036617187,
## 4018.92605454904, 1126.75208777039, 46419.7071212119, 40518.245535159,
## 1694.98665831936, 574.623671760581, 1702.4824798981, 5996.1109517865,
## 69953.0192338809, 12502.1228676848, 1410.80831010391, 11830.9871827426,
## 1843.51129605938, 5248.38510132467, 5056.79818301417, 3174.41602280327,
## 15098.2396638793, 19318.4675574803, 48152.5517316883, 12194.3636281641,
## 8992.818359375, 668.367607718484, 3458.29313092533, 1909.61892930731,
## 1405.01675087286, 7431.64056276737, 12252.2419252636, 409.685377442495,
## 50689.9927011895, 17386.3393801483, 23493.4748984427, 2126.71787352569,
## 360.693566954162, 5717.12039642913, 25118.4531381716, 3611.59123580456,
## 7607.05491085672, 7181.4146394809, 629.995666503906, 3157.77225744389,
## 51830.6266145746, 69631.6582374907, 928.785971496461, 960.039123535156,
## 5402.52217628701, 427.088746002285, 854.79541443012, 4262.72641245691,
## 13573.561381508, 3227.99929205846, 18837.1770296723, 654.941056595561,
## 3961.404296875, 59006.3585416509, 14591.5820259459, 1429.81292408762,
## 3028.38180760278, 3037.45660462072, 614.647305175744, 1472.25441084827
## ), `Life expectancy` = c(61.982, 76.463, 76.377, 61.643, 75.39,
## 72.043, 74.626, 83.3, 81.190243902439, 69.366, 71.598, 72.381,
## 77.571, 72.3706829268293, 81.790243902439, 70.47, 59.821, 71.815,
## 63.63, 75.3, 61.141, 72.75, 74.642, 71.4634146341463, 59.27,
## 61.663, 74.052, 69.584, 60.333, 81.5870731707317, 53.895, 52.525,
## 78.944, 78.211, 72.83, 63.417, 59.193, 63.519, 77.023, 58.598,
## 76.4243902439025, 81.203, 77.2219512195122, 81.4048780487805,
## 62.305, 72.814, 72.615, 73.67, 70.221, 70.748, 60.594, 76.9439024390244,
## 57.066, 64.975, 67.114, 81.8853658536585, 82.3243902439024, 65.821,
## 62.083, 71.694, 80.790243902439, 63.795, 80.0829268292683, 69.237,
## 58.892, 59.652, 65.673, 63.192, 70.123, 74.1634146341463, 83.1658536585366,
## 67.24, 67.57, 73.875, 70.378, 82.3536585365854, 82.5, 82.6463414634146,
## 70.5, 84.4456097560976, 74.256, 70.23, 61.427, 67.417, 83.5268292682927,
## 71.9, 68.061, 72.9804878048781, 75.047, 53.062, 60.747, 71.911,
## 74.0365853658537, 82.5975609756098, 64.485, 62.904, 74.884, 79.918,
## 58.941, 64.364, 73.680243902439, 70.213, 70.71, 68.846, 70.975,
## 73.8243902439024, 74.042, 59.325, 65.672, 59.269, 68.45, 81.309756097561,
## 82.2073170731707, 73.837, 61.576, 52.676, 73.2463414634146, 83.1634146341463,
## 72.541, 66.098, 76.223, 65.351, 70.262, 72.377, 69.266, 75.5024390243902,
## 81.3780487804878, 79.272, 72.809756097561, 69.900243902439, 66.072,
## 72.767, 67.591, 67.093, 72.780487804878, 73.3975609756097, 60.062,
## 83.0926829268293, 74.6146341463415, 80.6756097560976, 70.348,
## 55.28, 62.341, 83.2292682926829, 76.399, 71.111, 69.629, 65.267,
## 70.274, 83.0560975609756, 83.7512195121951, 71.594, 66.201, 78.715,
## 67.737, 61.619, 70.986, 72.971, 73.772, 74.587, 62.705, 69.6478048780488,
## 76.3292682926829, 75.436, 70.862, 70.449, 73.618, 61.223, 59.253
## ), `log Adjusted net national income per capita` = c(5.82801022262618,
## 8.56844006800524, 7.95389480867446, 7.11716351785511, 9.1042626447934,
## 8.29984566787093, 10.126422001509, 10.7300918515692, 10.6754859653926,
## 8.28536816384218, 10.0937990962625, 7.8405930489443, 9.60305864657272,
## 8.66768376530854, 10.6392873936562, 8.48125131020158, 7.09910833461803,
## 7.98535415910059, 7.87295804369541, 8.69471031408697, 8.52941207128763,
## 8.63890342370246, 10.0458577483636, 9.22172764235211, 6.50257769585033,
## 5.04930636760194, 7.94479913394373, 7.22616165305576, 7.21727486052911,
## 10.6450970544735, 5.96880887307919, 6.18690381912103, 9.31432917282442,
## 9.10668373768068, 8.53594145008684, 7.27715043726229, 5.82581216287487,
## 6.9842876921444, 9.30280390157647, 7.75545473995384, 9.58558570878318,
## 10.1140516096179, 9.89564603331247, 10.9818352260584, 7.93246712853533,
## 8.7824763279726, 8.90487581489351, 8.38758320809665, 8.08154053649003,
## 8.20967557478321, 7.97775518389516, 10.0195965100871, 7.97970957992652,
## 6.68300934281017, 8.22948453366828, 10.7013595512172, 10.5047881604496,
## 8.41312850885483, 6.4117600599901, 8.31506980386921, 10.6685479137859,
## 7.50873394660428, 9.75629681067668, 8.36844097030465, 6.76376865410174,
## 6.51828912560715, 8.75765364865322, 7.45295172589235, 7.78615746250109,
## 9.60234082401507, 10.9005178142351, 7.55806058636537, 8.07845561016784,
## 8.00616404361605, 8.21397781407635, 10.8118907678816, 10.6803898182914,
## 10.2982811597255, 8.45622886918634, 10.3261397947854, 8.27024523834226,
## 8.7865669438223, 7.49376630873422, 7.93416132954103, 10.242058662338,
## 6.78045283023648, 7.52014836311428, 9.69516424671441, 8.03742165117424,
## 6.93569341024639, 5.83896737916284, 8.49982433499157, 9.90134265038918,
## 11.2616547233665, 6.06912808939666, 6.34124768348556, 8.97050259598851,
## 9.02438621274509, 6.4954772936992, 7.50563562078097, 9.0079422316644,
## 8.92251594646416, 8.21365990877693, 8.46658497679138, 7.93457326258353,
## 9.04582132367392, 8.12007554357792, 5.90333162795262, 6.94851678692307,
## 8.29876999531229, 7.02709451448263, 10.7454793704852, 10.6095076586878,
## 7.43543014858555, 6.35371534263463, 7.43984274716197, 8.6988663633519,
## 11.1555791419053, 9.43365373828579, 7.25191808869323, 9.37847740089094,
## 7.51942734516255, 8.56567570851948, 8.52848879184268, 8.06287896458389,
## 9.62233343745617, 9.8688167856496, 10.7821294112881, 9.40872912627607,
## 9.10418157770931, 6.50483833312692, 8.14853043157055, 7.55465898768355,
## 7.24780450402621, 8.91350191586441, 9.41346421353503, 6.0153894930262,
## 10.8334837874772, 9.76344008376071, 10.0644779976804, 7.6623351658412,
## 5.88802875285764, 8.65122053041412, 10.1313580397916, 8.19190373975095,
## 8.93683137343539, 8.87925166763273, 6.44571294079669, 8.05762207614286,
## 10.8557365009153, 11.1509746026392, 6.8338783263683, 6.86697403731393,
## 8.59462119319205, 6.05699182770462, 6.75086215893613, 8.35766423747304,
## 9.5158791636068, 8.07961781007707, 9.84358769771397, 6.48454524165547,
## 8.28435386180855, 10.9854004889733, 9.58820006784862, 7.26529889237962,
## 8.01578369894746, 8.01877580100213, 6.42104861578054, 7.29455011781155
## )), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA,
## -169L))
O próximo passo é rodar o modelo de regressão simples com as duas variáveis selecionadas.
#####Modelo de regressão 1
#rodando o modelo 1
reg1<-lm(`Life expectancy`~`Adjusted net national income per capita`,dados)
summary(reg1)
##
## Call:
## lm(formula = `Life expectancy` ~ `Adjusted net national income per capita`,
## data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.884 -3.805 1.067 3.707 10.174
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.700e+01 5.050e-01 132.67 <2e-16
## `Adjusted net national income per capita` 3.309e-04 2.534e-05 13.06 <2e-16
##
## (Intercept) ***
## `Adjusted net national income per capita` ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.358 on 167 degrees of freedom
## Multiple R-squared: 0.5052, Adjusted R-squared: 0.5022
## F-statistic: 170.5 on 1 and 167 DF, p-value: < 2.2e-16
À primeira vista, o modelo aparenta um resultado interessante com a variável explicativa apresentando significância estatística.
#plot do modelo 1
ggplot(dados, aes(x = `Adjusted net national income per capita`, y = `Life expectancy`)) +
geom_point(color = "blue", size = 2) + # Gráfico de dispersão
geom_smooth(method = "lm", color = "red", se = FALSE, linetype = "dashed") + # Linha de regressão linear
geom_smooth(method = "loess", color = "green", se = FALSE) + # Linha de suavização LOESS
labs(title = "Gráfico 1 de Dispersão com Regressão Linear e LOESS",
x = "Renda Per Capita",
y = "Expectativa de Vida") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Pelo plot do modelo logo se percebe que ele sofre com viés, tendo em vista que é evidente que o modelo não está conforme a hipótese Gauus-Markov de linearidade nos parâmetros.
Para lidar com o problema da não linearidade do modelo 1, foi feito um modelo alternativo com o uso do log da renda per capita dos países.
#modelo 2 com log
reg2<-lm(`Life expectancy`~`log Adjusted net national income per capita`,dados)
summary(reg2)
##
## Call:
## lm(formula = `Life expectancy` ~ `log Adjusted net national income per capita`,
## data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.6311 -2.1945 0.4386 2.2958 8.5732
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 32.6810 1.8715 17.46
## `log Adjusted net national income per capita` 4.5197 0.2188 20.66
## Pr(>|t|)
## (Intercept) <2e-16 ***
## `log Adjusted net national income per capita` <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.039 on 167 degrees of freedom
## Multiple R-squared: 0.7187, Adjusted R-squared: 0.7171
## F-statistic: 426.8 on 1 and 167 DF, p-value: < 2.2e-16
Percebe-se que a variável renda per capita é estatisticamente significativa e que este segundo modelo possui um R² maior que o do modelo 1.
#plot do modelo 2
ggplot(dados, aes(x = `log Adjusted net national income per capita`, y = `Life expectancy`)) +
geom_point(color = "blue", size = 2) + # Gráfico de dispersão
geom_smooth(method = "lm", color = "red", se = FALSE, linetype = "dashed") + # Linha de regressão linear
geom_smooth(method = "loess", color = "green", se = FALSE) + # Linha de suavização LOESS
labs(title = "Gráfico 2 de Dispersão com Regressão Linear e LOESS",
x = "Renda Per Capita",
y = "Expectativa de Vida") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
O scatter plot do modelo 2 ajuda a visualizar que os dados com o log da renda per capita apresentam, apesar de um desvio incial, um ajustamento mais fiel ao modelo linear com uma maior linearidade na tendência geral dos dados.
Em seguida, serão testados alguns dos pressupostos Gauss-Markov para garantir que modelo de MQO esteja adequado.
###pressuposto da exogeneidade
dados$residuos <- resid(reg2)
ggplot(dados, aes(x = `log Adjusted net national income per capita`, y = residuos)) +
geom_point(color = "blue") +
geom_smooth(method = "loess", color = "red", se = FALSE) +
labs(title = "Resíduos vs log(Renda per capita)",
x = "log(Renda per capita)",
y = "Resíduos") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
ggplot(dados, aes(x = fitted(reg2), y = residuos)) +
geom_point(color = "blue") +
geom_smooth(method = "loess", color = "red", se = FALSE) +
labs(title = "Resíduos vs Valores ajustado",
x = "Valores ajustados",
y = "Resíduos") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
Nenhum desses plots evidenciam algum padrão rígido entre os resíduos e o log da renda per capita ou os valores ajustados do modelo, além da maior parte dos resíduos estarem contornando a média zero em ambos os gráficos, o que ajuda a corroborar com a hipótese da exogeneidade do modelo. A exogeneidade é uma hipótese imprtante para que os estimadores dos parâmetros do modelo não sejam viesados.
resettest(reg2)
##
## RESET test
##
## data: reg2
## RESET = 0.87544, df1 = 2, df2 = 165, p-value = 0.4186
Por último, foi feito o teste RESET de Ramsey, que entregou p-value = 0.4186. Nesse caso, a não rejeição de H0 evidencia que não há a presença de varável relevante omitida no modelo, o que também contribui com a exogeneidade.
#Testar heterocedasticidade
bptest(reg2)
##
## studentized Breusch-Pagan test
##
## data: reg2
## BP = 6.7668, df = 1, p-value = 0.009287
Foi realizado o teste de Breusch-Pagan, cujo p-value = 0.009287 indica a rejeição de H0. Com este teste, conclui-se a presença de heterocedasticidade nos resíduos, o que interfere na eficiência dos parâmetros
Outrossim, para garantir a viabilidade dos testes de hipótese dos parâmetros, faz-se necessário verificar a normalidade dos resíduos.
#Histograma dos resíduos
hist(resid(reg2), breaks = 20, col = "skyblue", freq = FALSE,
main = "Histograma dos resíduos com curva normal")
curve(dnorm(x, mean = mean(resid(reg2)), sd = sd(resid(reg2))),
col = "red", lwd = 2, add = TRUE)
O histograma mostra uma distribuição com uma certa simetria, mas com caudas pesadas e que em sua distribuição geral não se adequam perfeitamente a uma normal.
# Teste de normalidade dos resíduos (Shapiro-Wilk)
shapiro.test(resid(reg2))
##
## Shapiro-Wilk normality test
##
## data: resid(reg2)
## W = 0.96879, p-value = 0.0007476
O teste de Shapiro-Wilk com um p-value = 0.0007476, rejeita-se H0, o que evidencia estatisticamente que os resíduos não possuem uma distribuição normal.
# Plot QQ-plot dos resíduos
qqnorm(resid(reg2))
qqline(resid(reg2), col="red")
Por fim, o QQ-plot, que em boa parte do gráfico apresenta um bom ajustamente, evidencia nos valores extremos, principalmente nos inferiores, um desvio. Isso também corrobora com a hipótese de não normalidade.
Diante do problema de heterocedasticidade e de não normalidade dos resíduos, será aplicado o método de erros-padrão robustos no modelo.
# Erros-padrão robustos do tipo HC1 (o mais usado, similar ao Stata)
coeftest(reg2, vcov = vcovHC(reg2, type = "HC1"))
##
## t test of coefficients:
##
## Estimate Std. Error t value
## (Intercept) 32.6810 1.7872 18.286
## `log Adjusted net national income per capita` 4.5197 0.1951 23.166
## Pr(>|t|)
## (Intercept) < 2.2e-16 ***
## `log Adjusted net national income per capita` < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
stargazer(reg2, reg2, type = "text",
se = list(NULL, sqrt(diag(vcovHC(reg2, type = "HC1")))),
column.labels = c("OLS", "Robust SE"))
##
## ==========================================================================
## Dependent variable:
## ----------------------------
## `Life expectancy`
## OLS Robust SE
## (1) (2)
## --------------------------------------------------------------------------
## `log Adjusted net national income per capita` 4.520*** 4.520***
## (0.219) (0.195)
##
## Constant 32.681*** 32.681***
## (1.872) (1.787)
##
## --------------------------------------------------------------------------
## Observations 169 169
## R2 0.719 0.719
## Adjusted R2 0.717 0.717
## Residual Std. Error (df = 167) 4.039 4.039
## F Statistic (df = 1; 167) 426.775*** 426.775***
## ==========================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01