R - How to calculate the density from the accumulated curve based on the Cox regression model

Asked

Viewed 410 times

5

I calculated the accumulated probability (cdf) of my data, based on the probability of exceedance (Edf), using the Cox regression model. So far ok, no problem at all.

However, does anyone know if there is a command to turn this data into probability density (pdf)?

I have tested using the histogram function, but it does not work properly.

dado1<-c(128.1072, 124.2218, 127.5064, 143.5201, 121.6476, 121.4071, 133.5725, 127.9324, 115.7151, 131.6176, 113.7500, 122.2064, 133.9970, 125.4781, 122.9766, 132.7081, 124.9619, 134.4549, 127.4127, 121.9021, 111.9924, 122.4483, 132.1261, 129.7735,124.7136, 118.2293, 120.5072, 129.5527, 125.7787)

dado2<-c(174.07874, 132.74495, 84.52224, 82.93248, 113.13792, 112.87297, 163.48032, 170.10432, 184.41215, 169.30945, 152.35201, 127.44576, 130.62528, 123.20640, 59.61600, 48.75264, 77.10335, 113.93281, 83.99231, 164.27521, 111.81314, 72.06912, 169.04448, 229.45537, 79.48800,  57.23136,  72.33408, 95.38560, 136.18944)

dado3<-c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)

dado4<-c(1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010)

dados<-data.frame(cbind(dado4,dado2,dado1,dado3))

require(survival)

curva <- coxph(Surv(dados[,2], dados[,4]) ~ dados[,3], dados)
a<-summary(curva)
coef<-as.numeric(data.frame(a$coef[1]))

edf<-survfit(curva)$surv
edf<-append(1,edf)
cdf<-1-edf

ano<-12
prevcox<-edf^exp(coef*dados[ano,3])

I want to find the pdf, based on the cdf and Edf found, for the prevcox variable.

r
  • 2

    For empirical densities perhaps core densities, which can be calculated with density(x). See the help page ?density for the various available cores. To obtain the values of x and of y can be with pdf <- density(x); pdf$x; pdf$y.

  • Thanks @Ruibarradas! However, when I went to test your suggestion, you presented an error. "pdf$x error: closure object cannot be divided into subsets". Any suggestions?

  • It must be because pdf is already the name of a base R function and sometimes there are name conflicts. Just try f, thus: f <- density(x); f$x; f$y.

  • Now yes, thank you @Ruibarradas!

  • because you are using survival analysis if your data presents no censorship?

  • Are you sure you want to find the pdf based on cdf? For in survival models, survival functions and accumulated risk are usually used to interpret the model

Show 1 more comment

1 answer

1

If such values are the result of your F(x), then you can approximate your fdp(x) from the following equation:

fórmula

where x2 and X1 are the values of "x" you used to calculate your F(x), where x2>X1, but the difference between x2 and X1 must be close to zero, so you can assume, for example that the above equation results in fdp(x2).

In code R, you can do something like this:

F_para_fdp <- function(x1, x2, Fx2, Fx1){
  return((Fx2-Fx1)/(x2-x1))
}

  • Thanks @Guilhermeparreira! However, I don’t think I could make it work. I indicated who would be my X1 and x2, and I circled exactly as you wrote above. Apparently it runs, but I get no results. And when I type "F_para_fdp", it brings back exactly all the lines I wrote.

  • You’re welcome, Iara, for nothing! The function I put in is just to start the process... Note that the values you put in the question are the result of F(x), not the values of "x". For me to provide you with a complete answer, I need you to provide the values of "x".

  • Guilherme Parreira, I answered here in the post to be more detailed. Thank you!

  • Hello @Iara, you should edit the text and the title of your question, not put as an answer. Because the way the question was asked, it’s pretty generic, and then my answer fits, but in your case you need it in the context of Cox’s model. Next week I will probably study such a model, and I will be able to provide an answer here for you!!

  • Thank you Guilherme Parreira! I believed I could find out that way. But I rewrote the question to make it clear.

Browser other questions tagged r

You are not signed in. Login or sign up in order to post.