Laboratorio_Renato
Objetivo:
- Ver quais variáveis impactam mais a taxa de latrocinio
Dicionário de dados:
tx_latrc - taxa de latrocÃnios por 100 mil hab
gini_ibge - coeficiente de gini
raz_1040 - razão de renda entre os 10% mais ricos e os 40% mais pobres
raz_2020 - razão de renda entre os 20% mais ricos e os 20% mais pobres
Obs-1: (a), (b) e (c) medem o mesmo fenômeno, a desigualdade. A variável principal é o Gini, mas coloquei outras 2 para ter mais opções. Obs-2: O gini é a opção principal porque, para ele os valores chegam até 2015, enquanto para a razão de renda chegam apenas até 2014. Obs-3: O IBGE muda, a partir de 2015, a metodologia para calcular renda e, por isso, a forma de cálculo do gini (de 2016 em diante) e das razões de renda (de 2015 em diante) são levemente discrepantes. Por isso, encerrei o painel em 2015.taxa de desligamentos - número de desligamentos / população total Obs: imagino que possua um alto nÃvel colinearidade com as variáveis que medem desigualdade.
tx_pres - presos por 100 mil hab
perc_jov (1) - percentual de jovens (15-24 anos) na população
perc_hom - percentual de homens na população
pbf - gasto per capita com o programa bolsa famÃlia
densidade urbana (1) - população total / área urbana da UF
densidade urbana (2) - (população total * taxa de urbanização) / área urbana Obs: (j) e (k) medem o mesmo fenômeno, a densidade urbana.
taxa de casamentos - número de casamentos / população total
require(DT) #para criar tabelas
require(tidyverse) #conjunto de pacotes no R
require(readxl) #para ler xls & xlsx
#install.packages("factoextra")
require(factoextra) #para fazer o PCA para selecionar as variaveis importantes
require(lme4) #modelo de efeitos mistos
require(redres) #para analise dos residuos no modelo misto
require(car)Planilha de dados
dados=read_excel("dados_renato.xlsx",sheet=1)
names(dados)=c("Ano","UF","Cod_UF","tx_latrc","tx_pres","gini_ibge","perc_jov_15_24","perc_hom","pbf","densidade_urbana1","densidade_urbana2","taxa_casamentos","Taxa_desligamentos","raz_2020","raz_1040")
datatable(dados,options = list(pageLength = 5))dados_comp = dados[complete.cases(dados), ] %>%
filter(tx_latrc != "NA") %>%
mutate(tx_latrc = as.numeric(tx_latrc),
Cod_UF=as.factor(Cod_UF))
dados_PCA = dados_comp %>%
select(-c("UF","Cod_UF"))- 2 dimensões respondem por ~ 60% dos dados
res.pca <- prcomp(dados_PCA, scale = TRUE)
fviz_eig(res.pca)
fviz_pca_biplot(res.pca, repel = TRUE,
col.var = "#2E9FDF", # Variables color
col.ind = "#696969", # Individuals color
label = "var"
)
# Variaveis Eveliny tx_pres + perc_jov_15_24 + densidade_urbana1 + raz_1040 x1 + x3 + x6 + x11 + (1 | id))
=> Falta Ano
Variaveis principais: Ano, gini_ibge, Cod_UF + perc_jov_15_24 + perc_hom + densidade_urbana2 + Taxa_desligamentos +
Escolhendo as variáveis principais - há problema de normalidade
model1 = lm(tx_latrc ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1 + raz_1040, data=dados_comp)
car::S(model1)## Call: lm(formula = tx_latrc ~ Ano + tx_pres + perc_jov_15_24 +
## densidade_urbana1 + raz_1040, data = dados_comp)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.247e+01 3.470e+01 -1.512 0.132
## Ano 2.626e-02 1.718e-02 1.529 0.128
## tx_pres 1.502e-04 4.255e-04 0.353 0.724
## perc_jov_15_24 1.841e-02 2.574e-02 0.715 0.475
## densidade_urbana1 6.221e-04 1.690e-03 0.368 0.713
## raz_1040 2.286e-02 1.395e-02 1.638 0.103
##
## Residual standard deviation: 0.6375 on 219 degrees of freedom
## Multiple R-squared: 0.02547
## F-statistic: 1.145 on 5 and 219 DF, p-value: 0.3377
## AIC BIC
## 443.85 467.76car::residualPlots(model1)
## Test stat Pr(>|Test stat|)
## Ano 1.2684 0.205998
## tx_pres -0.8960 0.371215
## perc_jov_15_24 0.9129 0.362279
## densidade_urbana1 1.4316 0.153704
## raz_1040 3.3145 0.001075 **
## Tukey test 1.4211 0.155295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::vif(model1)## Ano tx_pres perc_jov_15_24 densidade_urbana1
## 1.459538 1.502251 1.373669 1.621390
## raz_1040
## 1.391901plot(model1)
shapiro.test(residuals(model1))##
## Shapiro-Wilk normality test
##
## data: residuals(model1)
## W = 0.88091, p-value = 2.645e-12#Transformação de Box cox apropriada
a=boxCox(model1)
summary(a)## Length Class Mode
## x 100 -none- numeric
## y 100 -none- numericmaximo=cbind(a$x,a$y)
datatable(maximo)maximo[which.max(maximo[,2])]## [1] 0.3434343#valor para lambda=0.3434343Verificando para o Estado de SP
- Não rejeita a normalidade mas os resÃduos parecem não ser homoscedásticos
dados_sp = dados_comp %>%
filter(Cod_UF=="20")
model_sp = lm(tx_latrc ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1 + raz_1040, data=dados_sp)
car::S(model_sp)## Call: lm(formula = tx_latrc ~ Ano + tx_pres + perc_jov_15_24 +
## densidade_urbana1 + raz_1040, data = dados_sp)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.616e+02 1.873e+02 -2.464 0.0905 .
## Ano 2.317e-01 9.498e-02 2.440 0.0925 .
## tx_pres -8.944e-03 5.443e-03 -1.643 0.1989
## perc_jov_15_24 6.893e-02 1.521e-01 0.453 0.6812
## densidade_urbana1 -4.407e-02 4.396e-02 -1.003 0.3899
## raz_1040 1.546e-01 9.132e-02 1.693 0.1891
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 0.05786 on 3 degrees of freedom
## Multiple R-squared: 0.898
## F-statistic: 5.282 on 5 and 3 DF, p-value: 0.1006
## AIC BIC
## -21.64 -20.26car::residualPlots(model_sp)
## Test stat Pr(>|Test stat|)
## Ano -4.5379 0.04529 *
## tx_pres -0.6522 0.58121
## perc_jov_15_24 1.9031 0.19736
## densidade_urbana1 -0.7927 0.51106
## raz_1040 3.5047 0.07265 .
## Tukey test -8.6695 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(model_sp)
shapiro.test(residuals(model_sp))##
## Shapiro-Wilk normality test
##
## data: residuals(model_sp)
## W = 0.85972, p-value = 0.09521hist(dados_sp$tx_latrc)
Modelo com Box Cox & \(\lambda = 0.3434343\)
- ainda com problemas de normalidade
model2 = lm((tx_latrc^0.3434343-1)/0.3434343 ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1 + raz_1040, data=dados_comp)
car::S(model2)## Call: lm(formula = (tx_latrc^0.3434343 - 1)/0.3434343 ~ Ano + tx_pres +
## perc_jov_15_24 + densidade_urbana1 + raz_1040, data = dados_comp)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.400e+01 3.227e+01 -1.983 0.0486 *
## Ano 3.162e-02 1.597e-02 1.980 0.0490 *
## tx_pres 1.264e-04 3.957e-04 0.319 0.7497
## perc_jov_15_24 2.065e-03 2.393e-02 0.086 0.9313
## densidade_urbana1 9.824e-04 1.571e-03 0.625 0.5324
## raz_1040 1.784e-02 1.297e-02 1.375 0.1706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 0.5928 on 219 degrees of freedom
## Multiple R-squared: 0.02992
## F-statistic: 1.351 on 5 and 219 DF, p-value: 0.244
## AIC BIC
## 411.10 435.02car::residualPlots(model2)
## Test stat Pr(>|Test stat|)
## Ano 1.3167 0.18932
## tx_pres -0.4603 0.64575
## perc_jov_15_24 0.7998 0.42471
## densidade_urbana1 1.6523 0.09991 .
## raz_1040 2.4547 0.01488 *
## Tukey test 1.1298 0.25857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(model2)
shapiro.test(residuals(model2))##
## Shapiro-Wilk normality test
##
## data: residuals(model2)
## W = 0.98259, p-value = 0.007206hist((dados_comp$tx_latrc^0.3-1)/0.3)
Transformação Box Cox com dois parametros: https://www.statisticshowto.com/box-cox-transformation
Testa-se a transformação para um modelo linear simples (com apenas Ano na variável explicativa)
require("geoR")## Loading required package: geoR## Warning: package 'geoR' was built under R version 4.0.3## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.8-1 (built on 2020-02-08) is now loaded
## --------------------------------------------------------------bc2 <- boxcoxfit(dados_comp$tx_latrc,dados_comp$Ano, lambda2 = TRUE)
lambda1 <- bc2$lambda[1]
lambda2 <- bc2$lambda[2]
lambda1## lambda
## -2.437589lambda2## lambda2
## 3.728376ytrans = ((dados_comp$tx_latrc + lambda2)^lambda1-1)/lambda1
hist(ytrans)
Modelo com Box Cox de dois parametros
- Verificou-se Normalidade dos residuos & aparente homoscedasticidade
l1=-2.437589
l2=3.728376
dados_comp = dados_comp %>%
mutate(ytrans = ((tx_latrc + l2)^l1-1)/l1)
model3 = lm(ytrans ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1 + raz_1040, data=dados_comp)
par(mfrow=c(3,3))
car::S(model3)## Call: lm(formula = ytrans ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1
## + raz_1040, data = dados_comp)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.083e-02 1.387e-01 0.655 0.5133
## Ano 1.532e-04 6.867e-05 2.232 0.0267 *
## tx_pres 3.128e-07 1.701e-06 0.184 0.8543
## perc_jov_15_24 7.028e-06 1.029e-04 0.068 0.9456
## densidade_urbana1 2.815e-06 6.754e-06 0.417 0.6772
## raz_1040 9.687e-05 5.577e-05 1.737 0.0838 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 0.002548 on 219 degrees of freedom
## Multiple R-squared: 0.03613
## F-statistic: 1.642 on 5 and 219 DF, p-value: 0.15
## AIC BIC
## -2041.08 -2017.17car::residualPlots(model3)
## Test stat Pr(>|Test stat|)
## Ano 1.3133 0.19047
## tx_pres -0.2780 0.78126
## perc_jov_15_24 0.8649 0.38803
## densidade_urbana1 1.5477 0.12315
## raz_1040 2.3094 0.02186 *
## Tukey test 1.2245 0.22077
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#AIC = -2524.768
plot(model3)
shapiro.test(residuals(model3))##
## Shapiro-Wilk normality test
##
## data: residuals(model3)
## W = 0.99226, p-value = 0.2837hist(dados_comp$ytrans)
AIC(model3)## [1] -2041.078
Verificando para o Estado de SP
- transformação Box Cox 2 parametros
- Gráficos de normalidade & homoscedasticidade estão ruins mas deve ser pelo fato de poucas observações.
l1=-2.437589
l2=3.728376
dados_sp = dados_sp %>%
mutate(ytrans = ((tx_latrc + l2)^l1-1)/l1)
model2_sp = lm(ytrans ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1 + raz_1040, data=dados_sp)
car::S(model2_sp)## Call: lm(formula = ytrans ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1
## + raz_1040, data = dados_sp)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.5409047 1.1461188 -2.217 0.1134
## Ano 0.0014753 0.0005811 2.539 0.0848 .
## tx_pres -0.0000594 0.0000333 -1.784 0.1725
## perc_jov_15_24 0.0003361 0.0009308 0.361 0.7420
## densidade_urbana1 -0.0002858 0.0002690 -1.063 0.3659
## raz_1040 0.0009927 0.0005587 1.777 0.1737
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 0.000354 on 3 degrees of freedom
## Multiple R-squared: 0.8931
## F-statistic: 5.011 on 5 and 3 DF, p-value: 0.1075
## AIC BIC
## -113.38 -112.00car::residualPlots(model2_sp)
## Test stat Pr(>|Test stat|)
## Ano -5.0221 0.03744 *
## tx_pres -0.5959 0.61172
## perc_jov_15_24 1.8304 0.20868
## densidade_urbana1 -0.7205 0.54605
## raz_1040 3.1381 0.08830 .
## Tukey test -7.4659 8.274e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::vif(model2_sp)## Ano tx_pres perc_jov_15_24 densidade_urbana1
## 221.57083 235.19445 62.50251 14.68077
## raz_1040
## 32.95654plot(model2_sp)
shapiro.test(residuals(model2_sp))##
## Shapiro-Wilk normality test
##
## data: residuals(model2_sp)
## W = 0.85692, p-value = 0.08874hist(dados_sp$ytrans)
Modelo Linear misto
- selecionei apenas as variáveis de interesse no complete cases
- O efeito aleatorio aparece no intercepto & inclinaçao da reta do tempo.
l1=-2.437589
l2=3.728376
dados_comp2 = dados %>%
select(Cod_UF,tx_latrc,Ano,tx_pres,perc_jov_15_24,densidade_urbana1,raz_1040) %>%
filter(tx_latrc != "NA") %>%
mutate(tx_latrc = as.numeric(tx_latrc),
ytrans = ((tx_latrc + l2)^l1-1)/l1,
Cod_UF = as.factor(Cod_UF))
dados_comp2 = dados_comp2[complete.cases(dados_comp2), ]
modelo_misto <- lmer(ytrans ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1 + raz_1040 + (Ano | Cod_UF), dados_comp2)## boundary (singular) fit: see ?isSingularmodelo_misto## Linear mixed model fit by REML ['lmerMod']
## Formula: ytrans ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1 +
## raz_1040 + (Ano | Cod_UF)
## Data: dados_comp2
## REML criterion at convergence: -2035.502
## Random effects:
## Groups Name Std.Dev. Corr
## Cod_UF (Intercept) 1.811e-03
## Ano 1.839e-06 -1.00
## Residual 1.811e-03
## Number of obs: 225, groups: Cod_UF, 27
## Fixed Effects:
## (Intercept) Ano tx_pres perc_jov_15_24
## 1.545e-01 1.217e-04 1.067e-06 3.733e-05
## densidade_urbana1 raz_1040
## 7.880e-06 -1.600e-06
## convergence code 0; 0 optimizer warnings; 1 lme4 warningsrc_resids <- compute_redres(modelo_misto)## Loading required namespace: testthatplot(rc_resids,main="resÃduos condicionais versus Ãndices")
plot_resqq(modelo_misto)
shapiro.test(rc_resids)##
## Shapiro-Wilk normality test
##
## data: rc_resids
## W = 0.99493, p-value = 0.6587# creates normal quantile plots for each random effect
random <- ranef(modelo_misto)
aleatorio1 = random[["Cod_UF"]][["(Intercept)"]]
plot(aleatorio1,main="efeitos aleatórios versus Ãndices")
shapiro.test(aleatorio1)##
## Shapiro-Wilk normality test
##
## data: aleatorio1
## W = 0.98141, p-value = 0.8927aleatorio2 = random[["Cod_UF"]][["Ano"]]
plot(aleatorio2,main="efeitos aleatórios2 versus Ãndices")
shapiro.test(aleatorio2)##
## Shapiro-Wilk normality test
##
## data: aleatorio2
## W = 0.98141, p-value = 0.8927plot_ranef(modelo_misto)
random## $Cod_UF
## (Intercept) Ano
## 1 -6.047734e-04 6.142843e-07
## 2 9.039128e-04 -9.181277e-07
## 3 1.456990e-04 -1.479904e-07
## 4 2.944495e-03 -2.990790e-06
## 5 -3.633221e-03 3.690344e-06
## 6 -1.777237e-03 1.805178e-06
## 7 -6.640034e-05 6.744447e-08
## 8 -2.768174e-04 2.811698e-07
## 9 3.108264e-03 -3.157135e-06
## 10 -3.509963e-04 3.565159e-07
## 11 1.284615e-03 -1.304812e-06
## 12 1.273771e-03 -1.293795e-06
## 13 -1.039033e-03 1.055373e-06
## 14 -1.338967e-03 1.360017e-06
## 15 -2.114846e-04 2.148072e-07
## 16 1.978564e-05 -2.009728e-08
## 17 2.566769e-03 -2.607123e-06
## 18 1.331081e-03 -1.352011e-06
## 19 -8.537882e-04 8.672140e-07
## 20 9.758565e-04 -9.911991e-07
## 21 1.499520e-03 -1.523095e-06
## 22 1.127105e-03 -1.144827e-06
## 23 -4.516287e-04 4.587300e-07
## 24 2.992658e-04 -3.039727e-07
## 25 -2.335953e-03 2.372680e-06
## 26 -1.483854e-03 1.507182e-06
## 27 -3.055985e-03 3.104034e-06
##
## with conditional variances for "Cod_UF"AIC(modelo_misto)## [1] -2015.502#-2717.188Valor-p do modelo linear misto
- Verificamos problema de matriz singular - não verifiquei como resolver
- O ano e a taxa de desligamentos apresentaram significancia estatÃstica
https://www.r-bloggers.com/2014/02/three-ways-to-get-parameter-specific-p-values-from-lmer/ Para Obter os valores p
require(lmerTest)## Loading required package: lmerTest##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step# re-fit model
modelo_misto2 <- lmer(ytrans ~ Ano + tx_pres + perc_jov_15_24 + densidade_urbana1 + raz_1040 + (Ano | Cod_UF),dados_comp2, REML = FALSE)## boundary (singular) fit: see ?isSingular# get Satterthwaite-approximated degrees of freedom
a <- coef(summary(modelo_misto2))[, 3]
# get approximate p-values
b <- coef(summary(modelo_misto2))[, 5]
resultados=cbind(a,b)
names(resultados)=c("coeficiente","valor-p")
resultados## a b
## (Intercept) 176.43522 0.22840345
## Ano 173.56139 0.05230941
## tx_pres 41.96871 0.68387814
## perc_jov_15_24 153.41814 0.72347021
## densidade_urbana1 49.48465 0.51102764
## raz_1040 202.02655 0.99772149
## attr(,"names")
## [1] "coeficiente" "valor-p" NA NA NA
## [6] NA NA NA NA NA
## [11] NA NAPendências no assunto Modelos Lineares com Efeitos Mistos:
- O pacote redres exibe vários tipos de resÃduos - verificar as propriedades destes resÃduos e se todos precisam exibir normalidade e homogeneidade das variâncias;
- Fazer os gráficos da variável resposta versus covariáveis
Fazendo o gráfico de dispersão & correlograma
- Correlações baixÃssimas variáveis não correlacionadas com a resposta
- No entanto para o estado de SP verificamos relações interessantes e não tão baixas. Isto nos leva a crer que pode haver interação entre Estado & variáveis, e também pensar que o estado não apresenta efeito aleatório de indivÃduo (dados em painel), e sim um efeito principal.
- O próximo passo poderia ser agrupar estados com caracterÃsticas similares - em macrorregiões territoriais ou por distâncias - k-means; para então ajustar um modelo linear.
correl = dados_comp %>%
select(-c("Ano","Cod_UF","UF"))
pairs(correl)
require(corrplot)## Loading required package: corrplot## corrplot 0.84 loadedcorrplot.mixed(cor(correl), lower.col = "black", number.cex = .7)
## Para o estado de SP
correl = dados_comp %>%
filter(Cod_UF=="20") %>%
select(-c("Ano","Cod_UF","UF"))
pairs(correl)
require(corrplot)
corrplot.mixed(cor(correl), lower.col = "black", number.cex = .7)