1 Conceitos básicos

Este trabalho utiliza técnicas de Análise de Sobrevivência para dados dos alunos do bacharelado em Ciência e Tecnologia da UFERSA. São apresentados, nesta seção, alguns conceitos importantes que permitirão a compreensão do estudo.

  • Tempo de falha: variável resposta que representa o tempo até a ocorrência de um evento de interesse. Existem três elementos constituintes do tempo de falha: (i) o tempo inicial; (ii) a escala de medida e (iii) o evento de interesse.
  • Censura: observação parcial de uma resposta. Acontece em situações onde, por alguma razão, o acompanhamento do paciente (no caso de estudos clínicos) foi interrompido. O paciente pode morrer de uma causa diferente da estudada, pode mudar de cidade, etc. Existem três tipos de censura. No Tipo I, o estudo tem um prazo fixo e o tempo de falha de alguns indivíduos é superior a este prazo. No Tipo II, o estudo é encerrado após o evento de interesse acontecer em um dado número de indivíduos. Por fim, o terceiro mecanismo de censura é o tipo aleatório, em que o indivíduo sai do estudo antes da ocorrência do evento de interesse. Em aplicações médicas, isto também pode ocorrer se o paciente morre por uma razão diferente da estudada.

A Análise de Sobrevivência é, portanto, um conjunto de técnicas estatísticas desenvolvidas para estudar dados censurados. Seu principal campo de aplicação é em situações médicas envolvendo dados censurados. No entanto, ocorrem situações similares em que a aplicação destas técnicas é adequada, especialmente em engenharia, ciências sociais, demografia.

Nosso objetivo é utilizar a análise de sobrevivência em dados dos alunos do Bacharelado em Ciência e Tecnologia da UFERSA. O evento de interesse é a conclusão do curso, considerando a data de integralização. A retenção, concebida aqui como o prolongamento do tempo de formatura dos alunos para além do planejado, é o objeto de estudo deste trabalho. Considera-se também a evasão, isto é, o abandono do curso por parte do aluno, como um fator retratado neste estudo.

Para que o tempo inicial fosse o mesmo para todos os indivíduos, foi selecionado o semestre 2015.1 para o estudo. Sobre a escala de tempo, são contados os dias corridos desde o ingresso no curso até o evento de interesse ou a censura. Os alunos que evadiram são considerados como dados censurados, isto é, indivíduos que saíram do estudo antes que o evento de interesse (integralização) ocorresse. Note que fica caraterizada a censura aleatória.

1.1 Função de sobrevivência

A função de sobrevivência é definida como a probabilidade de um indivíduo não falhar até o tempo \(t\), ou seja, que o seu tempo de vida seja superior a \(t\):

\[S(t) = P(T\geq t)\]

1.2 Função de Taxa de Falha ou de Risco

Seja um intervalo de tempo dado por \([t_1, t_2)\). A probabilidade de que a falha ocorra neste intervalo de tempo é dada por \(S(t_1) - S(t_2)\). Assim, definimos a taxa de falha para o intervalo como a probabilidade de que a falha ocorra neste intervalo, dado que não ocorreu antes de \(t_1\), dividida pelo comprimento do intervalo:

\[\frac{S(t_1)-S(t_2)}{(t_2-t_1)S(t_1)}\] Definindo o intervalo como \([t, t+\Delta t)\), temos a função taxa de falha:

\[\lambda(t) = \frac{S(t) - S(t+\Delta t)}{\Delta t S(t)}\]

Para obter uma taxa de falha instantânea, assumimos um valor arbitrariamente pequeno para a variação de tempo \(\Delta t\), obtendo então

\[\lambda(t) = lim_{\Delta t \rightarrow 0}\frac{P(t \leq T \leq t + \Delta t)}{\Delta t}\]

2 Tratamento dos dados

Os dados recebidos correspondem a 16166 registros anônimos de estudantes dos cursos de Bacharelado em Ciência e Tecnologia da Ufersa em 20 variáveis: “sexo”, “raça”, “Necessidade Especial”, “Município Nascimento”, “uf”, “curso”, “campus”, “turno”, “Status discente”, “ano_ingresso”, “periodo_ingresso”, “Data Ingresso”, “Data integralização”, “Data Colacao de grau”, “Forma Ingresso”, “Data movimentação” , “Tipo Movimentação”, “ira”, “classe” e “Reprovações”.

As seguintes manipulações foram feitas nos dados:

  • As variáveis Raça, Necessidades especiais, classe foram binarizadas, isto é, foram agrupadas em duas classes para facilitar a interpretação dos resultados,
  • Foi selecionada a turma de 2015.1 para o estudo. O tempo de falha foi computado da seguinte maneira: para os alunos com status “CONCLUÍDO” foi contabilizado o número de dias desde a data de ingresso até a integralização do curso; para os demais alunos, exceto os ATIVOS, foi computado o número de dias desde o ingresso até a data da última movimentação. Os alunos ATIVOS foram excluídos do estudo.
## Análise de sobrevivencia sobre o curso 
## de Ciência e Tecnologia da Ufersa

# Pacotes -----------------------------------------------------------------

library(readxl)
library(summarytools)
library(dplyr)
library(ggplot2)
library(survival)
library(flexsurv)
library(survminer)


# Abrindo os dados --------------------------------------------------------

dados <- read_excel("./dados.xlsx", 
                    col_types = c("text", "text", "text", 
                                  "text", "text", "text", "text", "text", 
                                  "text", "numeric", "numeric", "date", 
                                  "date", "date", "text", "text", "text", 
                                  "text", "text", "numeric"))

dados$`Data Ingresso` = as.Date(dados$`Data Ingresso`)
dados$`Data integralização` = as.Date(dados$`Data integralização`)
dados$`Data movimentação` = as.Date(dados$`Data movimentação`)


# print(dfSummary(dados, valid.col = FALSE, varnumbers = FALSE, graph.col = FALSE), col.widths=c(200,200,200,200), 
#       max.tbl.height = 300, max.tbl.width = 500, method = "render")

# Abrindo os dados  -----------------------------------------------

dados %>% filter(ano_ingresso == 2015 & periodo_ingresso == 1 & sexo %in% c("M","F"), 
                 turno=="Matutino e Vespertino", 
                 `Status discente` %in% c("CONCLUÍDO", "CANCELADO", "TRANCADO")) %>% 
  mutate(tempo = ifelse(is.na(`Data integralização`),
                        `Data movimentação` - `Data Ingresso`,
                        `Data integralização` - `Data Ingresso`), 
         censura = ifelse(`Status discente` == "CONCLUÍDO", 1, 0)) -> surv_data

#View(surv_data)

colnames(surv_data) = c("sexo", "raça", "pne", "Município Nascimento",
                        "uf", "curso", "campus", "turno", "Status discente", "ano_ingresso", 
                        "periodo_ingresso", "Data Ingresso", "Data integralização", 
                        "Data Colacao de grau", "Forma Ingresso", "Data movimentação",
                        "Tipo Movimentação", "ira", "classe", "Reprovações", "tempo", "censura")

# binarizar os regressores

surv_data$raça <- ifelse(surv_data$raça == "Branco", "Branco", "Nao_branco")
surv_data$pne <- ifelse(surv_data$pne == "NENHUMA", "NAO", "SIM")
surv_data$classe <- ifelse(surv_data$classe %in% c("CLASSE C", "CLASSE D","CLASSE E"),
                           "DE", "ABC")
surv_data$ira <- as.numeric(surv_data$ira)

surv_data %>% filter(!is.na(surv_data$tempo)) -> surv_data

3 Estimador de Kaplan-Meier

Grosso modo, a análise descritiva consiste em apresentar medidas de tendência central e variabilidade dos dados utilizados no estudo. No entanto, a presença de censura nos dados inviabiliza este procedimento. Assim, em dados de tempo de vida, o principal método de análise exploratória é a representação da função de sobrevivência.

Para dados não censurados, a função de sobrevivência empírica é dada pela divisão entre o número de observações que não falharam até o tempo \(t\) e o total de observações no estudo. Para dados censurados, podemos utilizar a técnica não-paramétrica conhecida como Estimador de Kaplan-Meier:

\[\hat{S}(t) = \prod_{j:t_j<t}\left(\frac{n_j-d_j}{n_j}\right) = \prod_{j:t_j<t}\left(1-\frac{d_j}{n_j}\right)\]

em que

  • \(t_1 < \ldots < t_k\) são os tempos distintos e ordenados de falha
  • \(d_j\) é o número de falhas em \(t_j\), \(j=1,\ldots,k\), e
  • \(n_j\) é o número de indivíduos que não falharam e não foram censurados até o instante imediatamente anterior a \(t\).

Para os alunos de Ciência e Tecnologia ingressantes em 2015.1, temos a seguinte função de sobrevivência obtida por meio do estimador de Kaplan-Meier:

## Call: survfit(formula = surv_object ~ 1, data = surv_data)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   465    354       1   0.9972 0.00282       0.9917        1.000
##   630    338       3   0.9883 0.00580       0.9770        1.000
##   644    335       3   0.9795 0.00768       0.9645        0.995
##   815    268       1   0.9758 0.00848       0.9593        0.993
##   966    257       1   0.9720 0.00925       0.9541        0.990
##   970    240       4   0.9558 0.01214       0.9323        0.980
##   995    235       1   0.9518 0.01275       0.9271        0.977
##  1133    230       1   0.9476 0.01335       0.9218        0.974
##  1138    222       1   0.9433 0.01395       0.9164        0.971
##  1140    221       3   0.9305 0.01560       0.9005        0.962
##  1141    218       4   0.9135 0.01750       0.8798        0.948
##  1143    214      21   0.8238 0.02437       0.7774        0.873
##  1145    193       2   0.8153 0.02486       0.7680        0.865
##  1148    190      55   0.5793 0.03212       0.5196        0.646
##  1153    127       1   0.5747 0.03219       0.5150        0.641
##  1289    111       1   0.5695 0.03231       0.5096        0.637
##  1299    108      13   0.5010 0.03355       0.4394        0.571
##  1301     95       1   0.4957 0.03361       0.4340        0.566
##  1302     94      15   0.4166 0.03389       0.3552        0.489
##  1304     79       7   0.3797 0.03364       0.3192        0.452
##  1305     71       1   0.3744 0.03359       0.3140        0.446
##  1477     62       3   0.3562 0.03355       0.2962        0.428
##  1484     55       9   0.2979 0.03321       0.2395        0.371
##  1485     46       2   0.2850 0.03301       0.2271        0.358
##  1486     44       3   0.2656 0.03261       0.2088        0.338
##  1487     40       4   0.2390 0.03194       0.1839        0.311
##  1625     33       8   0.1811 0.03005       0.1308        0.251
##  1631     21       1   0.1724 0.02983       0.1228        0.242
##  1634     20       1   0.1638 0.02956       0.1150        0.233
##  1639     17       1   0.1542 0.02935       0.1062        0.224
##  1795     16       1   0.1445 0.02906       0.0975        0.214
##  1800     15       2   0.1253 0.02820       0.0806        0.195
##  1801     13       1   0.1156 0.02763       0.0724        0.185
##  1803     12       2   0.0964 0.02617       0.0566        0.164
##  1813      8       2   0.0723 0.02455       0.0371        0.141
##  1814      6       1   0.0602 0.02323       0.0283        0.128

3.1 Comparando grupos

Podemos particionar os dados dos estudantes de acordo com os valores apresentados por uma variável de interesse, por exemplo, sexo ou raça. Com isso, é possível obter as funções de sobrevivência para diferentes e fazer comparações entre os tempos de vida dos dois grupos.

Esta comparação entre os grupos pode ser feita por testes de hipóteses, sendo o teste de logrank o mais conhecido e utilizado para dados de tempo de vida. Testamos a igualdade entre as funções de sobrevivência dos grupos sexo (Masculino e Feminino), raça (autodeclarados Brancos e Não-Brancos) e PNE. Para o nível de confiança de 10%, há diferença nas curvas de sobrevivência em relação ao sexo, com indícios de que os homens apresentam um maior índice de rentenção, e em relação ao grupo de portadores de necessidades especiais.

4 Ajustando modelos paramétricos

É possível utilizar modelos probabilísticos para descrever o tempo de vida a partir dos dados observados. Neste caso, a análise foi feita para os dados dos alunos de 2015.1 sem nenhum tipo de partição. Ajustamos o modelo gama generalizado e seus modelos encaixados, isto é, modelos obtidos a partir da distribuição gama generalizada por meio de restrições nos parâmetros. A estimação dos parâmetros é feita pelo método da máxima verossimilhança. Para selecionar um destes modelos, utiliza-se a razão de verossimilhança para testar a hipótese de que os modelos são adequados. A estatística de teste é dada por

\[TRV = -2 \log{\left[\frac{L(\hat{\theta}_M)}{L(\hat{\theta}_G)} \right]} = 2[\log{L(\hat{\theta}_G)} - \log{L(\hat{\theta}_M)}]\]

em que \(L(\hat{\theta}_G)\) é o logaritmo da função de verossimilhança do modelo generalizado e \(L(\hat{\theta}_M)\) é o logaritmo da função de verossimilhança do modelo avaliado. que sob \(H_0\) tem distribuição \(\chi^2\) com graus de liberdade igual a diferença do número de parâmetros \((\hat{\theta}_G e \hat{\theta}_M)\) dos modelos sendo comparados. Portanto, utilizamos o modelo gama generalizado como referência, devemos escolher o modelo em que o valor-p do teste da razão de verossimilhança indique que não há diferença estatística entre o modelo avaliado e o modelo de referência.

A seguir, são apresentados os modelos, por meio da sua função densidade de probabilidade, as estimativas dos parâmetros e as função de sobrevivência associada. Em seguida, apresentam-se os resultados do teste da razão de verossimilhança, que indicaram o modelo gama como o mais apropriado. Finalmente, é apresentado o gráfico do modelo.

4.1 Modelo Exponencial

  • Função densidade de Probabilidade:

\[f(t) = \frac{1}{\alpha} exp\left\{-\left(\frac{t}{\alpha}\right)\right\}\]

  • Função de sobrevivência:

\[S(t) = \exp\left\{-\left(\frac{t}{\alpha} \right)\right\}\]

## Call:
## flexsurvreg(formula = Surv(tempo, censura) ~ 1, data = surv_data, 
##     dist = "exponential")
## 
## Estimates: 
##       est       L95%      U95%      se      
## rate  4.21e-04  3.64e-04  4.87e-04  3.13e-05
## 
## N = 457,  Events: 181,  Censored: 276
## Total time at risk: 430046
## Log-likelihood = -1587.94, df = 1
## AIC = 3177.88

4.2 Modelo Lognormal

  • Função densidade de probabilidade:

\[f(t) = \frac{1}{\sqrt{2\pi t\sigma}} exp\left\{ -\frac{1}{2} \left(\frac{\log t - \mu}{\sigma} \right)^2 \right\}\]

  • Função de sobrevivência:

\[S(t) = \Phi\left(\frac{-\log{t}+\mu}{\sigma}\right)\]

## Call:
## flexsurvreg(formula = Surv(tempo, censura) ~ 1, data = surv_data, 
##     dist = "lognorm")
## 
## Estimates: 
##          est     L95%    U95%    se    
## meanlog  7.1936  7.1630  7.2242  0.0156
## sdlog    0.2309  0.2086  0.2556  0.0120
## 
## N = 457,  Events: 181,  Censored: 276
## Total time at risk: 430046
## Log-likelihood = -1332.319, df = 2
## AIC = 2668.637

4.3 Modelo Weibull

  • Funcão densidade de probabilidade:

\[f(t) = \frac{\gamma}{\alpha^\gamma}t^{\gamma-1}exp\left\{ - \left(\frac{t}{\alpha}\right)^\gamma \right\}\]

  • Função de sobrevivência:

\[S(t) = \exp\left\{-\left(\frac{t}{\alpha} \right)^\gamma\right\}\]

## Call:
## flexsurvreg(formula = Surv(tempo, censura) ~ 1, data = surv_data, 
##     dist = "weibull")
## 
## Estimates: 
##        est       L95%      U95%      se      
## shape     5.412     4.884     5.997     0.284
## scale  1464.879  1425.957  1504.865    20.128
## 
## N = 457,  Events: 181,  Censored: 276
## Total time at risk: 430046
## Log-likelihood = -1335.912, df = 2
## AIC = 2675.825

4.4 Modelo Gama

  • Função densidade de probabilidade:

\[f(t) = \frac{1}{\Gamma(k)\alpha^k}t^{k-1}\exp\left\{-\left(\frac{t}{\alpha}\right)\right\}, \,\,\, t >0.\] em que \(\Gamma(k) = \int_0^\infty x^{k-1}\exp\{-x\}dx\) é a função gama.

  • Função de sobrevivência

\[S(t) = \int_t^\infty \frac{1}{\Gamma(k)\alpha^k}u^{k-1}\exp\left\{-\left(\frac{u}{\alpha}\right)\right\}du\]

## Call:
## flexsurvreg(formula = Surv(tempo, censura) ~ 1, data = surv_data, 
##     dist = "gamma")
## 
## Estimates: 
##        est       L95%      U95%      se      
## shape  20.90986  17.18076  25.44835   2.09562
## rate    0.01539   0.01253   0.01891   0.00161
## 
## N = 457,  Events: 181,  Censored: 276
## Total time at risk: 430046
## Log-likelihood = -1328.883, df = 2
## AIC = 2661.766

4.5 Gama generalizada

\[f(t) = \frac{\gamma}{\Gamma(k)\alpha^{\gamma k}} t^{\gamma k - 1} \exp\left\{- \left(\frac{t}{\alpha}\right)^\gamma \right\}\]

## Call:
## flexsurvreg(formula = Surv(tempo, censura) ~ 1, data = surv_data, 
##     dist = "gengamma")
## 
## Estimates: 
##        est     L95%    U95%    se    
## mu     7.2307  7.1924  7.2690  0.0195
## sigma  0.2087  0.1854  0.2350  0.0126
## Q      0.3968  0.1265  0.6672  0.1379
## 
## N = 457,  Events: 181,  Censored: 276
## Total time at risk: 430046
## Log-likelihood = -1328.021, df = 3
## AIC = 2662.042

5 Modelos com regressores

## Call:
## flexsurvreg(formula = Surv(tempo, censura) ~ raça + sexo + pne + 
##     classe, data = surv_data, dist = "exponential")
## 
## Estimates: 
##                 data mean  est        L95%       U95%       se       
## rate                   NA   0.000467   0.000246   0.000884   0.000152
## raçaNao_branco   0.566832  -0.101280  -0.394907   0.192348   0.149813
## sexoM            0.705446  -0.303415  -0.605890  -0.000941   0.154326
## pneSIM           0.034653  -1.787876  -3.755509   0.179756   1.003913
## classeDE         0.878713   0.243623  -0.366308   0.853553   0.311195
##                 exp(est)   L95%       U95%     
## rate                   NA         NA         NA
## raçaNao_branco   0.903680   0.673743   1.212092
## sexoM            0.738292   0.545589   0.999059
## pneSIM           0.167315   0.023389   1.196926
## classeDE         1.275863   0.693289   2.347975
## 
## N = 404,  Events: 181,  Censored: 223
## Total time at risk: 409557
## Log-likelihood = -1573.651, df = 5
## AIC = 3157.302
## Call:
## flexsurvreg(formula = Surv(tempo, censura) ~ raça + sexo + pne + 
##     classe, data = surv_data, dist = "lognorm")
## 
## Estimates: 
##                 data mean  est       L95%      U95%      se        exp(est)
## meanlog               NA    7.15629   7.02582   7.28677   0.06657        NA
## sdlog                 NA    0.22740   0.20539   0.25178   0.01181        NA
## raçaNao_branco   0.56683    0.01501  -0.04443   0.07446   0.03033   1.01513
## sexoM            0.70545    0.04121  -0.02129   0.10371   0.03189   1.04207
## pneSIM           0.03465    0.41941   0.11558   0.72325   0.15502   1.52107
## classeDE         0.87871   -0.00438  -0.12818   0.11941   0.06316   0.99563
##                 L95%      U95%    
## meanlog               NA        NA
## sdlog                 NA        NA
## raçaNao_branco   0.95654   1.07730
## sexoM            0.97893   1.10928
## pneSIM           1.12252   2.06113
## classeDE         0.87970   1.12683
## 
## N = 404,  Events: 181,  Censored: 223
## Total time at risk: 409557
## Log-likelihood = -1327.251, df = 6
## AIC = 2666.501
## Call:
## flexsurvreg(formula = Surv(tempo, censura) ~ raça + sexo + pne + 
##     classe, data = surv_data, dist = "gamma")
## 
## Estimates: 
##                 data mean  est       L95%      U95%      se        exp(est)
## shape                 NA   21.59657  17.73139  26.30430   2.17290        NA
## rate                  NA    0.01674   0.01321   0.02121   0.00202        NA
## raçaNao_branco   0.56683   -0.00794  -0.06556   0.04967   0.02940   0.99209
## sexoM            0.70545   -0.05137  -0.11170   0.00897   0.03078   0.94993
## pneSIM           0.03465   -0.41273  -0.71878  -0.10667   0.15615   0.66184
## classeDE         0.87871   -0.00723  -0.12650   0.11204   0.06085   0.99280
##                 L95%      U95%    
## shape                 NA        NA
## rate                  NA        NA
## raçaNao_branco   0.93655   1.05093
## sexoM            0.89431   1.00901
## pneSIM           0.48735   0.89882
## classeDE         0.88117   1.11856
## 
## N = 404,  Events: 181,  Censored: 223
## Total time at risk: 409557
## Log-likelihood = -1322.982, df = 6
## AIC = 2657.963
## Call:
## flexsurvreg(formula = Surv(tempo, censura) ~ raça + sexo + pne + 
##     classe, data = surv_data, dist = "gengamma")
## 
## Estimates: 
##                 data mean  est        L95%       U95%       se       
## mu                     NA   7.172357   7.052544   7.292169   0.061130
## sigma                  NA   0.200535   0.177921   0.226024   0.012242
## Q                      NA   0.472570   0.200648   0.744493   0.138739
## raçaNao_branco   0.566832   0.000648  -0.055223   0.056519   0.028506
## sexoM            0.705446   0.060235   0.001954   0.118516   0.029736
## pneSIM           0.034653   0.406926   0.094413   0.719439   0.159448
## classeDE         0.878713   0.019689  -0.095134   0.134512   0.058584
##                 exp(est)   L95%       U95%     
## mu                     NA         NA         NA
## sigma                  NA         NA         NA
## Q                      NA         NA         NA
## raçaNao_branco   1.000648   0.946274   1.058147
## sexoM            1.062086   1.001956   1.125825
## pneSIM           1.502193   1.099013   2.053281
## classeDE         1.019884   0.909252   1.143978
## 
## N = 404,  Events: 181,  Censored: 223
## Total time at risk: 409557
## Log-likelihood = -1321.162, df = 7
## AIC = 2656.324
##              Modelo Verossimilhanca        TRV valor_p
## 1 Gama Generalizado       -1321.162   0.000000    1.00
## 2       Exponencial       -1573.651 504.978103    0.00
## 3        Log-Normal       -1327.251  12.177167    0.00
## 4           Weibull       -1335.912  29.501028    0.00
## 5             gamma       -1322.982   3.639603    0.16

6 Modelo de Cox

O modelo de Cox é amplamente utilizada no ajuste do tempo de resposta até a ocorrência de um evento utilizando covariáveis. A principal suposição do modelo é a de que as funções taxa de risco dos grupos definidos pelas covariáveis são proporcionais. Supondo que temos dois grupos (0 e 1), temos a proporção entre as funções taxa de falha constante ao longo do tempo:

\[\frac{\lambda_1(t)}{\lambda_0(t)} = K\] Definindo \(X\) como variável indicadora dos grupos, temos:

\[\begin{equation} X = \begin{cases} 1 & \text{se grupo $0$}\\ 2 & \text{se grupo $1$}. \end{cases} \end{equation}\]

Além disso, definimos \(K = \exp\{\beta x\}\). Daí,

\[\lambda(t) = \lambda_0(t)\exp\{\beta x\}\]

ou seja, a função taxa de falha assume a seguinte forma:

\[\begin{equation} \lambda(t) = \begin{cases} \lambda_1(t) = \lambda_0(t)\exp\{\beta\} & \text{se $x=1$}\\ \lambda_0(t) & \text{se $x=0$} \end{cases} \end{equation}\]

Considerando agora o caso multivariado, suponha que tenhamos \(p\) covariáveis, representadas pelo vetor \(\mathbf{x} = (x_1, \ldots, x_p)\). Uma função \(g(.)\) deve ser especificada, com a condição de \(g(0) = 1\).

\[\lambda(t) = \lambda_0\,g(\mathbf{x}^\top\mathbf{\beta})\] A função mais utilizada, neste caso é a exponencial, conforme equação abaixo. O ajuste se dá pela estimação dos \(p\) betas.

\[g(\mathbf{x}^\top\mathbf{\beta}) = \exp{(\mathbf{x}^\top\mathbf{\beta})} = \exp{(\beta_1 x_1 + \ldots + \beta_p x_p)}\]

O modelo de Cox também é chamado de modelo de riscos proporcionais, uma vez que a razão entre as funções taxa de falha se mantém constante, em outras palavras, não depende do tempo:

\[\frac{\lambda_i(t)}{\lambda_j(t)} = \frac{\lambda_0(t)\exp{(\mathbf{x}_i^\top\mathbf{\beta})}}{\lambda_0(t)\exp{(\mathbf{x}_j^\top\mathbf{\beta})}} = \exp{(\mathbf{x}_i^\top\mathbf{\beta} - \mathbf{x}_j^\top\mathbf{\beta})},\]

A suposição de que a razão das taxas de falha são constates é central no modelo de Cox. Portanto, é necessário, após o ajuste, verificar se a suposição está sendo cumprida. Caso as taxas de falha não sejam proporcionais, a suposição é violada e o modelo torna-se inadequado para descrever os dados. Dentre os métodos de verificação da suposição de riscos proporcionais encontram-se os métodos gráficos e um teste definido por Cox.

A partir da expressão acima, observa-se que os coeficientes do modelo de Cox atuam acelerando ou desacelerando a função de risco. Os coeficientes devem ser interpretados à luz da propriedade de riscos proporcionais do modelo. A interpretação se dá como a razão de riscos ou risco relativo instantâneo no momento \(t\) para uma covariável, considerando todas as outras covariáveis como constantes. Por exemplo, suponha dois indivíduos (\(i\) e \(j\)) com todos os valores observados de covariáveis iguais, exceto pela variável \(l\), daí

\[\frac{\lambda_i(t)}{\lambda_j(t)} = \exp\left\{\beta_l(x_{il}-x_{jl})\right\}\]

No caso aqui tratado, os indivíduos seriam dois alunos de C&t com todos os valores de covariáveis iguais exceto, por exemplo, sexo. Então o risco de “morte” entre os alunos do sexo masculino seria \(\exp{(\beta_l)}\) vezes o risco de estudantes do sexo feminino, mantidas fixas as outras covariáveis. Desta forma, podemos ter ideia do efeito individual de cada covariável no tempo de resposta até o evento de interesse.

Abaixo apresentamos o ajuste do modelo de Cox, o teste para verificar a suposição de riscos proporcionais e os gráficos de resíduos do modelo.

## Call:
## coxph(formula = Surv(tempo, censura) ~ raça + sexo + pne + classe, 
##     data = surv_data)
## 
##   n= 404, number of events= 181 
##    (53 observations deleted due to missingness)
## 
##                   coef exp(coef) se(coef)      z Pr(>|z|)  
## raçaNao_branco -0.0149    0.9852   0.1522 -0.098   0.9221  
## sexoM          -0.3473    0.7066   0.1579 -2.199   0.0279 *
## pneSIM         -2.0703    0.1262   1.0062 -2.058   0.0396 *
## classeDE       -0.1170    0.8896   0.3139 -0.373   0.7093  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## raçaNao_branco    0.9852      1.015   0.73103    1.3278
## sexoM             0.7066      1.415   0.51846    0.9629
## pneSIM            0.1262      7.927   0.01756    0.9064
## classeDE          0.8896      1.124   0.48089    1.6456
## 
## Concordance= 0.551  (se = 0.024 )
## Likelihood ratio test= 13.28  on 4 df,   p=0.01
## Wald test            = 9.12  on 4 df,   p=0.06
## Score (logrank) test = 10.57  on 4 df,   p=0.03

O ajuste mostra que as variáveis sexo e PNE foram significativas (\(p < 0.05\)). Como consequência, observe-se que o intervalo de confiança para a razão de riscos, obtida por \(\exp\{\beta_i\}\), nestas duas variáveis não inclui o valor 1.

A interpretação dos coeficientes é, portanto, a seguinte: para a variável sexo, o risco de ocorrência do evento de interesse (integralização) para o sexo masculino é 0.7066 vezes o risco para o sexo feminino. Em outras palavras, os estudantes do sexo masculino levam mais tempo para se formar. Para a variável PNE, o risco de ocorrência do evento de interesse para os portadores de deficiência é 0.1262 vezes o risco dos estudantes não portadores de deficiência. Estas duas informações mostram que há evidência de que estes fatores (sexo e PNE) contribuem para a retenção no Bacharelado em Ciência e Tecnologia no semestre 2015.1.

Para verificar a suposição do modelo, executa-se o teste de hipóteses de Cox para verificar a existência de riscos proporcionais. Os valores-p altos mostram que não há coeficientes dependentes do tempo no modelo de regressão, portanto, a suposição de riscos proporcionais é verificada.

O resultado acima pode ser visualizado nos gráficos dos resíduos de Schoenfeld contra o tempo, onde se não são observadas tendências.

7 Conclusão

Neste trabalho, utilizamos elementos de Análise de Sobrevivência para modelar o tempo, em dias, necessário para os alunos integralizarem a carga horária necessária para colação de grau. O objetivo é entender quais fatores atuam no prolongamento do tempo para integralização, fornecendo elementos para explicar a retenção estudantil. A evasão dos alunos foi considerado como uma censura nos dados. Para comparação, foi selecionada a turma de 2015.1 e algumas das variáveis foram binarizadas para facilitar a interpretação dos coeficientes,

Os resultados mostram que, dentre as variáveis incluídas no modelo, o sexo e a presença de necessidades especiais são significativas no sentido de aumentar o tempo necessário para integralização. Em termos de razão de riscos, foram obtidas evidências de que os estudantes do sexo masculino e portadores de deficiência apresentam maior retenção.