Introdução

Os dados utilizados contém informações sobre variáveis potencialmente associadas à presença de doença cardiovascular (DCV) em homens adultos residentes em uma determinada comunidade. O quadro abaixo descreve as variáveis disponíveis, contendo informações aferidas na linha de base do estudo, isto é, sem considerar o seguimento dos indivíduos ao longo do tempo.

VARIÁVEL	DESCRIÇÃO
age	Idade (anos)
totchol	Colesterol total (mmol/l)
bmi	Índice de massa corporal (\(kg/m^2\))
systol	Pressão arterial sistólica (mmHg)
smoking	Categorias do tabagismo: 1 = Nunca fumou 2 = Ex-tabagista 3 = Tabagista
activity	Autorrelato do nível de atividade física: 1 = Ativo 2 = Atividade média 3 = Inativo
chd	Incidência de doença cardiovascular (DCV) 0 = Não 1 = Sim

Visualização do banco de dados
id	age	totchol	bmi	systol	smoking	activity
1	46	8.660333	30.83653	128	1	1
2	42	6.243496	21.67126	129	3	2
3	45	7.748424	26.54321	121	3	2
4	41	7.300000	19.03811	123	3	2
5	52	6.629518	24.96801	132	1	2
6	48	4.833674	22.64738	147	3	1

Análise Descritiva

A seguir, iremos realizar uma análise descritiva para uma melhor compreensão dos dados que estamos lidando. Iremos primeiro analisar as variáveis numéricas utilizando gráficos de box-plot, e as variáveis categóricas utilizando graficos de barra. Equivale à Tabela 1 pedida no item c) do exercício proposto.

Variáveis numéricas

As variáveis numericas presentes são: ‘age’, ‘totchol’, ‘bmi’ e ‘systol’.

##### SEPARANDO AS VARIÁVEIS NUMÉRICAS
var_num <- dados_tabagismo %>% dplyr::select("age", "totchol", "bmi", "systol")


#### DESCRITIVA
aux <- stat.desc(var_num)

aux <- aux[c(4, 5, 8, 9, 12:14), ]
rownames(aux) <- c("Mínimo", "Máximo", "Mediana", "Média", "Variância", "Desvio Padrão", "Coeficiente de Variação")
kable(aux, digits = 3, align = "c", caption = "Principais Estatísticas descritivas do banco de dados")

Principais Estatísticas descritivas do banco de dados
	age	totchol	bmi	systol
Mínimo	40.000	2.848	15.549	82.000
Máximo	59.000	12.445	42.593	221.000
Mediana	49.000	6.294	25.826	131.000
Média	49.111	6.372	26.025	132.924
Variância	32.232	1.335	11.113	323.151
Desvio Padrão	5.677	1.155	3.334	17.976
Coeficiente de Variação	0.116	0.181	0.128	0.135

Abaixo, conseguimos visualizar o gráfico box-plot de cada variável em sua própria escala.

# creating a plot
aux2 <- gather(var_num)
colnames(aux2) <- c("Variável", "Valor")


ggplot(aux2) +
  geom_boxplot(aes(x = Variável, y = Valor, colour = Variável)) +
  facet_wrap(~Variável, scale = "free", switch = NULL) +
  ggtitle("Variáveis contínuas") +
  theme_ipsum() +
  theme(
    strip.text.x = element_blank()
  )

Variáveis categóricas

As variáveis categóricas presentes são: ‘smoking’, ‘activity’ e ‘age’.

Abaixo, temos uma gráfico de barras em que conseguimos visualizar a frequência da cada resposta.

var_cat <- c("smoking", "activity", "chd")

aux1 <- dados_tabagismo %>% dplyr::select(c(var_cat))

aux1 <- lapply(aux1, function(x) as.factor(x))

aux1 <- as.data.frame(aux1)


aux1 %>%
  gather(key = "Variável", value = "Valor", smoking, activity, chd) %>%
  ggplot() +
  geom_bar(mapping = aes(x = Valor, fill = Variável)) +
  facet_wrap(~Variável, scale = "free_x") +
  scale_fill_viridis(discrete = T) +
  ggtitle("Variáveis categóricas") +
  theme_ipsum() +
  xlab("")

VARIÁVEL	DESCRIÇÃO
smoking	Categorias do tabagismo: 1 = Nunca fumou 2 = Ex-tabagista 3 = Tabagista
activity	Autorrelato do nível de atividade física: 1 = Ativo 2 = Atividade média 3 = Inativo
chd	Incidência de doença cardiovascular (DCV) 0 = Não 1 = Sim

smoking

Categorias do tabagismo:

1 = Nunca fumou

2 = Ex-tabagista

3 = Tabagista

activity

Autorrelato do nível de atividade física:

1 = Ativo

2 = Atividade média

3 = Inativo

chd

Incidência de doença cardiovascular (DCV)

0 = Não

1 = Sim

Correlação entre as variáveis

Podemos começar nossas análises explorando a associação entre as variáveis através da representação gráfica da matriz de correlações de Pearson

M <- dados_tabagismo %>%
     dplyr::select(-c("id")) %>%
     cor()

corrplot(M, title = "Correlação entre as Variáveis", mar = c(0, 0, 2, 0), type = "full", order = "hclust", col = brewer.pal(n = 8, name = "GnBu"))

Podemos notar no gráfico (acima) que as variáveis preditoras não estão significativamente correlacionadas, o que nos dá indícios de que multicolinearidade não será um problema na parametrização dos nossos modelos.

Regressão Logística

A regressão logística é uma técnica estatística amplamente utilizada na análise de dados quando a variável dependente é categórica, binária ou multinomial. Ela é útil em situações em que estamos interessados em entender como uma variável categórica está relacionada a uma ou mais variáveis explicativas.

Como o objetivo é estimar associação entre a concentração de colesterol total e a presença de DCV, sendo esta útilma um variável binária e a ocorrência do desfecho de interesse sendo pouco comum (menor que 10%), podemos prosseguir com esse modelo, no qual utilizaremos a função “glm” (Modelo Linear Generalizado), junto da tabela de razão de chances em que estimaremos a associação entre as variáveis. Equivale em cada modelo à Tabela 2 pedida no item c) e ao item b) do exercício proposto.

Modelo 1

modelo1 <- glm(factor(chd) ~ totchol, family = binomial(logit), data = dados_tabagismo)
summary(modelo1)

## 
## Call:
## glm(formula = factor(chd) ~ totchol, family = binomial(logit), 
##     data = dados_tabagismo)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.29894    0.47834 -11.078  < 2e-16 ***
## totchol      0.35101    0.06899   5.088 3.61e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1138.7  on 2950  degrees of freedom
## Residual deviance: 1114.0  on 2949  degrees of freedom
## AIC: 1118
## 
## Number of Fisher Scoring iterations: 6

Razão de Chances

O modelo traz os resultados dos estimadores na forma logarítma, ou seja, o log das chances da variável. No entanto, para uma melhor interpretação da relação do colesterol com o CHD é necessária a transformação deste coeficiente efetuando a exponenciação da variavel da regressão. Assim, obtém-se a razão das chances (OR - Odds Ratio em inglês) para as variáveis independentes.

aux3 <- cbind(RC = exp(modelo1$coefficients), exp(confint(modelo1)))
aux3 <- cbind(aux3, RC_AJUSTADO = lapply(data.frame(aux3)$RC, function(x) (x - 1) * 100))

kable(aux3, align = "c")

	RC	2.5 %	97.5 %	RC_AJUSTADO
(Intercept)	0.00499688918404385	0.00194355819868628	0.0126957747710203	-99.5003110815956
totchol	1.42050801727679	1.23955866971147	1.62496243301668	42.0508017276789

O resultado acima evidencia que para uma alteração em 1 (uma) unidade em \(totchol\) (colesterol), a chance de que y seja igual a 1 aumenta em 42% ((1,42-1)*100), assim como podemos ver na coluna “RC_AJUSTADO”.

Modelo 2

modelo2 <- glm(factor(chd) ~ totchol + age, family = binomial(logit), data = dados_tabagismo)
summary(modelo2)

## 
## Call:
## glm(formula = factor(chd) ~ totchol + age, family = binomial(logit), 
##     data = dados_tabagismo)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -6.53052    0.91447  -7.141 9.24e-13 ***
## totchol      0.35444    0.06940   5.107 3.27e-07 ***
## age          0.02444    0.01529   1.598     0.11    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1138.7  on 2950  degrees of freedom
## Residual deviance: 1111.4  on 2948  degrees of freedom
## AIC: 1117.4
## 
## Number of Fisher Scoring iterations: 6

Razão de Chances

aux4 <- cbind(RC = exp(modelo2$coefficients), exp(confint(modelo2)))
aux4 <- cbind(aux4, RC_AJUSTADO = lapply(data.frame(aux4)$RC, function(x) (x - 1) * 100))
kable(aux4, align = "c")

	RC	2.5 %	97.5 %	RC_AJUSTADO
(Intercept)	0.00145824375136498	0.000237474692738171	0.00858890731745356	-99.8541756248635
totchol	1.42538168745875	1.24282885064722	1.63186469645102	42.5381687458753
age	1.02473789929805	0.994542044802386	1.05606068780656	2.4737899298052

O resultado acima evidencia que para uma alteração em 1 (uma) unidade em \(totchol\) (colesterol), a chance de que y seja igual a 1 aumenta em 42% ((1,42-1)**100). Já para a variável \(age\) (idade), a cada unidade aumentado, temos um aumento de 2,4% a chance de que de ocorrência de doença caridovascular (y = 1). As demais probabilidades podemos ver na coluna “RC_AJUSTADO”

Modelo 3

modelo3 <- glm(factor(chd) ~ totchol + age + systol + factor(smoking) + factor(activity), family = binomial(logit), data = dados_tabagismo)
summary(modelo3)

## 
## Call:
## glm(formula = factor(chd) ~ totchol + age + systol + factor(smoking) + 
##     factor(activity), family = binomial(logit), data = dados_tabagismo)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -9.234853   1.020510  -9.049  < 2e-16 ***
## totchol            0.311544   0.070182   4.439 9.03e-06 ***
## age                0.005742   0.015968   0.360  0.71913    
## systol             0.024619   0.004425   5.563 2.65e-08 ***
## factor(smoking)2   0.432998   0.296431   1.461  0.14410    
## factor(smoking)3   0.792264   0.260006   3.047  0.00231 ** 
## factor(activity)2 -0.069441   0.214563  -0.324  0.74621    
## factor(activity)3  0.103823   0.272123   0.382  0.70281    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1138.7  on 2950  degrees of freedom
## Residual deviance: 1070.3  on 2943  degrees of freedom
## AIC: 1086.3
## 
## Number of Fisher Scoring iterations: 6

- Razão de Chances

aux5 <- cbind(RC = exp(modelo3$coefficients), exp(confint(modelo3)))
aux5 <- cbind(aux5, RC_AJUSTADO = lapply(data.frame(aux5)$RC, function(x) (x - 1) * 100))
kable(aux5, align = "c")

	RC	2.5 %	97.5 %	RC_AJUSTADO
(Intercept)	9.75785331127465e-05	1.28727810588835e-05	0.000706084736420388	-99.9902421466887
totchol	1.36553248866415	1.18908927275423	1.56614025693593	36.5532488664147
age	1.0057589959597	0.974787559342455	1.03783690054484	0.575899595970042
systol	1.02492406414737	1.01599224044113	1.03379526792123	2.49240641473718
factor(smoking)2	1.54187255805916	0.871110024710928	2.80443981370446	54.1872558059163
factor(smoking)3	2.20839149396134	1.35682382707733	3.78093773743459	120.839149396134
factor(activity)2	0.932915373697504	0.618258959384974	1.43762386079243	-6.70846263024963
factor(activity)3	1.10940367560657	0.646414422702671	1.88750939736712	10.9403675606567

Podemos ver pela tabela acima o aumento (ou redução) da probabilidade de ocorrência de y = 1 na coluna “RC_AJUSTADO”, a partir do aumento aumento em uma unidade de cada variável em questão.

Modelo 3 com Interação

modelo3_int <- glm(factor(chd) ~ totchol + age + systol + factor(smoking) + factor(activity) + totchol * systol, family = binomial(logit), data = dados_tabagismo)
summary(modelo3_int)

## 
## Call:
## glm(formula = factor(chd) ~ totchol + age + systol + factor(smoking) + 
##     factor(activity) + totchol * systol, family = binomial(logit), 
##     data = dados_tabagismo)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -12.557873   3.344759  -3.754 0.000174 ***
## totchol             0.809033   0.482116   1.678 0.093330 .  
## age                 0.005677   0.015966   0.356 0.722148    
## systol              0.048113   0.022896   2.101 0.035608 *  
## factor(smoking)2    0.435434   0.296146   1.470 0.141472    
## factor(smoking)3    0.788548   0.259778   3.035 0.002402 ** 
## factor(activity)2  -0.061472   0.214654  -0.286 0.774590    
## factor(activity)3   0.116168   0.272344   0.427 0.669708    
## totchol:systol     -0.003517   0.003375  -1.042 0.297352    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1138.7  on 2950  degrees of freedom
## Residual deviance: 1069.2  on 2942  degrees of freedom
## AIC: 1087.2
## 
## Number of Fisher Scoring iterations: 6

Razão de Chances

aux6 <- cbind(RC = exp(modelo3_int$coefficients), exp(confint(modelo3_int)))
aux6 <- cbind(aux6, RC_AJUSTADO = lapply(data.frame(aux6)$RC, function(x) (x - 1) * 100))
kable(aux6, align = "c")

	RC	2.5 %	97.5 %	RC_AJUSTADO
(Intercept)	3.5171030401279e-06	5.16899834472897e-09	0.00256815781316754	-99.999648289696
totchol	2.245735569747	0.867549853120531	5.74496417850452	124.5735569747
age	1.00569331955937	0.974729768044457	1.03776431470892	0.569331955936803
systol	1.04928926362151	1.00272550988944	1.09699034580784	4.92892636215083
factor(smoking)2	1.5456330624385	0.87371377535939	2.80973724201806	54.56330624385
factor(smoking)3	2.20019985922277	1.35238210335861	3.76526468392048	120.019985922277
factor(activity)2	0.940379484096669	0.623107124191283	1.44940864533774	-5.96205159033314
factor(activity)3	1.12318463702314	0.654197214736472	1.9118814413823	12.3184637023139
totchol:systol	0.996488944900573	0.989941552144605	1.00314012522968	-0.351105509942717

A inclusão da interação entre a variável totchol (colesterol) e systol não foi significativa para o modelo. Podemos ver pela tabela acima o aumento (ou redução) da probabilidade de ocorrência de y = 1 na coluna “RC_AJUSTADO”, a partir do aumento aumento em uma unidade de cada variável em questão.

Modelo 3 com a Variável Colesterol Categorizada

Conforme pedido no exercício, a variável colesterol foi categorizada de acordo com os quantis teóricos.

quantis <- quantile(dados_tabagismo$totchol)

dados_tabagismo <- mutate(dados_tabagismo, totchol_CAT = ifelse(totchol <= quantis[2], 1,
                                                          ifelse(totchol <= quantis[3], 2,
                                                           ifelse(totchol <= quantis[4], 3, 4)
                                                           )
                                                           ))

modelo3_cat <- glm(factor(chd) ~ factor(totchol_CAT) + age + systol + factor(smoking) + factor(activity), family = binomial(logit), data = dados_tabagismo)
summary(modelo3_cat)

## 
## Call:
## glm(formula = factor(chd) ~ factor(totchol_CAT) + age + systol + 
##     factor(smoking) + factor(activity), family = binomial(logit), 
##     data = dados_tabagismo)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -7.864445   0.933777  -8.422  < 2e-16 ***
## factor(totchol_CAT)2  0.625548   0.308529   2.028  0.04261 *  
## factor(totchol_CAT)3  0.749093   0.300235   2.495  0.01259 *  
## factor(totchol_CAT)4  1.161250   0.283836   4.091 4.29e-05 ***
## age                   0.004369   0.015969   0.274  0.78440    
## systol                0.024895   0.004420   5.632 1.78e-08 ***
## factor(smoking)2      0.413332   0.296492   1.394  0.16329    
## factor(smoking)3      0.779760   0.259988   2.999  0.00271 ** 
## factor(activity)2    -0.060556   0.214621  -0.282  0.77783    
## factor(activity)3     0.135497   0.271477   0.499  0.61770    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1138.7  on 2950  degrees of freedom
## Residual deviance: 1069.5  on 2941  degrees of freedom
## AIC: 1089.5
## 
## Number of Fisher Scoring iterations: 6

Razão de Chances

aux7 <- cbind(RC = exp(modelo3_cat$coefficients), exp(confint(modelo3_cat)))
aux7 <- cbind(aux7, RC_AJUSTADO = lapply(data.frame(aux7)$RC, function(x) (x - 1) * 100))
kable(aux7, align = "c", caption = "Razão de Chance")

Razão de Chance
	RC	2.5 %	97.5 %	RC_AJUSTADO
(Intercept)	0.000384162398223441	6.02658975015615e-05	0.0023531828234168	-99.9615837601776
factor(totchol_CAT)2	1.8692692626856	1.03366007662204	3.49379925598762	86.9269262685603
factor(totchol_CAT)3	2.11508122088432	1.19268463773291	3.90013888222246	111.508122088432
factor(totchol_CAT)4	3.19392244373777	1.87123244816547	5.73316971265598	219.392244373777
age	1.00437833884887	0.973441878178512	1.03640655817231	0.437833884886785
systol	1.02520713570187	1.01628613821358	1.03407246892538	2.52071357018739
factor(smoking)2	1.51184760048303	0.854025818031547	2.75010827556651	51.1847600483026
factor(smoking)3	2.18094927737806	1.33996260843053	3.73374018384052	118.094927737806
factor(activity)2	0.941240900468348	0.623724335387774	1.45065964554431	-5.87590995316524
factor(activity)3	1.1451054490519	0.668174387925964	1.94615385471416	14.5105449051897

Podemos ver pela tabela acima o aumento (ou redução) da probabilidade de ocorrência de y = 1 na coluna “RC_AJUSTADO”, a partir do aumento aumento em uma unidade de cada variável em questão.

# - **Gráfico de associação a partir da Categorização**
# 
# ggplot(dados_tabagismo, aes(x = totchol, y = chd, colour = factor(totchol_CAT))) +
#   geom_point() +
#   geom_smooth(method = lm, color = "red", fill = "#69b3a2", se = TRUE) +
#   theme_ipsum()
# 
# A alternativa de trabalharmos com a variável colesterol categorizada iria produzir resultados viesados, pois estaríamos considerando uma mesma probabilidade de ocorrência para os elementos presentes em cada um dos grupos. Podemos ver, a partir da linha vermelha presente no gráfico, o comportamento contínuo da variável em questão.

Pressupostos

Abaixo, iremos verificar alguns dos pressupostos solicitados para a validação do modelo.

Multicolinearidade

vif(modelo2)

##  totchol      age 
## 1.001216 1.001216

vif(modelo3)

##                      GVIF Df GVIF^(1/(2*Df))
## totchol          1.010141  1        1.005058
## age              1.089439  1        1.043762
## systol           1.084502  1        1.041394
## factor(smoking)  1.017926  2        1.004452
## factor(activity) 1.014303  2        1.003557

vif(modelo3_int)

##                       GVIF Df GVIF^(1/(2*Df))
## totchol          47.953409  1        6.924840
## age               1.092760  1        1.045352
## systol           29.256065  1        5.408888
## factor(smoking)   1.018635  2        1.004626
## factor(activity)  1.015745  2        1.003913
## totchol:systol   75.346017  1        8.680208

vif(modelo3_cat)

##                         GVIF Df GVIF^(1/(2*Df))
## factor(totchol_CAT) 1.009372  3        1.001556
## age                 1.088325  1        1.043228
## systol              1.085220  1        1.041739
## factor(smoking)     1.021625  2        1.005363
## factor(activity)    1.011910  2        1.002964

Como podemos ver acima, todos os índices foram menores do que 10, o que indica ausência de multicolinearidade, assim como havíamos concluído anteriormente a partir do gráfico de correlações.

Ajuste Global

hnp(modelo1)

## Binomial model

hnp(modelo2)

## Binomial model

hnp(modelo3)

## Binomial model

hnp(modelo3_int)

## Binomial model

hnp(modelo3_cat)

## Binomial model

Pelos gráficos acima, podemos ver que todos os resíduos dos modelos se adequaram ao ajuste proposto (“binomial(logit)”), estando todos os pontos dentro do envelope de confiança.

Conclusão

Para definirmos qual foi o melhor modelo, iremos compará-los a partir do Critério de Akaike (AIC).

aux8 <- cbind(modelo1$aic, modelo2$aic, modelo3$aic, modelo3_int$aic, modelo3_cat$aic)
colnames(aux8) <- c("Modelo 1", "Modelo 2", "Modelo 3", "Modelo 3 Interação", "Modelo 3 Categorização")
rownames(aux8) <- "AIC"

kable(aux8, align = "c")

	Modelo 1	Modelo 2	Modelo 3	Modelo 3 Interação	Modelo 3 Categorização
AIC	1117.974	1117.41	1086.254	1087.174	1089.533

O AIC é calculado como a soma do logaritmo da verossimilhança do modelo e um termo de penalidade proporcional ao número de parâmetros do modelo. O modelo com o menor valor de AIC é geralmente considerado o melhor modelo candidato. Portanto, temos que o melhor candidato é o Modelo 3.

Referente às limitações dos modelos, como as variáveis disponíveis contém informações aferidas na linha de base do estudo, isto é, sem considerar o seguimento dos indivíduos ao longo do tempo, nos foi possível caluclar a associação entre as variáveis pela Razão de Chances, enquanto que para obtermos estimativas mais exatas, seria recomendado o uso da Razão de Prevalência.

Bioestatística - Avaliação 1

João Lucas Simon

18/02/2023