Otimização de um modelo de risco proporcional de Cox usando Optimx ()

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br

[ad_1]

[This article was first published on R | Joshua Entrop, and kindly contributed to R-bloggers]. (Você pode relatar problemas sobre o conteúdo desta página aqui)


Quer compartilhar seu conteúdo em R-bloggers? clique aqui se você tiver um blog, ou aqui se não tiver.

Nesta postagem do blog, vamos otimizar um modelo de risco proporcional de Cox usando um método de estimativa de probabilidade máxima (MLE). Para isso, vamos primeiro definir a função de verossimilhança de nosso modelo de Cox e suas primeiras derivadas parciais, às vezes chamadas de função de pontuação. Mais tarde, passaremos essas funções para diferentes algoritmos de otimização usando optimx :: optimx () para obter as estimativas mais prováveis ​​dos parâmetros do nosso modelo. Você pode encontrar o R-script que usei para esta postagem como .TXT arquivo aqui.

O modelo de risco proporcional de Cox é provavelmente o modelo mais comumente aplicado para análise de sobrevivência em pesquisas epidemiológicas. Os modelos Cox oferecem uma maneira eficiente de se ajustar à escala de tempo subjacente (t ) que de outra forma requer modelos paramétricos mais complexos. Nesses modelos, por exemplo, no modelo de Poisson ou Weibull, existem duas maneiras de ajustar para a escala de tempo subjacente: qualquer um usa uma abordagem de divisão de tempo onde o tempo (t ) é dividido em vários intervalos de tempo que são adicionados como covariáveis ​​ao modelo ou um usa uma função spline para modelar o efeito da escala de tempo subjacente. Ambas as abordagens precisam de suposições sobre a forma do efeito da escala de tempo subjacente, por exemplo, os pontos de tempo usados ​​para a divisão ou os graus de liberdade da função spline. Além disso, ambas as abordagens adicionam vários parâmetros ao modelo. Ajustar esses modelos paramétricos mais complexos pode ser um desafio ao enfrentar dados limitados. Em vez disso, o modelo de Cox só precisa de uma suposição para poder se ajustar à escala de tempo subjacente, ou seja, o efeito da exposição é considerado proporcional ao longo do tempo. Essa suposição é conhecida como suposição de risco proporcional e está sujeita a muitos debates. Além disso, o modelo de Cox não precisa de nenhum parâmetro adicional para ajustar para a escala de tempo (t ). Portanto, ele também pode ser usado para analisar dados mais limitados.

Leia Também  Inoapps Anuncia Nova Divisão Global de Mudança de Empresa

Como em minhas postagens anteriores, usaremos o conjunto de dados de câncer de pulmão incluído no {sobrevivência} pacote como um exemplo. Para obter mais informações sobre este conjunto de dados, dê uma olhada no arquivo de ajuda ? sobrevivência :: pulmão Especificamente, modelaremos o efeito do sexo e da idade na sobrevida de pacientes com câncer de pulmão neste conjunto de dados. Para estimar esse efeito, usaremos um modelo de Cox.

O modelo Cox segue a forma geral de

[ lambda(t|X) = lambda_0

Em nosso caso, ajustaremos o seguinte modelo de Cox incluindo as variáveis ​​independentes femininas ((x_1) ) e idade ((x_2) ).

[ lambda(t|X) = lambda_0

Portanto, agora vamos começar a carregar o conjunto de dados e configurar as variáveis.

# 1. Prefix -------------------------------------------------------------------
# Remove all files from ls
rm(list = ls())
# Loading packages
require(survival)
require(optimx)
require(numDeriv)
require(purrr)
require(dplyr)
require(tibble)
require(broom)
# 2. Loading data set ---------------------------------------------------------
#Reading the example data set lung from the survival package
lung 

The Cox model in its basic form requires unique event time i.e. no ties, as we will see later. To deal with this assumption we will randomly add or subtract small amounts of time from each event time in the data set, as seen in the code below. In an upcoming blog post I will address different approaches to deal with ties in Cox models.

#Removes time ties in data set
set.seed(2687153)
lung$time %
count(time) %>%
arrange(desc(n)) %>%
head(5)
## time n
## 1 5.042231 1
## 2 10.917351 1
## 3 11.016322 1
## 4 11.020302 1
## 5 11.989413 1

Excelente! Não temos mais nenhum vínculo em nosso conjunto de dados. Agora podemos prosseguir com a definição da função log-verossimilhança do nosso modelo de Cox. A função de log-verossimilhança do modelo de Cox geralmente segue a forma

[ ln L(beta) = sum d_i bigg( X_i beta - ln sum_{j:t_jgeq t_i} theta_j bigg) ]

Onde ( theta = exp ( beta X) ) e (d_i ) é o indicador de evento para o (i ^ {th} ) tema. Se inserirmos nossas variáveis ​​independentes de cima, obtemos ( theta = exp ( beta_1 x_1 + beta_2 x_2) ) para o nosso caso específico.

O ( sum_ {j: t_j geq t_i} theta_j ) parte da fórmula acima torna o cálculo um pouco mais complicado, pois requer uma certa ordem das observações. Precisamos somar ( theta_j ) para todas as observações (j ) que tem um evento posterior ou exatamente ao mesmo tempo que nosso (i ^ th ) observação. Basicamente, essa é a soma cumulativa de ( theta_j ) ao longo dos tempos do evento (t ) por ordem decrescente. Para calcular essa quantidade em R, podemos classificar o conjunto de dados por tempos de eventos decrescentes e depois usar o base :: cumsum () função para calcular a soma cumulativa.

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br

Vamos juntar tudo isso e definir nossa função de log-verossimilhança em R.

# 3. Define log-likelihood function for Cox regression model ------------------
negll %
arrange(desc
mutate(thetaj = cumsum(theta))
#Estimate negative log likelihood value
val 

To improve our optimisation we should also pass the gradient functions for our model to the optimx() later. The gradient function for the Cox model in general follows

[ ln L'(beta) = sum d_i bigg(X_i - frac{sum_{j:t_jgeq t_i} theta_j X_j}{sum_{j:t_jgeq t_i} theta_j} bigg)]

In our case we yield the following two gradient functions for (beta_1) and (beta_2).

[ ln L'(beta_1) = sum d_i bigg(x_{1i} - frac{sum_{j:t_jgeq t_i} theta_j x_{1j}}{sum_{j:t_jgeq t_i} theta_j} bigg)]

[ ln L'(beta_2) = sum d_i bigg(x_{2i} - frac{sum_{j:t_jgeq t_i} theta_j x_{2j}}{sum_{j:t_jgeq t_i} theta_j} bigg)]

We can use this function to get the following gradient function for our Cox model. Note that we have to use the base::cumsum() twice in the code below to calculate the cumulative sum of (theta_j x_{1j}) and (theta_j x_{2j}).

# 4. Define gradient function for Cox regression model ------------------------
negll_grad %
arrange(desc
mutate(thetaj = cumsum(theta),
thetajx1 = cumsum(theta * x1),
thetajx2 = cumsum(theta * x2))
#Calculate partial gradient functions
gg[1] 

Lets just check if our gradient function is correct by comparing it with the approximation of the gradient function calculated with the numDerive::grad() function.

# 4.1 Compare gradient function with numeric approximation of gradient ========
# compare gradient at 1, 0, 0, 0
mygrad 

Looks like we get the same numbers and our gradient functions works fine.

Now we pass both our log-likelihood and gradient function on to our optimx() call.

# 5. Find minimum of log-likelihood function ----------------------------------
# Passing names to the values in the par vector improves readability of results
opt %
rownames_to_column("algorithm") %>%
filter(convcode != 9999) %>%
arrange(value) %>%
select(algorithm, beta_female, beta_age, value) %>%
head(7)
## algorithm beta_female beta_age value
## 1 Rcgmin -0.5133011 0.01703863 742.7912
## 2 nlminb -0.5133013 0.01703864 742.7912
## 3 nlm -0.5133009 0.01703866 742.7912
## 4 L-BFGS-B -0.5133009 0.01703860 742.7912
## 5 BFGS -0.5133055 0.01703866 742.7912
## 6 Nelder-Mead -0.5134599 0.01702275 742.7912
## 7 CG -0.4557679 0.01735647 742.8464

Seis dos algoritmos de otimização implementados no {optimx} pacote rendeu valores de máxima verossimilhança iguais até quatro decimais. Isso sugere diferenças desprezíveis entre esses modelos. Se imprimirmos mais decimais dos valores de log da verossimilhança máximos, provavelmente veremos algumas pequenas diferenças entre eles. No entanto, todos os modelos sugerem que as mulheres têm cerca de metade do risco de morrer de câncer de pulmão em comparação com os homens e o risco de morrer aumenta com o aumento da idade.

Vamos verificar quais estimativas obteríamos para sexo e idade se ajustássemos um modelo de Cox usando o sobrevivência :: coxph () função e compare-os com nossas estimativas.

# 6. Estimate regression coefficients using coxph ----------------------------
cox_model %
bind_rows() %>%
filter(!is.na(mean_dif)) %>%
mutate(mean_dif = abs(mean_dif)) %>%
arrange(mean_dif)
## opt_name mean_dif
## 1 Rcgmin 3.497854e-08
## 2 nlminb 1.140406e-07
## 3 L-BFGS-B 2.315549e-07
## 4 nlm 2.650035e-07
## 5 BFGS 4.366428e-06
## 6 Nelder-Mead 1.587632e-04
## 7 CG 5.753333e-02
## 8 Rvmmin 5.133012e-01

Podemos ver que a diferença média entre nossas estimativas e as estimativas produzidas com o sobrevivência :: coxph () modelo é negligenciável para a maioria de nossos modelos. Parece que tudo funcionou bem. No entanto, como um bom pesquisador, naturalmente também gostaríamos de obter algumas estimativas de incerteza. Então, vamos pegar o modelo que equipamos com o Rcgmin algoritmo de nosso optimx () saída e calcular o erro padrão para nossas estimativas usando a matriz de hessian. Se você estiver mais interessado no cálculo, pode dar uma olhada em meu post anterior no blog.

# 8. Estimate the standard error ----------------------------------------------
#Extract hessian matrix for the Rcgmin optimisation
hessian_m %
print()
## se_rcgmin se_coxph
## 1 0.167457336 0.167457337
## 2 0.009224444 0.009224444
all.equal(ses[,"se_rcgmin"], ses[, "se_coxph"])
## [1] TRUE

Com base no erro padrão, podemos calcular os intervalos de confiança para nossas estimativas agora.

# 9. Estimate 95%CIs using estimation of SE -----------------------------------
# Extracting estimates from the Rcgmin optimisaiton
coef_test %
round(4)
## Estimate CI_lower CI_upper se
## beta_female -0.5133 -0.8415 -0.1851 0.1675
## beta_age 0.0170 -0.0010 0.0351 0.0092

Excelente! Obtivemos nosso próprio modelo de Cox com intervalos de confiança.

Para resumir, especificamos nossa função de log-verossimilhança e sua função gradiente, e a otimizamos usando o optimx :: optimx () função. Com base na saída do optimx () chamada, fomos capazes de obter o erro padrão e os intervalos de confiança para nossas estimativas.



[ad_2]

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br