Perfil dos jovens de 15 a 17 anos assassinados em Minas Gerais no ano de 2013

Estudo via Regressão Logística

1 Introdução

De acordo com Trindade et al. (2015), diversos estudos têm apontado para a existência de um crescimento real da violência no Brasil, em particular das mortes por homicídios, desde o final da década de 1970. As regiões geográficas e seus respectivos municípios, principalmente as grandes cidades, apresentam um aumento na mortalidade por causas externas a partir da década de 1990.

Neste contexto, diversas são as notícias de violência cometida pelos jovens e contra os jovens no país, principalmente na faixa de 15 a 19 anos. De acordo com (“Homicídio é Principal Causa de Mortes de Jovens de 16 E 17 No País,” n.d.) e Waiselfisz (2015), o Brasil ocupa o terceiro lugar em relação a 85 países no ranking de mortes de adolescentes de 15 a 19 anos, perdendo apenas para México e El Salvador. São 54,9 mortes a cada 100 mil jovens. Na faixa de 16 e 17 anos, a taxa de mortalidade ficou em 54,1 homicídios por 100 mil adolescentes em 2013, um aumento de 2,7% em relação a 2012 e de 38,3% na década.

De acordo com Waiselfisz (2015), o homicídio é uma das principais causas de mortes de adolecentes de 16 a 17 anos no país e quando se faz uma análise dos dados de Minas Gerais do sistema de mortalidade do Datasus, conforme Figura 1 é possível notar que o maior número de assassinatos ocorre aos 17 anos de idade. Neste sentido, o presente trabalho tem a finalidade de caracterizar o perfil quanto a estado civil, escolaridade, sexo, raça/cor e idade (15 a 17 anos) dos jovens que morreram assassinados no estado de Minas Gerais no ano de 2013.

2 Material e Métodos

A informação de mortalidade é uma das mais importantes na área da saúde, pois o óbito é um evento único e seu registro obrigatório. No Brasil, o Ministério da Saúde através da Fundação Nacional de Saúde e do Datasus são responsáveis por divulgar os dados de mortalidade do país por meio do sistema de informação da mortalidade (SIM).

A fonte de informação primária da base de dados são os atestados de óbito emitidos pelos cartórios civis. Esta base contém informações sobre a data do óbito, idade, sexo, estado civil, local de ocorrência, causa de mortalidade, município de residência, ocupação e escolaridade. Apesar da enorme quantidade de informações, o banco de dados apresenta problemas sérios de preenchimento de algumas variáveis como educação, estado civil, ocupação, entre outras, que dificultam o seu uso.

Utilizou-se para a análise apenas as variáveis consideradas prioritárias pelo Ministério da Saúde, idade, sexo, estado cívil, escolaridade e causa de mortalidade, nas quais os valores não preenchidos foram retirados do estudo. A causa de mortalidade está codificada segundo a Classificação Internacional de Doenças através da CID10.

As descrições das codificações estão na tabela abaixo:

\[\begin{array}{c|c|c} \hline \textrm{Variáveis} &\textrm{ Categorias} &\textrm{ Descrição} \\ \hline \textrm{Y}& 0& \textrm{Morte registrada como causas distintas de homicídio}\\ & 1& \textrm{Morte registrada como homicídio}\\ \hline \textrm{S} & 1*& \textrm{Masculino}\\ & 2& \textrm{Feminino}\\ \hline \textrm{R} & 1*& \textrm{Raça/Cor Branca}\\ & 2& \textrm{Raça/Cor Negra}\\ & 4& \textrm{Raça/Cor Parda}\\ & 5& \textrm{Raça/Cor Indígena}\\ \hline \textrm{ESC} & 1*& \textrm{Nenhum estudo}\\ & 2& \textrm{1 a 3 anos de estudo}\\ & 3& \textrm{4 a 7 anos de estudo}\\ & 4& \textrm{8 a 11 anos de estudo}\\ & 5& \textrm{12 ou mais anos de estudo}\\ \hline \textrm{I} & \textrm{Idade (Contínua)} &\textrm{15 a 17 anos} \\ \hline \textrm{*Categorias de Referência} \end{array}\]

Como havia poucos indivíduos na categoria distinta de solteiro para a variável estado civil, foram considerados somente indivíduos solteiros na análise. Não houve nenhum indivíduo de 15 a 17 anos caracterizado com a raça/cor amarela.

Foram coletados os dados de 5.418 jovens de 15 a 17 anos que morreram no ano de 2013 no estado de Minas Gerais, para a análise de regressão logística foram retirados 62 indivíduos (aproximadamente 1,15% dos dados) por não terem todas as informações completas e por terem informações cuja categoria era caracterizada como ``ignorado’’ no banco de dados. Dessa forma, restaram 5356 indivíduos para análise.

A análise será feita com o objetivo de analisar a relação entre uma variável dependente e quatro outras variáveis independentes como mencionado na tabela anterior, por meio da teoria de modelos de regressão. Como a variável resposta é binária (0 ou 1), usa-se conceitos de modelos lineares generalizados considerando a distribuição binomial como componente aleatório do modelo. o componente sistemático é dado pela combinação linear das variáveis explicativas e para função de ligação considera-se como possível as funções: logit, probit, complemento log-log e cauchit.

Assim, temos:

\[\begin{eqnarray*} Y_i &\sim& Binomial (n, \pi_i)\\ g(\pi_i) &=& \beta_{0} + \beta_{1}x_{1i} + \beta_{2}x_{2i} + \dots + \beta_{p}x_{pi} \end{eqnarray*}\]

Sendo \(Y\) a variável resposta, \(x_{1i}, x_{2i}, \dots, x_{ni}\) as i-ésimas realizações das respectivas variáveis explicativas \(X_1, X_2, \dots, X_n\) e \(g(\pi_i)\) e a função de ligação, funções citadas acima. Como houve indivíduos com o mesmo conjunto de covariáveis considera-se o número de repetições igual a \(n\) para a distribuição de \(Y_{i}.\)

2.1 Estatística Descritiva dos Dados

No gráfico abaixo temos a proporção de pessoas assassinadas por idade, dos 10 aos 60 anos. Note que o número de homicídios aumenta até os 17 anos e a partir dai começa a diminuir lentamente, haja visto que a morte por assassinato dos 15 aos 17 anos é superior a morte causada por outras causas, como pode ser verificado no segundo gráfico.

dados2=read.table("Dados_Completos_a_partir_dos_dez.txt",header=TRUE)
attach(dados2)

par(mfrow=c(1,1))
a2=table(I,Y)
prop2=a2[,2]/(a2[,1]+a2[,2])
Idade2 <- seq(10, 60, 1)
plot(Idade2,prop2,ylab="Proporção", xlab="Idade", main="Proporção por idade de pessoas \n assassinadas no ano de 2013 em Minas Gerais \n em relação ao total de pessoas que morreram")
lines(Idade2,prop2)
title(sub="Figura 1")

rm(dados2)

Os dados abaixo apresentam a quantidade de indíviduos em cada categoria das respectivas variáveis estudadas.

setwd("~/GitHub/webpage/fsbmat.github.io/blog_posts/17-09-2016")
dados=read.table("15_17.txt",header=TRUE)
attach(dados)
#S=sexo, R=Raça/cor, ESC=Escolaridade, I=Idade, Y=1: Morrer assassinado e 0: morrer de outra forma
dados=transform(R=factor(R),ESC=factor(ESC),S=factor(S),dados)
#head(dados)

#Estatística Descritiva dos dados
kable(summary(dados),format = "markdown",digits = 2, padding = 3)

	Y	S	R	ESC	I
	Min. :0.0000	1:4666	1:1503	1: 64	Min. :15.00
	1st Qu.:0.0000	2: 690	2: 398	2: 890	1st Qu.:16.00
	Median :1.0000	NA	4:3428	3:3030	Median :16.00
	Mean :0.5743	NA	5: 27	4:1351	Mean :16.23
	3rd Qu.:1.0000	NA	NA	5: 21	3rd Qu.:17.00
	Max. :1.0000	NA	NA	NA	Max. :17.00

Na Figura 2 e Figura 3 podemos constatar o que se pode chamar de uma tragédia estadual, morrem mais jovens com idade entre 15 e 17 anos assassinados no estado de Minas Gerais do que devido a outras causas.

tb <- table(dados$Y)
names(tb)<- c("Outras causas", "Assasinados")
bp <- barplot(tb, beside = T, las = 1, xlab = 'discriminação das mortes',
              ylab = 'Frequência', , ylim = c(0, 4000), main="Quantidade de jovens com idades entre 15 e 17 anos \n que morreram em 2013 no Estado de Minas Gerais")
text(x = c(bp), y = c(tb), labels = tb, pos = 3)
title(sub="Figura 2")

#Note que as mortes por assassinato superam as mortes por outras causas

bp=barplot(table(Y,I), beside=T, leg=c("Morte por outra causa", "Homicidío" ),
           args.legend = list(x = "topleft",bty = "n"),ylim=c(0, 2000), ylab="Quantidade", xlab="Idade", main="Quantidade de jovens com idades entre 15 e 17 anos \n que morreram em 2013 no Estado de Minas Gerais")
text(bp, table(Y,I)+100, table(Y,I))
title(sub="Figura 3")

A Figura 4 mostra o número de jovens de 15 a 17 anos, discriminados por sexo, que morreram por motivos distintos de homicídio e por homicídio. Note que o número de jovens, do sexo masculino, nesta faixa etária que morrem assassinados, é maior do que o número de jovens que morrem por outras causas.

bp=barplot(table(Y,S),ylim=c(0, 3500), beside=T,leg=c("Morte por outra causa", "Homicidío" ),
           ylab="Quantidade", xlab="Gênero",names.arg = c("masculino","feminino"), main="Quantidade de jovens com idades entre 15 e 17 anos \n que morreram em 2013 no Estado de Minas Gerais \n discriminados por Gênero")
text(bp, table(Y,S)+150, table(Y,S))
title(sub="Figura 4")

A Figura 5 mostra o número de jovens de 15 a 17 anos, discriminados por raça/cor, que morreram por motivos distintos de homicídio e por homicídio. Note que o número de jovens, de cor negra e/ou parda, nesta faixa etária que morrem assassinados, é maior do que o número de jovens que morrem por outras causas.

#barplot(table(Y,R))
bp=barplot(table(Y,R),args.legend = list(x = "topleft",bty = "n"), ylim=c(0, 2500), beside=T,leg=c("Morte por outra causa", "Homicidío" ),
           ylab="Quantidade", xlab="Raça/Cor",names.arg = c("Branco","Negro","Pardo","Indígena"), main="Quantidade de jovens com idades entre 15 e 17 anos \n que morreram em 2013 no Estado de Minas Gerais \n discriminados por Raça/Cor")
text(bp, table(Y,R)+150, table(Y,R))
title(sub="Figura 5")

A grande maioria dos jovens de 15 a 17 anos que morreram assassinadas em Minas Gerais no ano de 2013 tinham somente de 1 a 7 anos de estudo, como mostra a Figura 6.

bp=barplot(table(Y,ESC), args.legend = list(x = "topleft",bty = "n"), ylim=c(0, 2500), beside=T,leg=c("Morte por outra causa", "Homicidío" ),
           ylab="Quantidade", xlab="Escolaridade",names.arg = c("Nenhum","1 a 3","4 a 7","8 a 11","12 ou mais"), main="Quantidade de jovens com idades entre 15 e 17 anos \n que morreram em 2013 no Estado de Minas Gerais \n discriminados por escolaridade")
text(bp, table(Y,ESC)+150, table(Y,ESC))
title(sub="Figura 6")

3 Ajuste do Melhor Modelo de Regressão Logística e Seleção de Variáveis

Após a descrição e exploração dos dados é necessário buscar a seleção de variáveis que melhor explicam a variável dependente, a seleção será realizada pelo algoritmo stepwise considerando como critério de seleção o AIC (Critério de Informação de Akaike), poderia ser usado também o método forward ou o método backward.

Considerou-se como modelo completo o modelo aditivo com todos os efeitos principais e todas as interações duplas somando ao todo 10 parâmetros. Os modelos com interações três a três e com a interação quatro a quatro tiveram problemas de convergência para o método stepwise e por isso não foram considerados.

##Stepwise selection
step(glm(formula = Y ~ S+R+ESC+I+S*R+S*ESC+S*I+R*ESC+R*I+ESC*I, family=binomial(link=logit), data = dados), direction = "both")

## Start:  AIC=6785.71
## Y ~ S + R + ESC + I + S * R + S * ESC + S * I + R * ESC + R * 
##     I + ESC * I
## 
##         Df Deviance    AIC
## - R:ESC 10   6723.3 6773.3
## - S:R    3   6717.8 6781.8
## - R:I    3   6720.5 6784.5
## - S:ESC  4   6722.6 6784.6
## - ESC:I  4   6722.6 6784.6
## - S:I    1   6716.9 6784.9
## <none>       6715.7 6785.7
## 
## Step:  AIC=6773.3
## Y ~ S + R + ESC + I + S:R + S:ESC + S:I + R:I + ESC:I
## 
##         Df Deviance    AIC
## - S:R    3   6726.1 6770.1
## - R:I    3   6727.2 6771.2
## - S:ESC  4   6729.7 6771.7
## - ESC:I  4   6730.1 6772.1
## - S:I    1   6724.5 6772.5
## <none>       6723.3 6773.3
## + R:ESC 10   6715.7 6785.7
## 
## Step:  AIC=6770.1
## Y ~ S + R + ESC + I + S:ESC + S:I + R:I + ESC:I
## 
##         Df Deviance    AIC
## - R:I    3   6730.1 6768.1
## - S:ESC  4   6732.3 6768.3
## - ESC:I  4   6733.1 6769.1
## - S:I    1   6727.2 6769.2
## <none>       6726.1 6770.1
## + S:R    3   6723.3 6773.3
## + R:ESC 10   6717.8 6781.8
## 
## Step:  AIC=6768.12
## Y ~ S + R + ESC + I + S:ESC + S:I + ESC:I
## 
##         Df Deviance    AIC
## - S:ESC  4   6736.4 6766.4
## - ESC:I  4   6736.6 6766.6
## - S:I    1   6731.2 6767.2
## <none>       6730.1 6768.1
## + R:I    3   6726.1 6770.1
## + S:R    3   6727.2 6771.2
## + R:ESC 10   6722.6 6780.6
## - R      3   6816.6 6848.6
## 
## Step:  AIC=6766.39
## Y ~ S + R + ESC + I + S:I + ESC:I
## 
##         Df Deviance    AIC
## - S:I    1   6737.0 6765.0
## - ESC:I  4   6743.7 6765.7
## <none>       6736.4 6766.4
## + S:ESC  4   6730.1 6768.1
## + R:I    3   6732.3 6768.3
## + S:R    3   6733.6 6769.6
## + R:ESC 10   6729.4 6779.4
## - R      3   6823.4 6847.4
## 
## Step:  AIC=6764.96
## Y ~ S + R + ESC + I + ESC:I
## 
##         Df Deviance    AIC
## - ESC:I  4   6744.2 6764.2
## <none>       6737.0 6765.0
## + S:I    1   6736.4 6766.4
## + R:I    3   6732.9 6766.9
## + S:ESC  4   6731.2 6767.2
## + S:R    3   6734.2 6768.2
## + R:ESC 10   6730.0 6778.0
## - S      1   6801.8 6827.8
## - R      3   6824.0 6846.0
## 
## Step:  AIC=6764.21
## Y ~ S + R + ESC + I
## 
##         Df Deviance    AIC
## <none>       6744.2 6764.2
## + ESC:I  4   6737.0 6765.0
## + S:I    1   6743.7 6765.7
## + S:ESC  4   6737.7 6765.7
## + R:I    3   6740.7 6766.7
## + S:R    3   6741.4 6767.4
## + R:ESC 10   6737.1 6777.1
## - I      1   6781.2 6799.2
## - S      1   6809.5 6827.5
## - R      3   6831.1 6845.1
## - ESC    4   6995.5 7007.5

## 
## Call:  glm(formula = Y ~ S + R + ESC + I, family = binomial(link = logit), 
##     data = dados)
## 
## Coefficients:
## (Intercept)           S2           R2           R4           R5  
##    -3.54342     -0.71483      0.64436      0.54040     -1.38578  
##        ESC2         ESC3         ESC4         ESC5            I  
##     0.33147      0.09528     -0.87388     -2.14116      0.22666  
## 
## Degrees of Freedom: 5355 Total (i.e. Null);  5346 Residual
## Null Deviance:       7306 
## Residual Deviance: 6744  AIC: 6764

step(glm(formula = Y ~ S+R+ESC+I+S*R+S*ESC+S*I+R*ESC+R*I+ESC*I, family=binomial(link=probit), data = dados), direction = "both")

## Start:  AIC=6786.08
## Y ~ S + R + ESC + I + S * R + S * ESC + S * I + R * ESC + R * 
##     I + ESC * I
## 
##         Df Deviance    AIC
## - R:ESC 10   6723.9 6773.9
## - S:R    3   6718.0 6782.0
## - R:I    3   6720.7 6784.7
## - ESC:I  4   6723.1 6785.1
## - S:I    1   6717.2 6785.2
## - S:ESC  4   6723.5 6785.5
## <none>       6716.1 6786.1
## 
## Step:  AIC=6773.86
## Y ~ S + R + ESC + I + S:R + S:ESC + S:I + R:I + ESC:I
## 
##         Df Deviance    AIC
## - S:R    3   6726.4 6770.4
## - R:I    3   6727.7 6771.7
## - ESC:I  4   6730.7 6772.7
## - S:ESC  4   6730.9 6772.9
## - S:I    1   6725.0 6773.0
## <none>       6723.9 6773.9
## + R:ESC 10   6716.1 6786.1
## 
## Step:  AIC=6770.45
## Y ~ S + R + ESC + I + S:ESC + S:I + R:I + ESC:I
## 
##         Df Deviance    AIC
## - R:I    3   6730.4 6768.4
## - S:ESC  4   6733.3 6769.3
## - ESC:I  4   6733.5 6769.5
## - S:I    1   6727.5 6769.5
## <none>       6726.4 6770.4
## + S:R    3   6723.9 6773.9
## + R:ESC 10   6718.0 6782.0
## 
## Step:  AIC=6768.39
## Y ~ S + R + ESC + I + S:ESC + S:I + ESC:I
## 
##         Df Deviance    AIC
## - ESC:I  4   6736.9 6766.9
## - S:ESC  4   6737.3 6767.3
## - S:I    1   6731.4 6767.4
## <none>       6730.4 6768.4
## + R:I    3   6726.4 6770.4
## + S:R    3   6727.7 6771.7
## + R:ESC 10   6722.5 6780.5
## - R      3   6816.5 6848.5
## 
## Step:  AIC=6766.94
## Y ~ S + R + ESC + I + S:ESC + S:I
## 
##         Df Deviance    AIC
## - S:I    1   6738.0 6766.0
## - S:ESC  4   6744.7 6766.7
## <none>       6736.9 6766.9
## + ESC:I  4   6730.4 6768.4
## + R:I    3   6733.5 6769.5
## + S:R    3   6734.1 6770.1
## + R:ESC 10   6729.0 6779.0
## - R      3   6822.9 6846.9
## 
## Step:  AIC=6765.99
## Y ~ S + R + ESC + I + S:ESC
## 
##         Df Deviance    AIC
## - S:ESC  4   6745.2 6765.2
## <none>       6738.0 6766.0
## + S:I    1   6736.9 6766.9
## + ESC:I  4   6731.4 6767.4
## + R:I    3   6734.5 6768.5
## + S:R    3   6735.2 6769.2
## + R:ESC 10   6730.1 6778.1
## - I      1   6774.5 6800.5
## - R      3   6824.2 6846.2
## 
## Step:  AIC=6765.24
## Y ~ S + R + ESC + I
## 
##         Df Deviance    AIC
## <none>       6745.2 6765.2
## + ESC:I  4   6737.8 6765.8
## + S:ESC  4   6738.0 6766.0
## + S:I    1   6744.7 6766.7
## + R:I    3   6741.7 6767.7
## + S:R    3   6742.5 6768.5
## + R:ESC 10   6737.7 6777.7
## - I      1   6782.1 6800.1
## - S      1   6810.2 6828.2
## - R      3   6831.8 6845.8
## - ESC    4   6995.1 7007.1

## 
## Call:  glm(formula = Y ~ S + R + ESC + I, family = binomial(link = probit), 
##     data = dados)
## 
## Coefficients:
## (Intercept)           S2           R2           R4           R5  
##    -2.15891     -0.43592      0.39290      0.33093     -0.83534  
##        ESC2         ESC3         ESC4         ESC5            I  
##     0.20255      0.05935     -0.53803     -1.25556      0.13830  
## 
## Degrees of Freedom: 5355 Total (i.e. Null);  5346 Residual
## Null Deviance:       7306 
## Residual Deviance: 6745  AIC: 6765

step(glm(formula = Y ~ S+R+ESC+I+S*R+S*ESC+S*I+R*ESC+R*I+ESC*I, family=binomial(link=cloglog), data = dados), direction = "both")

## Start:  AIC=6788.08
## Y ~ S + R + ESC + I + S * R + S * ESC + S * I + R * ESC + R * 
##     I + ESC * I
## 
##         Df Deviance    AIC
## - R:ESC 10   6724.5 6774.5
## - S:ESC  4   6722.3 6784.3
## - R:I    3   6721.3 6785.3
## - S:R    3   6721.3 6785.3
## - ESC:I  4   6724.7 6786.7
## - S:I    1   6718.7 6786.7
## <none>       6718.1 6788.1
## 
## Step:  AIC=6774.54
## Y ~ S + R + ESC + I + S:R + S:ESC + S:I + R:I + ESC:I
## 
##         Df Deviance    AIC
## - S:ESC  4   6728.3 6770.3
## - R:I    3   6727.4 6771.4
## - S:R    3   6728.6 6772.6
## - ESC:I  4   6731.1 6773.1
## - S:I    1   6725.2 6773.2
## <none>       6724.5 6774.5
## + R:ESC 10   6718.1 6788.1
## 
## Step:  AIC=6770.27
## Y ~ S + R + ESC + I + S:R + S:I + R:I + ESC:I
## 
##         Df Deviance    AIC
## - R:I    3   6731.2 6767.2
## - S:R    3   6732.1 6768.1
## - S:I    1   6728.7 6768.7
## - ESC:I  4   6735.5 6769.5
## <none>       6728.3 6770.3
## + S:ESC  4   6724.5 6774.5
## + R:ESC 10   6722.3 6784.3
## 
## Step:  AIC=6767.17
## Y ~ S + R + ESC + I + S:R + S:I + ESC:I
## 
##         Df Deviance    AIC
## - S:R    3   6735.0 6765.0
## - S:I    1   6731.6 6765.6
## - ESC:I  4   6738.0 6766.0
## <none>       6731.2 6767.2
## + R:I    3   6728.3 6770.3
## + S:ESC  4   6727.4 6771.4
## + R:ESC 10   6725.5 6781.5
## 
## Step:  AIC=6765.02
## Y ~ S + R + ESC + I + S:I + ESC:I
## 
##         Df Deviance    AIC
## - S:I    1   6735.4 6763.4
## - ESC:I  4   6741.9 6763.9
## <none>       6735.0 6765.0
## + S:R    3   6731.2 6767.2
## + R:I    3   6732.1 6768.1
## + S:ESC  4   6731.5 6769.5
## + R:ESC 10   6728.5 6778.5
## - R      3   6820.4 6844.4
## 
## Step:  AIC=6763.41
## Y ~ S + R + ESC + I + ESC:I
## 
##         Df Deviance    AIC
## - ESC:I  4   6742.2 6762.2
## <none>       6735.4 6763.4
## + S:I    1   6735.0 6765.0
## + S:R    3   6731.6 6765.6
## + R:I    3   6732.5 6766.5
## + S:ESC  4   6732.1 6768.1
## + R:ESC 10   6728.9 6776.9
## - S      1   6802.4 6828.4
## - R      3   6820.7 6842.7
## 
## Step:  AIC=6762.19
## Y ~ S + R + ESC + I
## 
##         Df Deviance    AIC
## <none>       6742.2 6762.2
## + ESC:I  4   6735.4 6763.4
## + S:I    1   6741.9 6763.9
## + S:R    3   6738.4 6764.4
## + R:I    3   6739.7 6765.7
## + S:ESC  4   6738.2 6766.2
## + R:ESC 10   6735.7 6775.7
## - I      1   6779.5 6797.5
## - S      1   6809.7 6827.7
## - R      3   6828.2 6842.2
## - ESC    4   6995.3 7007.3

## 
## Call:  glm(formula = Y ~ S + R + ESC + I, family = binomial(link = cloglog), 
##     data = dados)
## 
## Coefficients:
## (Intercept)           S2           R2           R4           R5  
##    -2.71434     -0.52265      0.41674      0.37029     -1.16139  
##        ESC2         ESC3         ESC4         ESC5            I  
##     0.22393      0.08965     -0.60256     -1.69845      0.14817  
## 
## Degrees of Freedom: 5355 Total (i.e. Null);  5346 Residual
## Null Deviance:       7306 
## Residual Deviance: 6742  AIC: 6762

step(glm(formula = Y ~ S+R+ESC+I+S*R+S*ESC+S*I+R*ESC+R*I+ESC*I, family=binomial(link=cauchit), data = dados), direction = "both")

## Start:  AIC=6782.69
## Y ~ S + R + ESC + I + S * R + S * ESC + S * I + R * ESC + R * 
##     I + ESC * I
## 
##         Df Deviance    AIC
## - R:ESC 10   6720.7 6770.7
## - S:ESC  4   6717.4 6779.4
## - S:R    3   6715.7 6779.7
## - ESC:I  4   6719.8 6781.8
## - S:I    1   6714.5 6782.5
## <none>       6712.7 6782.7
## - R:I    3   6719.2 6783.2
## 
## Step:  AIC=6770.66
## Y ~ S + R + ESC + I + S:R + S:ESC + S:I + R:I + ESC:I
## 
##         Df Deviance    AIC
## - S:ESC  4   6724.3 6766.3
## - S:R    3   6724.6 6768.6
## - ESC:I  4   6727.3 6769.3
## - R:I    3   6725.3 6769.3
## - S:I    1   6722.4 6770.4
## <none>       6720.7 6770.7
## + R:ESC 10   6712.7 6782.7
## 
## Step:  AIC=6766.32
## Y ~ S + R + ESC + I + S:R + S:I + R:I + ESC:I
## 
##         Df Deviance    AIC
## - S:R    3   6728.3 6764.3
## - R:I    3   6729.0 6765.0
## - ESC:I  4   6731.5 6765.5
## - S:I    1   6725.5 6765.5
## <none>       6724.3 6766.3
## + S:ESC  4   6720.7 6770.7
## + R:ESC 10   6717.4 6779.4
## 
## Step:  AIC=6764.33
## Y ~ S + R + ESC + I + S:I + R:I + ESC:I
## 
##         Df Deviance    AIC
## - R:I    3   6733.0 6763.0
## - ESC:I  4   6735.3 6763.3
## - S:I    1   6729.5 6763.5
## <none>       6728.3 6764.3
## + S:R    3   6724.3 6766.3
## + S:ESC  4   6724.6 6768.6
## + R:ESC 10   6720.4 6776.4
## 
## Step:  AIC=6762.99
## Y ~ S + R + ESC + I + S:I + ESC:I
## 
##         Df Deviance    AIC
## - ESC:I  4   6739.6 6761.6
## - S:I    1   6734.0 6762.0
## <none>       6733.0 6763.0
## + R:I    3   6728.3 6764.3
## + S:R    3   6729.0 6765.0
## + S:ESC  4   6729.3 6767.3
## + R:ESC 10   6726.8 6776.8
## - R      3   6821.1 6845.1
## 
## Step:  AIC=6761.58
## Y ~ S + R + ESC + I + S:I
## 
##         Df Deviance    AIC
## - S:I    1   6740.6 6760.6
## <none>       6739.6 6761.6
## + ESC:I  4   6733.0 6763.0
## + R:I    3   6735.3 6763.3
## + S:R    3   6735.7 6763.7
## + S:ESC  4   6735.4 6765.4
## + R:ESC 10   6733.8 6775.8
## - R      3   6827.0 6843.0
## - ESC    4   6997.0 7011.0
## 
## Step:  AIC=6760.58
## Y ~ S + R + ESC + I
## 
##         Df Deviance    AIC
## <none>       6740.6 6760.6
## + S:I    1   6739.6 6761.6
## + ESC:I  4   6734.0 6762.0
## + R:I    3   6736.3 6762.3
## + S:R    3   6736.6 6762.6
## + S:ESC  4   6736.9 6764.9
## + R:ESC 10   6734.7 6774.7
## - I      1   6777.8 6795.8
## - S      1   6806.0 6824.0
## - R      3   6827.9 6841.9
## - ESC    4   6998.0 7010.0

## 
## Call:  glm(formula = Y ~ S + R + ESC + I, family = binomial(link = cauchit), 
##     data = dados)
## 
## Coefficients:
## (Intercept)           S2           R2           R4           R5  
##    -3.26795     -0.64598      0.59745      0.48396     -1.44031  
##        ESC2         ESC3         ESC4         ESC5            I  
##     0.29408      0.06591     -0.77829     -2.53332      0.20849  
## 
## Degrees of Freedom: 5355 Total (i.e. Null);  5346 Residual
## Null Deviance:       7306 
## Residual Deviance: 6741  AIC: 6761

O algoritmo com diferentes funções de ligação resultou no mesmo conjunto de variáveis, são elas: idade do indivíduo, raça/cor, escolaridade e sexo. O resultado do algoritmo é coerente com os gráficos apresentados anteriormente, pois todas as variáveis diferem dentre suas categorias em relação a discriminação da morte. Portanto, no decorrer deste trabalho a análise seguirá com o modelo abaixo.

\[\begin{eqnarray} Y_i &\sim& Binomial (n, \widehat{\pi_i})\\ g(\widehat{\pi_i}) &=& \widehat{\beta_{0}} + \widehat{\beta_{1}}sexo_{i} + \widehat{\beta_{2}}raça_{i}+\widehat{\beta_{3}}Escolaridade_{i}+\widehat{\beta_{4}}Idade_{i} \end{eqnarray}\]

Sendo as categorias determinadas conforme tabela 1 apresentada anteriormente.

Outras formas de escolha do melhor modelo são:

• Usando o pacote glmulti

library (glmulti)
search.1.aicc<-glmulti(y=Y ~.,data=dados, fitfunction="glm",level=1,method="h",crit="aicc",
                       family=binomial(link="logit"))

## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
## Completed.

print(search.1.aicc)

## glmulti.analysis
## Method: h / Fitting: glm / IC used: aicc
## Level: 1 / Marginality: FALSE
## From 16 models:
## Best IC: 6764.24890015859
## Best model:
## [1] "Y ~ 1 + S + R + ESC + I"
## Evidence weight: 0.99999997453726
## Worst IC: 7308.25373637169
## 1 models within 2 IC units.
## 0 models to reach 95% of evidence weight.

aa<-weightable(search.1.aicc)
cbind(model=aa[1:5,1],round(aa[1:5,2:3],digits=3))

##                     model     aicc weights
## 1 Y ~ 1 + S + R + ESC + I 6764.249       1
## 2     Y ~ 1 + S + R + ESC 6799.221       0
## 3     Y ~ 1 + R + ESC + I 6827.530       0
## 4     Y ~ 1 + S + ESC + I 6845.074       0
## 5         Y ~ 1 + R + ESC 6871.225       0

Note que não houve diferença em relação ao modelo definido.

Com os componentes aleatório e sistemático do modelo já definidos conforme discussões anteriores, é necessário fazer a escolha da função de ligação. Dentre as funções de ligação definidas para uma variável resposta binária, faz-se-a um comparativo conforme medidas descritas nos gráficos e tabela abaixo.

source("misc-fun.R")
model.log <- glm(formula = Y ~ S+R+ESC+I,
                 family = binomial(link = "logit"), data = dados)
dev1<-model.log$null.deviance
dev2<-model.log$deviance
r2_log<-1-(dev2/dev1)
r2_log

## [1] 0.07692661

model.prob <- glm(formula = Y ~ S+R+ESC+I,
                  family = binomial(link = "probit"), data = dados)
dev1<-model.prob$null.deviance
dev2<-model.prob$deviance
r2_prob<-1-(dev2/dev1)
r2_prob

## [1] 0.07678509

model.clog <- glm(formula = Y ~ S+R+ESC+I,
                  family = binomial(link = "cloglog"), data = dados)
dev1<-model.clog$null.deviance
dev2<-model.clog$deviance
r2_clog<-1-(dev2/dev1)
r2_clog

## [1] 0.0772023

model.cauc <- glm(formula = Y ~ S+R+ESC+I,
                  family = binomial(link = "cauchit"), data = dados)
dev1<-model.cauc$null.deviance
dev2<-model.cauc$deviance
r2_cauc<-1-(dev2/dev1)
r2_cauc

## [1] 0.07742312

par(mfrow=c(1,4))
envelope(model.log)
title("Ligação Logit")
envelope(model.prob)
title("Ligação Probit")
envelope(model.clog)
title("Ligação Complemento Log-Log")
envelope(model.cauc)
title("Ligação Cauchit")

Para verificar a qualidade do modelo ajustado, foi realizada a análise gráfica da curva ROC, como mostra as figuras abaixo.

library(pROC)
roc(model.log$y,model.log$fitted.values,plot=T,ci=T,identity=TRUE,print.auc=TRUE)

## 
## Call:
## roc.default(response = model.log$y, predictor = model.log$fitted.values,     ci = T, plot = T, identity = TRUE, print.auc = TRUE)
## 
## Data: model.log$fitted.values in 2280 controls (model.log$y 0) < 3076 cases (model.log$y 1).
## Area under the curve: 0.6802
## 95% CI: 0.6657-0.6947 (DeLong)

roc(model.prob$y,model.prob$fitted.values,plot=T,ci=T,identity=TRUE,print.auc=TRUE)

## 
## Call:
## roc.default(response = model.prob$y, predictor = model.prob$fitted.values,     ci = T, plot = T, identity = TRUE, print.auc = TRUE)
## 
## Data: model.prob$fitted.values in 2280 controls (model.prob$y 0) < 3076 cases (model.prob$y 1).
## Area under the curve: 0.6802
## 95% CI: 0.6657-0.6947 (DeLong)

roc(model.clog$y,model.clog$fitted.values,plot=T,ci=T,identity=TRUE,print.auc=TRUE)

## 
## Call:
## roc.default(response = model.clog$y, predictor = model.clog$fitted.values,     ci = T, plot = T, identity = TRUE, print.auc = TRUE)
## 
## Data: model.clog$fitted.values in 2280 controls (model.clog$y 0) < 3076 cases (model.clog$y 1).
## Area under the curve: 0.6796
## 95% CI: 0.6651-0.6941 (DeLong)

roc(model.cauc$y,model.cauc$fitted.values,plot=T,ci=T,identity=TRUE,print.auc=TRUE)

## 
## Call:
## roc.default(response = model.cauc$y, predictor = model.cauc$fitted.values,     ci = T, plot = T, identity = TRUE, print.auc = TRUE)
## 
## Data: model.cauc$fitted.values in 2280 controls (model.cauc$y 0) < 3076 cases (model.cauc$y 1).
## Area under the curve: 0.6802
## 95% CI: 0.6657-0.6947 (DeLong)

O valor correspondente à área abaixo da curva ROC semelhante para todos os modelos, e de acordo com os níveis conhecidos indicam que o modelo detém uma capacidade de discriminação aceitável.

\[\begin{array}{lccccc} \hline \textrm{Ligação} & \textrm{gl} & \textrm{AIC} & \textrm{Deviance} & \textrm{Area ROC} & \textrm{Pseudo R^2} \\ \hline \textrm{Logit }& 5355 & 6764 & 6744 & 0.6802 & 0.0769 \\ \textrm{Probit }& 5355 & 6765 & 6745 & 0.6802 & 0.0768 \\ \textrm{Comp. log-log }& 5355 & 6762 & 6742 & 0.6796 & 0.0772 \\ \textrm{Cauchit }& 5355 & 6761 & 6741 & 0.6802 & 0.0774 \\ \hline \end{array}\]

O teste de Hosmer-Lemeshow é um teste que avalia o modelo ajustado comparando as frequências observadas e as esperadas. O teste associa os dados as suas probabilidades estimadas da mais baixa a mais alta, então faz um teste qui-quadrado para determinar se as frequências observadas estão próximas das frequências esperadas. Já o teste de Pearson, é utilizado para fazer análise dos resíduos para modelos logísticos, trata-se de uma medida útil para avaliar o quão bem o modelo selecionado ajustou-se aos dados.

As hipóteses testadas foram:

• Ho: O modelo é adequado
• H1: O modelo não é adequado

library(ResourceSelection)
hoslem.test(model.log$y, fitted(model.log))

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model.log$y, fitted(model.log)
## X-squared = 13.027, df = 8, p-value = 0.1109

hoslem.test(model.prob$y, fitted(model.prob))

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model.prob$y, fitted(model.prob)
## X-squared = 13.456, df = 8, p-value = 0.09711

hoslem.test(model.clog$y, fitted(model.clog))

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model.clog$y, fitted(model.clog)
## X-squared = 16.935, df = 8, p-value = 0.0308

hoslem.test(model.cauc$y, fitted(model.cauc))

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model.cauc$y, fitted(model.cauc)
## X-squared = 8.7688, df = 8, p-value = 0.3622

• P-valor do teste qui-quadrado de Pearson:

#Resíduo de Pearson
rp=residuals(model.log,type="pearson")
1-pchisq(sum(rp^2),model.log$df.residu)

## [1] 0.4151227

rp=residuals(model.prob,type="pearson")
1-pchisq(sum(rp^2),model.prob$df.residu)

## [1] 0.4165986

rp=residuals(model.clog,type="pearson")
1-pchisq(sum(rp^2),model.clog$df.residu)

## [1] 0.4382763

rp=residuals(model.cauc,type="pearson")
1-pchisq(sum(rp^2),model.cauc$df.residu)

## [1] 0.4344812

Os modelos especificados com diferentes funções de ligação apresentaram um comportamento muito parecido. Com base nos gráficos de resíduos simulados, não há problemas quanto a especificação do modelo nas quatro diferentes funções de ligação propostas, todos os gráficos apresentaram resíduos dentro dos intervalos simulados. Já com base nas medidas de comparação exibidas na tabela, não nota-se desempenho distinto, em relação ao teste de Hosmer-Lemeshow somente o modelo de Clog-log não é adequado e de acordo com o resíduo de Pearson, todos são adequados, logo, em função da magnitude das medidas comparativas e pela vantagem interpretativa da especificação logit, dada em função de razão de chances, esta foi definida no modelo proposto.

4 Modelo ajustado

Com as análises realizadas, o modelo é definido da seguinte forma:

\[\begin{equation*} \begin{gathered} Y_i &\sim& Binomial (n, \widehat{\pi_i})\\ g(\widehat{\pi_i}) &=& \widehat{\beta_{0}} + \widehat{\beta_{1}}sexo_{i} + \widehat{\beta_{2}}raça_{i}+\widehat{\beta_{3}}Escolaridade_{i}+\widehat{\beta_{4}}Idade_{i} \end{gathered} \end{equation*}\]

A saída do R abaixo apresenta o conjunto de variáveis que compõem o modelo ajustado junto com os valores estimados dos coeficientes do modelo, o erro padrão dos coeficientes, os quantis Z da distribuição normal padrão e o p-valor de todos os parâmetros:

bestfit<- glm(formula = Y ~ S+R+ESC+I,
                 family = binomial(link = "logit"), data = dados)
summary(bestfit)

## 
## Call:
## glm(formula = Y ~ S + R + ESC + I, family = binomial(link = "logit"), 
##     data = dados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7491  -1.1673   0.8069   0.9369   2.2899  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.54342    0.66772  -5.307 1.12e-07 ***
## S2          -0.71483    0.08940  -7.996 1.28e-15 ***
## R2           0.64436    0.12215   5.275 1.33e-07 ***
## R4           0.54040    0.06591   8.199 2.42e-16 ***
## R5          -1.38578    0.55451  -2.499  0.01245 *  
## ESC2         0.33147    0.27559   1.203  0.22906    
## ESC3         0.09528    0.26815   0.355  0.72235    
## ESC4        -0.87388    0.27176  -3.216  0.00130 ** 
## ESC5        -2.14116    0.68340  -3.133  0.00173 ** 
## I            0.22666    0.03733   6.071 1.27e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 7306.3  on 5355  degrees of freedom
## Residual deviance: 6744.2  on 5346  degrees of freedom
## AIC: 6764.2
## 
## Number of Fisher Scoring iterations: 4

As categorias ESC2 e ESC3 foram não significativas, neste caso foi realizado o agrupamento das duas e avaliado um novo ajuste com a covariável escolaridade com somente três categorias, porém a categoria nova permaneceu não significativa, dessa forma, preferiu-se manter a covariável com quatro categorias, uma vez que não houve alteração na significância dos parâmetros.

5 Análise de Diagnóstico

Para subsidiar a avaliação da qualidade do modelo, a análise de diagnóstico verificará a adequação da distribuição proposta, da função de ligação, e do preditor linear. Ou seja, do modelo de regressão ajustado aos dados.

#Com relação a especificação das covariáveis no modelo podemos observar
#na figura que não há grandes evidências de má especificação das
#covariáveis no modelo, mesmo sendo observada a assimetria dos
#resíduos. Note que a interpretação dos gráficos de diagnóstico é mais
#flexível nestes casos, pois a limitação da variável resposta (suporte 0
#ou 1) interfere na interpretação gráfica.


## Preditor Linear
#a <- car::residualPlots(bestfit, layout=c(1,3))

Abaixo são apresentados 3 gráficos que auxiliam na identificação de possíveis fuga de suposições do modelo. No caso apresentado, não temos evidências gráficas para suspeitar de nenhuma suposição não atendida. No gráfico 1 a magnitude dos resíduos ultrapassa 2 somente para alguns resíduos. No segundo gráfico representando o resíduo vs. valores ajustados, temos uma disposição aparentemente centrada em zero, lembre-se que a natureza da variável resposta dificulta a interpretação. No terceiro e último gráfico deste figura, (c), temos o gráfico quantil-quantil com envelope simulado, onde resíduos dispostos dentro das bandas de confiança representam adequação dos dados ao modelo proposto.

## Pressuposições
diag.binom(bestfit)

Como entre as covariáveis, somente a idade é continua, porém com três valores distintos somente, 15, 16 e 17 anos, não foi verificado a existência de pontos influentes, pois caso, um ponto seja influente, vários também como mesma discriminação será.

6 Resultados

Com o modelo especificado e avaliado podemos realizar predições e interpretações. Utilizando as estimativas dos parâmetros, podemos encontrar os valores da Odds Ratio do modelo e os respectivos intervalos de confiança dos parâmetros do modelo. As categorias de referência foram, sexo masculino, Raça/Cor branca e nível de escolaridade nenhum estudo.

\[\begin{array}{ccccc} \hline \textrm{Categorias }& \textrm{Parâmetros }& \textrm{OR }& 2.5 \% & 97.5 \% \\ \hline \textrm{(Intercepto)}& -3.54 & 0.03 & -4.85 & -2.23 \\ \textrm{ S2 }& -0.71 & 0.49 & -0.89 & -0.54 \\ \textrm{ R2 }& 0.64 & 1.90 & 0.41 & 0.89 \\ \textrm{ R4 }& 0.54 & 1.72 & 0.41 & 0.67 \\ \textrm{ R5 }& -1.39 & 0.25 & -2.63 & -0.40 \\ \textrm{ ESC2 }& 0.33 & 1.39 & -0.22 & 0.86 \\ \textrm{ ESC3 }& 0.10 & 1.10 & -0.44 & 0.61 \\ \textrm{ ESC4 }& -0.87 & 0.42 & -1.42 & -0.35 \\ \textrm{ ESC5 }& -2.14 & 0.12 & -3.68 & -0.92 \\ \textrm{ I }& 0.23 & 1.25 & 0.15 & 0.30 \\ \hline \end{array}\]

Após o ajuste do modelo, pode-se usá-lo para estimar a probabilidade de um indivíduo que morreu, ter sido assassinado, a partir de:

\[\begin{eqnarray} \hat{\pi} = \frac{\exp(-3.54-0.72S2+0.64R2+\cdots-0.87ESC4-2.14ESC5+0.23I)}{1+exp(-3.54-0.72S2+0.64R2+\cdots-0.87ESC4-2.14ESC5+0.23I)}, \end{eqnarray}\]

Na tabela abaixo apresenta-se a probabilidade de um jovem que morreu, ter sido assassinado, dado algumas características. Note que, dado que está morto, um jovem de 17 anos, solteiro, do sexo masculino, de cor negra e com nenhuma escolaridade possui 72,19% de probabilidade de ter morrido assassinado.

\[\begin{array}{ccccc} \hline \textrm{ Sexo} &\textrm{ Raça/Cor}&\textrm{ Escolaridade} &\textrm{ Idade} & \textrm{Probabilidade} \\ \hline \textrm{Masculino }&\textrm{ Branca }&\textrm{ Nenhuma }& 17 & 0,58 \\ \textrm{ Masculino}&\textrm{ Negra }&\textrm{ Nenhum }& 17 & 0,72 \\ \textrm{ Masculino}&\textrm{ Negra }&\textrm{ 12 ou mais anos }& 17 & 0,23 \\ \textrm{ Feminina }&\textrm{ Branca }&\textrm{ 8 a 11 anos }& 16 & 0,18 \\ \textrm{ Feminina }&\textrm{ Branca }&\textrm{ 8 a 11 anos }& 17 & 0,22 \\ \textrm{ Feminina }&\textrm{ Parda }&\textrm{ 12 ou mais anos }& 17 & 0,12 \\ \hline \end{array}\]

Pode-se verificar ainda que:

• A chance de um indivíduo negro ser assassinado é 90% maior que a chance de um indivíduo branco;
• A chance de um indivíduo pardo ser assassinado é 72% maior que a chance de um indivíduo branco;
• A chance de um indivíduo indígena ser assassinado é 75% menor que a chance de um indivíduo branco;
• A chance de um indivíduo que tem de 8 a 11 anos de estudo ser assassinado é 58% menor que a chance de um indivíduo com nenhum estudo;
• A chance de um indivíduo que tem 12 ou mais anos de estudo ser assassinado é 88% menor que a chance de um indivíduo com nenhum estudo;
• A chance de um indivíduo do sexo feminino ser assassinado é 52% menor que a chance de um indivíduo do sexo masculino;
• A mudança de uma unidade na idade altera em 25% o logito do modelo.

7 Exemplo de poster deste trabalho em PDF e tex

• Poster em pdf

clique aqui

clique aqui

• Poster em tex

clique aqui

Referências

“Homicídio é Principal Causa de Mortes de Jovens de 16 E 17 No País.” n.d. http://g1.globo.com/politica/noticia/2015/06/homicidio-e-principal-causa-de-mortes-de-jovens-de-16-e-17-no-pais.html.

Trindade, Ruth França Cizino da, Flávia Azevedo de Mattos Moura Costa, Gustavo Bussi Caminiti, Claudia Benedita dos Santos, and others. 2015. “Mapa Dos Homicídios Por Arma de Fogo: Perfil Das Vítimas E Das Agressões.” Revista Da Escola de Enfermagem Da USP 49 (5): 748–55.

Waiselfisz, Julio Jacobo. 2015. “Mapa Da Violência 2015: Adolescentes de 16 E 17 Anos Do Brasil.” Centro Brasileiro de Estudos Latino-Americanos/Faculdade Latino-Americana de Ciências Sociais Brasil.