Aula 6: Script modelo logı́stico multinı́vel
Leo Bastos
9 de junho de 2015
O arquivo ’obitosCausaExternasRJ2012.csv’ contém o total de óbitos no
ano de 2012 por causas externas nos municı́pios do Rio de Janeiro por local de
residência, separados por sexo e raça.
Os dados sao provenientes do DATASUS. O desfecho é o óbito por causas
externa por sexo. A hipótese principal é de que os óbitos entre homens é maior
que das mulheres.
obitosRJ <- read.csv("obitosCausaExternasRJ2012.csv", na.strings = "-")
str(obitosRJ)
Note que esse dado é agregado, podemos ter os dados para uma análise de regressão logı́stica em diferentes formatos. Outro formato, é o formato longo, uma
linha por observação. Vou apresentar ambos formatos e ajustar a modelagem
de forma adequada.
O primeiro formato é o formato fornecido pelo DATASUS, que seria uma
observação por municı́pio (n = 92). Como eu criei quatro categorias de raça
então deverı́amos ter n = 92x4 = 368 linhas. Como houveram municı́pios que
não registraram mortes por causas externas nem de homens nem de mulheres
no ano de 2012 para alguma raça, então o número de linhas do banco se reduziu
para 304 linhas.
# Tratamento
# Definindo a raça branca como base.
obitosRJ$Raca <- relevel( obitosRJ$Raca, ref = "Branca")
# Assumindo que NA significa a n~
ao ocorr^
encia de
# casos naquele municı́pio.
obitosRJ$Masc[which(is.na(obitosRJ$Masc))] <- 0
obitosRJ$Fem[which(is.na(obitosRJ$Fem))] <- 0
# O modelo de regress~
ao logı́stica para esse formato de dados é dado por
glm( cbind(Masc, Fem) ~ 1, data = obitosRJ, family=binomial() )
glm( cbind(Masc, Fem) ~ 1 + Raca, data = obitosRJ, family=binomial() )
1
O segundo formato de dados, seria o formato longo, uma observação por
indivı́duo. O banco ’obitosCausaExternasRJ2012long.csv’ contém esses dados
desagregados, e é consideravelmente maior. Tem 12355 linhas, e três variáveis
(raça, municı́pio e sexo).
obitosRJ.long <- read.csv("obitosCausaExternasRJ2012long.csv")
# Y = 1 -> Masculino, Y = 0 -> Feminino.
str(obitosRJ.long)
# Mudando a base da raça
obitosRJ.long$Raca <- relevel( obitosRJ.long$Raca, ref = "Branca")
# O modelo de regress~
ao logı́stica para esse formato de dados é dado por
glm( Y ~ 1, data = obitosRJ.long, family=binomial() )
glm( Y ~ 1 + Raca, data = obitosRJ.long, family=binomial() )
Note que os coeficientes estimados são idênticos para a modelagem nos dois
formatos. Qual usar? Tanto faz, desde que use o comando apropriado.
Modelo multinı́vel
É razoável pensar que a mortalidade por causas externas por sexo seja diferente
em municı́pios distintos. Logo, vamos incluir no modelo os municı́pios como um
nı́vel hierárquico na modelagem.
require(lme4)
# Para os dados agregados
glmer(
glmer(
cbind(Masc, Fem) ~ (1|Munic.pio), data = obitosRJ, family=binomial() )
cbind(Masc, Fem) ~ Raca + (1|Munic.pio), data = obitosRJ, family=binomial() )
# Para os dados desagregados
glmer( Y ~ (1|Munic.pio), data = obitosRJ.long, family=binomial() )
glmer( Y ~ Raca + (1|Munic.pio), data = obitosRJ.long, family=binomial() )
Note novamente que os coeficientes são os mesmos, mas o tempo computacional é diferente para os dados agregados e desagregados.
Para fazer essa análise, vamos usar os dados agregados, e ajustar o modelo
logı́stico tradicional, com e sem raça, e o modelo logı́stico multinı́vel, com e sem
raça. E comparar os modelos usando um critério de comparação, no caso o AIC.
2
m0
m1
m2
m3
<<<<-
glm( cbind(Masc, Fem) ~ 1, data = obitosRJ, family=binomial() )
glm( cbind(Masc, Fem) ~ 1 + Raca, data = obitosRJ, family=binomial() )
glmer( cbind(Masc, Fem) ~ (1|Munic.pio), data = obitosRJ, family=binomial() )
glmer( cbind(Masc, Fem) ~ Raca + (1|Munic.pio), data = obitosRJ, family=binomial() )
extractAIC(m0)
extractAIC(m1)
extractAIC(m2)
extractAIC(m3)
Vamos escolher o modelo multinı́vel logı́stico com raça e olhar os coeficientes
para responder a pergunta da hipótese.
require(arm)
m3.ranef <- ranef(m3)
m3.ranef.se <- se.ranef(m3)
m3.coef <- m3.ranef$Munic.pio$'(Intercept)'
m3.coef.se <- m3.ranef.se$Munic.pio[,1]
coefplot(m3.coef, m3.coef.se)
View(cbind(m3.ranef$Munic.pio, m3.coef.se))
3