Fatores de Risco Sociais que Impactam em Suicídio Familiar - (2018)

1. Base de Dados e Pré-processamento

1.1 Carregando a base de dados

#carregando a base de dados em .csv
dados <- 
  data.table::fread("/Dados/Coorte100m/PBF_SHSM/DB_ORIGINAL/nih_sm/02_nih_sm_20220627_baseline_pop100_v3_base_sih_multimorbidade_base_sim.csv", select = c("cod_familiar_fam_eq","cod_sexo_pessoa_eq","cod_raca_cor_pessoa_eq","escolaridade_eq", "cod_local_domic_fam_eq","qtd_dias_receb_fam","cod_destino_lixo_domic_fam_eq", "cod_escoa_sanitario_domic_fam_eq","cod_abaste_agua_domic_fam_eq","cod_munic_ibge_fam_eq","causabas_sim", "dt_entrada_coorte", "dt_nasc_pessoa_eq")) 

|--------------------------------------------------|
|==================================================|

1.2 Pre-processando a base original

#funcao para calcular a moda
moda <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

#criando a variavel idade
dados <- dados %>% 
  mutate( 
    idade = as.integer((dt_entrada_coorte - dt_nasc_pessoa_eq) / 365)
  )

#retirando individuos com menos de 10 anos

dados <- dados %>% filter(idade > 10) #retirando criancas menores de 10 anos


#criando a variavel regiao

dados <- dados %>% 
  mutate(
      regiao = as.numeric(substr(cod_munic_ibge_fam_eq, 1, 1)), #selecionando o pirmeiro digito
      regiao = case_when( #rotulando numero por regiao
        regiao == 1~"Norte",
        regiao == 2~"Nordeste",
        regiao == 3~"Sudeste",
        regiao == 4~"Sul",
        regiao == 5~"Centro-Oeste"),
      
      regiao = case_when( #deixando a regiao nordeste como referencia (numero 1)
        regiao == "Norte"~1,
        regiao == "Nordeste"~2,
        regiao == "Sudeste"~3,
        regiao == "Sul"~5,
        regiao == "Centro-Oeste"~4)
  )

#########################################################################################

#criando as variaveis suicide e homicide
dados <- dados %>% 
  mutate( #separando a variavel causabas_sim em duas partes
    letra = as.character(substr(causabas_sim, 1, 1)), #letra
    numero = as.character(substr(causabas_sim, 2 ,4)) #numero com 3 digitos
)

#transformando a segunda parte em numero inteiro
dados$numero <- as.integer(dados$numero)

Warning: NAs introduced by coercion

#suicide
dados <- dados %>% 
  mutate(
    suicide = case_when( #verificando a letra e os range numerico
      letra == "X" & (numero >= 600 & numero <= 849)~1,
      letra == "X" & (numero >= 60 & numero <= 84)~1,
      TRUE~0
    )
  )

#homicide
dados <- dados %>% 
  mutate(
    homicide = case_when( #verificando a letra e os range numerico
      letra == "X" & (numero >= 850 & numero <= 999)~1,
      letra == "Y" & (numero >= 000 & numero <= 099)~1,
      
      letra == "X" & (numero >= 85 & numero <= 99)~1,
      letra == "Y" & (numero >= 00 & numero <= 09)~1,
      TRUE~0
    )
  )

#########################################################################################

#criando a variavel (urban) urbano ou rural
dados <- dados %>% 
  mutate( 
    urban = case_when( 
      cod_local_domic_fam_eq == 1~1,
      cod_local_domic_fam_eq == 2~0,
      TRUE~0)
  )


#criando a variavel raca
dados <- dados %>% 
  mutate( 
    raca = case_when( 
      cod_raca_cor_pessoa_eq == 0~NA_character_,
      cod_raca_cor_pessoa_eq == 1~"1",
      cod_raca_cor_pessoa_eq == 2~"2",
      cod_raca_cor_pessoa_eq == 3~"3",
      cod_raca_cor_pessoa_eq == 4~"4",
      cod_raca_cor_pessoa_eq == 5~"5",
      cod_raca_cor_pessoa_eq == 99~NA_character_)
  )

#categorizando a variavel idade
dados <- dados %>% 
  mutate(
    idade = case_when(
      idade < 19~1,
      idade >= 20 & idade < 40~2,
      idade >= 40 & idade < 60~3,
      idade >= 60~4)
  )


#criando a variavel homem_familia
dados <- dados %>% 
  mutate(
    homem_familia = case_when(
      cod_sexo_pessoa_eq == 1~1,
      cod_sexo_pessoa_eq == 2~0)
  )


#criando a variavel escolaridade
dados <- dados %>% 
  mutate(
    escolaridade = case_when(
      escolaridade_eq == 0~0,
      escolaridade_eq == 1~1,
      escolaridade_eq == 2~2,
      escolaridade_eq == 3~3,
      escolaridade_eq == 4~4,
      escolaridade_eq == 5 | escolaridade_eq == 6~5)
  )


#criando a variavel LivingConditions
dados <- dados %>% 
  mutate( #criando a variavel lixo
    lixo = case_when(
      cod_destino_lixo_domic_fam_eq == 1 | cod_destino_lixo_domic_fam_eq == 2~1,
      TRUE~0),
    
     escoamento = case_when( #criando a variavel escoamento
       
      cod_escoa_sanitario_domic_fam_eq == 1 | cod_escoa_sanitario_domic_fam_eq == 2 | cod_escoa_sanitario_domic_fam_eq == 3~1,
      TRUE~0),
    
    agua = case_when( #criando a variavel agua
      cod_abaste_agua_domic_fam_eq == 1~1,
      TRUE~0),
    
    #criando a variavel living_conditions
    living_conditions = (lixo + escoamento + agua)
    
  )

1.3 Realizando a mudanca de grao da base

#familias com pelo menos um suicidio
qtd_suic_familia <- dados %>% 
  group_by(cod_familiar_fam_eq) %>% #agrupando por familias
  summarise(
    qtd_pesoas_familia = n(),#quantidade de pessoas na familia
    qtd_suic_familia = sum(suicide), #quantidade de suicidios na familia
    qtd_homic_familia = sum(homicide), #quantidade de homicidio na familia
    homem_familia = (sum(homem_familia)/n())*100, #percentual de homens na familia
    raca = moda(raca), #raca predominante(moda) na familia
    regiao = moda(regiao), #regiao predominante (moda) da familia
    urban = moda(urban), #predominancia de membros da familia em area urbana
    idade = moda(idade), #predominancia (moda) da faixa etaria da familia
    escolaridade = moda(escolaridade), #predominancia (moda) da faixa etaria da familia
    living_conditions = moda(living_conditions)#condicao de vida predominante(moda) na familia
    ) %>% 
  filter(qtd_pesoas_familia > 1) #familias com mais de uma pessoa


#categorizar a quantidade de pessoas na familia
qtd_suic_familia <- qtd_suic_familia %>% 
  mutate( 
    pessoas_na_familia = case_when(
      qtd_pesoas_familia == 2~"1",
      qtd_pesoas_familia == 3~"2",
      qtd_pesoas_familia == 4~"2",
      qtd_pesoas_familia == 5~"3",
      qtd_pesoas_familia == 6~"3",
      qtd_pesoas_familia >7~"4")
  )
  


#criando a variavel desfecho
dados_familia <- qtd_suic_familia %>% 
  mutate( #se teve pelo menos um suicidio (1) e caso contrrario (0)
    desfecho = case_when(qtd_suic_familia >= 1~1,TRUE~0)
  )

1.4 Salvando a base agregada no grao familia

#salvando a base de dados em .csv
rm(dados)
write.csv(dados_familia,"/Dados/Coorte100m/PBF_SHSM/SHARED/Suicidio_Familia_2018.csv")
rm(dados_familia)
rm(qtd_suic_familia)
gc()

            used   (Mb) gc trigger    (Mb)   max used    (Mb)
Ncells   2646635  141.4  654226649 34939.5  817783311 43674.4
Vcells 226465661 1727.8 4665603180 35595.8 5814904718 44364.3

2. Análise Exploratória dos Dados

2.1 Carregando a base de dados agregada

dados <- 
  data.table::fread("/Dados/Coorte100m/PBF_SHSM/SHARED/Suicidio_Familia_2018.csv") 

#binarizando homicidio
dados <- dados %>% 
  mutate( 
    homicide = case_when(
      qtd_homic_familia < 1~"0",
      qtd_homic_familia >= 1~"1")
  )

#transformando em fator
#dados$desfecho <- factor(dados$desfecho) 
dados$raca <- factor(dados$raca) 
dados$regiao <- factor(dados$regiao) 
dados$urban <- factor(dados$urban) 
dados$idade <- factor(dados$idade) 
dados$escolaridade <- factor(dados$escolaridade) 
dados$living_conditions <- factor(dados$living_conditions) 
dados$pessoas_na_familia <- factor(dados$pessoas_na_familia) 
dados$homicide <- factor(dados$homicide) 

2.2 Tabela de contingência vs Variável binária suicídio familiar

#criando a tabela de contingência
dados %>% 
  select( #aqui ta selecionando as variáveis na tabela
    homicide,
    pessoas_na_familia,
    homem_familia,
    raca,
    regiao,
    urban,
    idade,
    escolaridade,
    living_conditions,
    desfecho) %>% 
  tbl_summary(by = desfecho,
              statistic = list(all_continuous2() ~ "{mean} ({sd})")) %>% 
  add_p() %>% #adiciona os p-valores dos testes estatísticos 
   bold_labels() 

Characteristic	0, N = 23,547,2801	1, N = 44,8521	p-value2
homicide			<0.001
0	23,270,095 (99%)	43,768 (98%)
1	277,185 (1.2%)	1,084 (2.4%)
pessoas_na_familia			<0.001
1	12,156,205 (52%)	13,921 (32%)
2	8,847,805 (38%)	20,276 (47%)
3	2,130,921 (9.1%)	8,364 (19%)
4	158,448 (0.7%)	980 (2.3%)
Unknown	253,901	1,311
homem_familia	50 (33, 50)	50 (50, 67)	<0.001
raca			<0.001
1	7,511,218 (35%)	16,360 (40%)
2	1,639,328 (7.5%)	2,906 (7.1%)
3	84,243 (0.4%)	119 (0.3%)
4	12,410,311 (57%)	21,300 (52%)
5	103,357 (0.5%)	532 (1.3%)
Unknown	1,798,823	3,635
regiao			<0.001
1	2,168,725 (9.2%)	3,422 (7.6%)
2	8,910,686 (38%)	16,280 (36%)
3	7,927,500 (34%)	12,806 (29%)
4	1,630,446 (6.9%)	3,209 (7.2%)
5	2,887,724 (12%)	9,082 (20%)
Unknown	22,199	53
urban			<0.001
0	6,389,463 (27%)	14,334 (32%)
1	17,157,817 (73%)	30,518 (68%)
idade			<0.001
1	6,363,885 (27%)	14,760 (33%)
2	10,371,942 (45%)	19,112 (43%)
3	4,515,290 (19%)	8,113 (18%)
4	1,933,978 (8.3%)	2,214 (5.0%)
Unknown	362,185	653
escolaridade			<0.001
0	2,894,171 (14%)	5,261 (14%)
1	34,516 (0.2%)	60 (0.2%)
2	137,415 (0.7%)	240 (0.6%)
3	7,189,028 (35%)	15,646 (41%)
4	6,659,739 (32%)	13,689 (36%)
5	3,880,640 (19%)	3,564 (9.3%)
Unknown	2,751,771	6,392
living_conditions			<0.001
0	2,050,291 (8.7%)	3,960 (8.8%)
1	1,880,562 (8.0%)	4,339 (9.7%)
2	4,743,067 (20%)	9,979 (22%)
3	14,873,360 (63%)	26,574 (59%)
1 n (%); Median (IQR)
2 Pearson's Chi-squared test; Wilcoxon rank sum test

NA
NA

3. Regressao Logistica

3.1 Construcao do modelo

#modelagem da regressao logistica
m1 <- glm(desfecho ~ homicide + pessoas_na_familia + homem_familia + raca + regiao + urban + idade + escolaridade + living_conditions, data = dados, family = binomial(link = "logit"))

summary(m1)

Call:
glm(formula = desfecho ~ homicide + pessoas_na_familia + homem_familia + 
    raca + regiao + urban + idade + escolaridade + living_conditions, 
    family = binomial(link = "logit"), data = dados)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.2948  -0.0677  -0.0522  -0.0442   4.1510  

Coefficients:
                      Estimate Std. Error  z value Pr(>|z|)    
(Intercept)         -7.4199077  0.0401221 -184.933  < 2e-16 ***
homicide1            0.3800119  0.0363318   10.459  < 2e-16 ***
pessoas_na_familia2  0.6418170  0.0130467   49.194  < 2e-16 ***
pessoas_na_familia3  1.1496823  0.0168863   68.084  < 2e-16 ***
pessoas_na_familia4  1.5703016  0.0387338   40.541  < 2e-16 ***
homem_familia        0.0118158  0.0002788   42.383  < 2e-16 ***
raca2               -0.1759061  0.0233661   -7.528 5.14e-14 ***
raca3               -0.2088755  0.1037400   -2.013  0.04407 *  
raca4               -0.1893797  0.0133215  -14.216  < 2e-16 ***
raca5                0.8278152  0.0519316   15.940  < 2e-16 ***
regiao2              0.1176746  0.0224504    5.242 1.59e-07 ***
regiao3              0.1549616  0.0234468    6.609 3.87e-11 ***
regiao4              0.3362446  0.0288597   11.651  < 2e-16 ***
regiao5              0.7091870  0.0254315   27.886  < 2e-16 ***
urban1              -0.0759177  0.0158873   -4.779 1.77e-06 ***
idade2               0.0953356  0.0133521    7.140 9.33e-13 ***
idade3              -0.0711234  0.0168181   -4.229 2.35e-05 ***
idade4              -0.4550333  0.0299688  -15.184  < 2e-16 ***
escolaridade1        0.1948047  0.1316932    1.479  0.13908    
escolaridade2        0.2092865  0.0688573    3.039  0.00237 ** 
escolaridade3        0.1549298  0.0204899    7.561 3.99e-14 ***
escolaridade4        0.0425457  0.0207425    2.051  0.04025 *  
escolaridade5       -0.5558247  0.0256448  -21.674  < 2e-16 ***
living_conditions1   0.0192234  0.0298007    0.645  0.51888    
living_conditions2  -0.0215566  0.0267848   -0.805  0.42093    
living_conditions3  -0.0367793  0.0275581   -1.335  0.18200    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 490582  on 18604863  degrees of freedom
Residual deviance: 475979  on 18604838  degrees of freedom
  (4987268 observations deleted due to missingness)
AIC: 476031

Number of Fisher Scoring iterations: 9

3.2 Tabela de significancia

#construindo tabela
tab_model(m1, file="tabela-suicidio.html")

Profiled confidence intervals may take longer time to compute. Use 'ci_method="wald"'
  for faster computation of CIs.

htmltools::includeHTML("tabela-suicidio.html")

	desfecho
Predictors	Odds Ratios	CI	p
(Intercept)	0.00	0.00 – 0.00	<0.001
homicide [1]	1.46	1.36 – 1.57	<0.001
pessoas na familia [2]	1.90	1.85 – 1.95	<0.001
pessoas na familia [3]	3.16	3.05 – 3.26	<0.001
pessoas na familia [4]	4.81	4.45 – 5.18	<0.001
homem familia	1.01	1.01 – 1.01	<0.001
raca [2]	0.84	0.80 – 0.88	<0.001
raca [3]	0.81	0.66 – 0.99	0.044
raca [4]	0.83	0.81 – 0.85	<0.001
raca [5]	2.29	2.06 – 2.53	<0.001
regiao [2]	1.12	1.08 – 1.18	<0.001
regiao [3]	1.17	1.12 – 1.22	<0.001
regiao [4]	1.40	1.32 – 1.48	<0.001
regiao [5]	2.03	1.93 – 2.14	<0.001
urban [1]	0.93	0.90 – 0.96	<0.001
idade [2]	1.10	1.07 – 1.13	<0.001
idade [3]	0.93	0.90 – 0.96	<0.001
idade [4]	0.63	0.60 – 0.67	<0.001
escolaridade [1]	1.22	0.93 – 1.56	0.139
escolaridade [2]	1.23	1.07 – 1.41	0.002
escolaridade [3]	1.17	1.12 – 1.22	<0.001
escolaridade [4]	1.04	1.00 – 1.09	0.040
escolaridade [5]	0.57	0.55 – 0.60	<0.001
living conditions [1]	1.02	0.96 – 1.08	0.519
living conditions [2]	0.98	0.93 – 1.03	0.421
living conditions [3]	0.96	0.91 – 1.02	0.182
Observations	18604864
R2 Tjur	0.001

NA

3.3 Plotando o grafico com Odds Ratio

theme_set(theme_sjplot())

plot_model(m1, colors = "bw", show.values = TRUE, value.offset = .5, title = "Suicide in the family", axis.labels = c("Score 3 - Adequate", "Score 2 - Bad", "Score 1 - Very bad", "≥ High school", "4-9 y of formal edu.", "<4 y of formal edu.", "Primary school", "Pre-school", "≥ 60", "40-59", "20-39", "Urban", "South", "Central-West", "Southeast", "North", "Indigenous", "Mixed/brown", "Asian descendent", "Black", "Male", "≥ 7 members", "5-6 members", "3-4 members", "Homicide"))

3.4 Verificando a qualidade do ajuste do modelo

3.4.1 Ausencia de outliers e verificando os residuos padronizados

#verificando a presenca de outliers
#plot(m1, which = 5)

#verificando os residuos padronizados
summary(stdres(m1))

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.21140 -0.04806 -0.03704 -0.00006 -0.03139 74.52982 

3.4.2 Ausencia de multicolinearidade

Para realizar o pressuposto de multiconearidade utilizamos a funcao “vif()” do pacote car (companion to applied regression). Quando VIF > 0 podemos dizer que existiu multicolinearidade entre variaveis.

# Mulricolinearidade

m1 %>% 
vif() %>% 
    kable(
    format = "markdown",
    caption = "Diagnosis of multicollinearity"
  )

Diagnosis of multicollinearity

	GVIF	Df	GVIF^(1/(2*Df))
homicide	1.015439	1	1.007690
pessoas_na_familia	1.196324	3	1.030326
homem_familia	1.024913	1	1.012380
raca	1.405512	4	1.043468
regiao	1.542113	4	1.055637
urban	1.714921	1	1.309550
idade	1.245282	3	1.037237
escolaridade	1.157303	5	1.014717
living_conditions	1.848228	3	1.107795

#vif(m1)

3.4.3 Aplicacao da ANOVA

# Contribution of each independent variable

m1 %>% 
anova(test = "Chisq")%>%
  kable(
    format = "markdown",
    caption = "Anova: Number of visits by diagnosis as outcome"
  )

Anova: Number of visits by diagnosis as outcome

	Df	Deviance	Resid. Df	Resid. Dev	Pr(>Chi)
NULL	NA	NA	18604863	490582.5	NA
homicide	1	354.644738	18604862	490227.8	0.0000000
pessoas_na_familia	3	7366.666842	18604859	482861.2	0.0000000
homem_familia	1	2151.681017	18604858	480709.5	0.0000000
raca	4	1195.314970	18604854	479514.2	0.0000000
regiao	4	1507.582868	18604850	478006.6	0.0000000
urban	1	183.599711	18604849	477823.0	0.0000000
idade	3	329.211529	18604846	477493.8	0.0000000
escolaridade	5	1508.789341	18604841	475985.0	0.0000000
living_conditions	3	6.287911	18604838	475978.7	0.0984127

3.4.3 Teste de hosmer lemeshow

library(ResourceSelection)

ResourceSelection 0.3-5      2019-07-22

h1 <- hoslem.test(m1$y, fitted(m1))
h1

    Hosmer and Lemeshow goodness of fit (GOF) test

data:  m1$y, fitted(m1)
X-squared = 26.39, df = 8, p-value = 0.0009004

library(generalhoslem) 

Loading required package: reshape

Attaching package: ‘reshape’

The following object is masked from ‘package:lubridate’:

    stamp

The following object is masked from ‘package:data.table’:

    melt

The following object is masked from ‘package:dplyr’:

    rename

The following objects are masked from ‘package:tidyr’:

    expand, smiths

h2 <- logitgof(m1$y, fitted(m1))
h2

    Hosmer and Lemeshow test (binary model)

data:  m1$y, fitted(m1)
X-squared = 26.39, df = 8, p-value = 0.0009004

library(performance)

performance::performance_hosmer(m1)

# Hosmer-Lemeshow Goodness-of-Fit Test

  Chi-squared: 26.390
           df:  8    
      p-value:  0.001

Summary: model does not fit well.

model_performance(m1)

# Indices of model performance

AIC       |       BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
--------------------------------------------------------------------------------------------------
4.760e+05 | 4.764e+05 | 9.600e-04 | 0.042 | 0.160 |    0.013 |      -Inf |       8.502e-06 | 0.996