Variável resposta

Vejamos o que acontece ao analisar o tempo de sobrevivência de cada jogador

A seguir, a distribuição da variável resposta player_survive_time :

plot_grid(pubg_tpp1 %>% 
            ggplot(aes(x=player_survive_time))+
            geom_histogram(aes(y = ..density..), bins = 30, fill="white", color="black")+
            geom_density(alpha=.2, fill="white")+
            scale_x_continuous(labels = scales::comma, limits = c(0,40), breaks = seq(0,40,5))+
            labs(x="",y="", title = "Tempo de sobrevivência dos jogadores selecionados")
          ,
          pubg_tpp1 %>% 
            ggplot(aes(x=" ", y=player_survive_time))+
            geom_boxplot()+
            labs(x="")+
            coord_flip()
          ,
          ncol = 1, nrow = 2, align = "v", rel_heights = c(3,1))

Note que possue uma assimetria positiva

Data Wrangling

Primeiramente, vejamos as variáveis se relacionam entre si e com a variável resposta com os coeficientes de correlação de Pearson:

# Correlations
pubg_tpp1 %>% 
  select_if(is.numeric) %>% 
  cor() %>% 
  corrplot::corrplot(method = "number",type = "upper",diag = F, order = "hclust",number.cex = 0.7, title = "Correlation correlated numerics", mar=c(0,0,1,0))

É possível notar que apenas a variável player_assists não correlaciona-se com a variável resposta nem com as demais variáveis e player_dmg e player_kills são fortemente correlacionadas, isso indica que pode ser interessante remover uma delas ou juntar toda essa informação em uma única variável, veremos…

Além disso nota-se que a distância percorrida a pé é fortemente correlacionada com a variável resposta enquanto que a distância de quem andou de carro não é tão correlacionada. Uma transformação na variável player_dist_ride para uma dummy drive indicando se o indivíduo dirigiu ou não pode representar melhor esta informação.

Vejamos algumas características peculiares:

pubg_tpp1 %>% 
  select(player_kills, player_dist_ride, player_assists) %>% 
  map_dfr(~quantile(.x,  probs = seq(0,1,0.25)) %>% round(2)) %>% 
  t  %>% tidy() %>% 
  `colnames<-`(c("variável",percent(seq(0,1,0.25)))) %>% 
  kable2()

Praticamente metade da amostra não registrou abates nem possui marcação de player_dist_ride. Como a variável player_dmg apresentou correlação com a variável resposta player_survive_time, vamos fazer algumas transformações:

1. Criar uma variável dummy drive se jogador usou carro
2. Somar a player_dist_ride e player_dist_walk em uma única variável: player_dist
3. Juntar player_kills, player_dmg e player_assists em uma única variável: player_performance

performance <- tibble(w_player_kills = log(pubg_tpp1$player_kills+0.5),
                      w_player_dmg = log(pubg_tpp1$player_dmg),
                      w_player_assists = pubg_tpp1$player_assists) %>% 
  mutate(player_performance = log(w_player_dmg + w_player_assists + w_player_kills))

grid.arrange(
  performance %>% 
    select(-player_performance) %>% 
    tidyr::gather() %>% 
    ggplot(aes(x=value))+
    geom_histogram(aes(y = ..density..), bins = 30, fill="white", color="black")+
    geom_density(alpha=.2, fill="white")+
    scale_x_continuous(labels = scales::comma, limits = c(-1.5,8), breaks = seq(-1,8,1))+
    labs(x="", y="")+
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank())+
    facet_wrap(~key, scales = "free")
  ,
  performance %>% 
    select(player_performance) %>% 
    tidyr::gather() %>% 
    ggplot(aes(x=value))+
    geom_histogram(aes(y = ..density..), fill="white", color="black",bins = 15)+
    geom_density(alpha=.2, fill="white")+
    scale_x_continuous(limits = c(-1.,2.5), breaks = seq(-1,3,0.5))+
    labs(x="", y="", title = "performance"),
  ncol=1
)

pubg_tpp1 <- 
  pubg_tpp1 %>% 
  mutate(player_dist = log(player_dist_ride + player_dist_walk)) %>%  
  mutate(player_assists_d = if_else(player_assists ==0, 0, 1)) %>% 
  mutate(player_performance = performance$player_performance )%>% 
  mutate(drive = ifelse(player_dist_ride==0, "no", "yes") %>% as.factor()) %>% 
  mutate(player_kills_d = ifelse(player_kills==0, "no", "yes") %>% as.factor()) 

g1 <- 
  pubg_tpp1 %>% 
  # select_if(~ !length(table(.x))==2 & is.numeric(.x)) %>% colnames() %>% 
  select(player_survive_time,player_performance,player_dist) %>% colnames() %>% 
  map2(c("Densidade", "", ""),
       ~ plot_grid(
         pubg_tpp1 %>% 
           ggplot(aes_string(x=.x)) + 
           geom_histogram(aes(y=..density..),colour="black", fill="white", bins = 15) +
           geom_density(alpha=.2, fill="lightgrey") +
           scale_x_continuous()+
           ggtitle(.x)+
           labs(x="", y=.y)+
           theme(axis.title.x=element_blank(),
                 axis.text.x=element_blank(),
                 axis.ticks.x=element_blank())
         ,
         pubg_tpp1 %>% 
           ggplot(aes_string(, y=.x))+
           geom_boxplot(aes(x=" "))+
           labs(x="", y="")+
           coord_flip()+
           theme(axis.title.x=element_blank(),
                 axis.text.x=element_blank(),
                 axis.ticks.x=element_blank()),
         
         ncol = 1, nrow = 2, align = "v", rel_heights = c(3,1)
       )
  )

dat <- 
  pubg_tpp1 %>% 
  select_if(~.x %>% table %>% length == 2) %>% 
  mutate_at(2,~if_else(.x==0, "no", "yes")) %>% 
  .[,-1]

g2 <- map2(colnames(dat),
           c( "Porcentagem", "",""),
           ~ dat[,.x] %>% 
             tidyr::gather() %>% 
             group_by(key, value) %>% 
             summarise(n = n()) %>% 
             mutate(prop = n/sum(n)) %>% 
             ggplot(aes(x = key, y = prop,fill = value)) + 
             geom_bar(position = "fill",stat = "identity", alpha=0.7) +
             scale_y_continuous(labels = percent_format())+
             labs(x="", y = .y)+
             scale_fill_manual(values = c("grey", "#FCC14B"), name = "Legenda:")
)

grid.arrange(g1[[1]], g1[[2]], g1[[3]],g2[[1]], g2[[2]], g2[[3]], ncol=3, heights=c(3/5, 2/5))

pubg_tpp1 %>% 
  select_if(~ is.numeric(.x) & length(table(.x))!=2) %>% 
  sumario_custom_num() %>% 
  kable2()

grafico_descritivo(x = pubg_tpp1,
                   colNames = c('player_survive_time', "player_performance", 'player_dist',
                                'player_assists_d',"player_kills_d", 'drive'),
                   color='drive',
                   colors = c("grey", "#FCC14B"))

Kaplan-Meier

Para isso existem algumas alternativas como o estimador de Kaplan-Meier, que utiliza os conceitos de independência e de probabilidade condicional para deduzir a probabilidade de sobreviver até o tempo \(t\).

Veja a seguir são ajustados os modelos univariados de Kaplan-Meier para cada uma das coivaráveis da amostra:

surv <- Surv(pubg_tpp1$player_survive_time)
resultado_km <-
  list(geral            = survfit(surv ~ 1 ,data = pubg_tpp1),
       player_assists_d = survfit(surv ~ player_assists_d ,data = pubg_tpp1),
       drive            = survfit(surv ~ drive,data = pubg_tpp1 ),
       player_kills_d   = survfit(surv ~ player_kills_d,data = pubg_tpp1))

Veja os resultados da função de sobrevivência sem levar em consideração nenhuma das coivaráveis:

surv_summary(resultado_km[[1]], pubg_tpp1) %>% .[1:5,-ncol(.)] %>% cbind(variable = "Geral") %>% 
  select(variable, everything())%>%
  kable2()

A função surv_summary() retorna um quadro de dados com as seguintes colunas:

• time: o tempo em que a curva tem um passo.
• n.risk: o número de sujeitos em risco em t.
• n.evento: o número de eventos que ocorrem no tempo t.
• n.censor: número de eventos censurados.
• surv: estimativa da probabilidade de sobrevivência.
• std.err: erro padrão de sobrevivência.
• superior: extremidade superior do intervalo de confiança
• inferior: extremidade inferior do intervalo de confiança
• estratos: indica a estratificação da estimativa de curvas. Os níveis de estratos (um fator) são os rótulos das curvas (se houver).

resultado_log_rank <- 
  c(geral = "",
    player_assists_d=round(1-pchisq(survdiff(surv~player_assists_d,data = pubg_tpp1)$chisq,1),5),
    drive=round(1-pchisq(survdiff(surv~drive,data=pubg_tpp1)$chisq,1),5),
    player_kills_d=round(1-pchisq(survdiff(surv~player_kills_d,data=pubg_tpp1)$chisq,1),5)
  )

survplot <- map2(resultado_km,
                 case_when(resultado_log_rank == '0' ~ "log-rank: \n p < 0,00001",
                           resultado_log_rank == "" ~ "log-rank não se aplica",
                           resultado_log_rank != '0' | resultado_log_rank != '' ~ 
                             paste0("log-rank: \n p =",as.numeric(resultado_log_rank))),
                 ~ autoplot(.x)+
                   ggtitle(stringr::str_remove_all(names(.x$strata)[1],"(=no|=yes)"))+
                   annotate("label",y = 0.20, x = 5,
                            label = .y,
                            size = 4, colour = "red",hjust=0.1)+ 
                   scale_fill_manual(values = c("grey", "#FCC14B"))+
                   scale_color_manual(values = c("grey", "#FCC14B"))+
                   theme(legend.position = c(0.85,0.7))+
                   scale_x_continuous(limits = c(0,30), breaks = seq(0,30,5))
                 
)
grid.arrange(survplot[[1]], survplot[[2]] ,survplot[[3]], survplot[[4]], ncol=2)

Função de risco (hazard) ou taxa de falha

Função de risco (hazard) ou taxa de falha é o risco “instantâneo” denotada por \(\lambda(t)\) é uma taxa, não uma probabilidade e pode assumir qualquer valor real maior que zero.

No exemplo representa a taxa de incidência ou risco acumulado para um indivíduo morrer até o momento \(t\), dado que sobreviveu até este momento. É muito informativa quando comparada com a função de sobrevivência pois diferentes \(S(t)\) podem ter formas semelhantes, enquanto que respectivas \(\lambda(t)\) podem diferir drasticamente.

survplot <-
  map(resultado_km  ,
      ~ ggsurvplot(.x, conf.int = TRUE, 
                   palette = c("grey", "#FCC14B"),
                   risk.table = F,break.time.by = 5,
                   fun = "cumhaz",title = stringr::str_remove_all(names(.x$strata)[1],"(=no|=yes)"))
  )
arrange_ggsurvplots(survplot, print = TRUE,
                    ncol = 2, nrow = 2)

O risco cumulativo \(H( t)\) pode ser interpretado como a força cumulativa da mortalidade. Em outras palavras, corresponde ao número de eventos que seriam esperados para cada indivíduo pelo tempo t se o evento fosse um processo repetitivo.

Modelo de cox

É caracterizado pela presença dos coeficientes \(\beta\)s que medem os efeitos (semelhantes à análise de regressão logística múltipla e linear múltipla) das variáveis explicativas sobre a função de risco. Em um modelo de regressão de riscos proporcionais de Cox, a medida do efeito é a taxa de risco, que é o risco de falha, dado que o participante sobreviveu até um tempo específico.

Algumas das suposições para o correto uso do modelo de regressão de riscos proporcionais de Co incluem:

• independência dos tempos de sobrevivência entre indivíduos distintos na amostra,
• relação multiplicativa entre os preditores e o risco,
• uma taxa de risco constante ao longo do tempo.

O modelo de riscos proporcionais de Cox é chamado de modelo semi-paramétrico , porque não há suposições sobre o formato da função de risco de linha de base. No entanto, existem outras suposições, como observado acima.

É possível utilizar as estatísticas de Wald, da razão de verossimilhança e escore para fazer inferências sobre os parâmetros do modelo

Veja a seguir a significância dos coeficiente estimado em modelos univariados para cada variável candidata ao modelo:

# Modelos univariados
covariates    <- c("player_kills","player_dist_ride","player_performance",
                   "player_dist_walk","player_dmg", "player_dist", 
                   "player_assists_d","drive", "player_kills_d")
univ_formulas <- map(covariates,~ as.formula(paste('Surv(player_survive_time) ~', .x)))
univ_models   <- map( univ_formulas, ~coxph(.x, data = pubg_tpp1))

# estrair resultados 
map2_df(univ_models,
        covariates,
        function(x,y){ 
          x                = summary(x)
          p.value          = signif(x$wald["pvalue"], digits=2)
          wald.test        = signif(x$wald["test"], digits=2)
          beta             = signif(x$coef[1], digits=2);#coeficient beta
          HR               = signif(x$coef[2], digits=2);#exp(beta)
          HR.confint.lower = signif(x$conf.int[,"lower .95"], 2)
          HR.confint.upper = signif(x$conf.int[,"upper .95"],2)
          HR               = paste0(HR, " (", HR.confint.lower, "-", HR.confint.upper, ")")
          res              = tibble(y,beta, HR, wald.test, p.value)
          colnames(res)    = c("covariates","beta", "HR (95% CI for HR)", "wald.test", "p.value")
          res
        }) %>% 
  kable2(linhas = 7)

Modelo de Cox usando uma variável categórica retorna uma razão de risco, que, acima de 1 indica uma covariável que está positivamente associada à probabilidade do evento e, portanto, negativamente associada ao tempo de sobrevida. O oposto vale para HR menor que um e HR = 1 indica que a covariável não tem efeito.

final_model  <- 
  coxph(Surv(player_survive_time) ~ player_performance+player_dist+drive,
        data = pubg_tpp1,x=T,method="breslow")

summary(final_model)

## Call:
## coxph(formula = Surv(player_survive_time) ~ player_performance + 
##     player_dist + drive, data = pubg_tpp1, x = T, method = "breslow")
## 
##   n= 200, number of events= 200 
## 
##                       coef exp(coef) se(coef)       z Pr(>|z|)    
## player_performance -0.7469    0.4738   0.1787  -4.179 2.92e-05 ***
## player_dist        -1.7599    0.1721   0.1150 -15.307  < 2e-16 ***
## driveyes            0.8832    2.4186   0.2091   4.225 2.39e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## player_performance    0.4738     2.1105    0.3338    0.6726
## player_dist           0.1721     5.8117    0.1374    0.2156
## driveyes              2.4186     0.4135    1.6055    3.6434
## 
## Concordance= 0.883  (se = 0.008 )
## Likelihood ratio test= 362.7  on 3 df,   p=<2e-16
## Wald test            = 274  on 3 df,   p=<2e-16
## Score (logrank) test = 407.8  on 3 df,   p=<2e-16

No modelo ajustado note-se que existe uma associação negativa entre player_performance e mortalidade e entre player_dist e mortalidade (ou seja, o risco de morte diminui para jogadores que percorrem maiores distâncias e possuem melhor performance).

As estimativas dos parâmetros representam o aumento no log esperado do risco relativo para cada aumento de uma unidade no preditor, mantendo os outros preditores constantes.

Para interpretabilidade, calcularemos as taxas de risco exponenciando das estimativas dos parâmetros. Para a player_performance, exp(-0.7469196)= 0.4738239. Isso implica que diminui para 47.38% do valor original do risco esperado em relação a um aumento de uma unidade na performance, mantendo as demais variáveis constantes. A interpretação de player_dist em escala logarítimica é feita de maneira semelhante.`

Já para os jogadores onde drive = 1 (que dirigiram durante a partida) existe uma relação positiva, como exp(0.8831835)= 2.4185871. O risco esperado corresponde à 2.4185871 do valor original nos que dirigiram em comparação aos que não dirigiram, mantendo as demais variáveis constantes.

map2_df(1:3,final_model$coefficients %>% names(),~
          tibble(
            variable = .y,
            beta             = signif(summary(final_model)$coef[.x,1], digits=2), #coeficient beta
            HR               = signif(summary(final_model)$coef[.x,2], digits=2), #exp(beta)
            HR.confint.lower = signif(summary(final_model)$conf.int[.x,"lower .95"], 2),
            HR.confint.upper = signif(summary(final_model)$conf.int[.x,"upper .95"],2)) %>% 
          mutate(HR= paste0(HR, " (", HR.confint.lower, "-", HR.confint.upper, ")")
          )
) %>% kable2()

res <- 
  tibble(residuo_deviance = resid(final_model,type="deviance") ,
         residuo_martingal = resid(final_model,type="martingal"),
         linear_predictors = final_model$linear.predictors)

# Graficos:
grid.arrange(
  ggplot(res, aes(x=linear_predictors, y=residuo_martingal))+ geom_point()+geom_hline(yintercept=0, color='coral')+ylab("Resíduos Martingual"),
  ggplot(res, aes(x=linear_predictors, y=residuo_deviance))+ geom_point()+geom_hline(yintercept=0, color='coral')+ylab("Deviance"),
  ncol=2
)

final_model %>% cox.zph %>% ggcoxzph

Considerações finais

Como era de se esperar, o risco de ser abatido diminui para jogadores que possuem melhor performance e também para os jogadores que percorrem maiores distâncias (o que mostra que ficar parado no jogo em uma zona pode não ser a melhor ideia, já é quanto mais se movimenta maior a quantidade de itens que podem ser coletados).

Interessante notar que a curva de sobrevivência para os jogadores que dirigiram apresenta resultado oposto ao risco esperado nos que dirigiram, isso ocorre pois esses dois modelos calculam medidas diferentes.