Determinar que riesgos comienza a aumentar para una variable continua

Question

Estoy interesado en una variable continua, es decir, de la presión arterial.

La mayor presión arterial, mayor es el riesgo de ataque cardíaco y accidente cerebrovascular. Sin embargo, los estudios con frecuencia informan de que también baja la presión arterial está asociada a resultados adversos.

La pregunta es: ¿cuál es el óptimo de la presión arterial? A lo que el valor de la presión arterial corre el riesgo de comenzar a aumentar?

En otras palabras, ¿cómo puedo modelar y visualizar de forma gráfica, lo que hazard ratio diversos niveles de la presión arterial se asocia con. Sospecho que algunos sugieren restringido splines cúbicos. ¿Tiene alguna sugerencia sobre R de los paquetes que me ayudará a visualizar el efecto de la presión arterial en el riesgo. Estoy bastante familiarizado con regresión de Cox y plan de uso de la RMS paquete. Las variables dependientes del tiempo son incluidos.

Datos de ejemplo (sin variables dependientes del tiempo):

event <- c(1,0,1,0,0,0,0,1,0,0,0,1,1,0,1,0,0,1,0,1,1,1,1,0,1,1,1,0,0,1)
survival <- c(4,29,24,29,29,29,29,19,29,29,29,3,9,29,15,29,29,11,29,5,13,20,22,29,16,21,9,29,29,15)
statin <- c(0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0)
bloodpressure <- c(160,120,150,140,135,110,139,140,153,129,149,163,179,129,144,119,100,115,145,150,130,120,122,129,116,171,129,126,159,150)
data <- data.frame(event, survival, statin, bloodpressure)
View(data)

require(rms)
fit <- coxph(Surv(survival, event) ~ statin + rcs(bloodpressure, 3), data=data)

He tenido algo como esto en mente: http://www.bmj.com/content/325/7372/1073/F1

http://www.nejm.org/doi/full/10.1056/NEJMoa1215740

Gracias

Answer 1

Creo que usted está en el camino correcto con el rcs de la rms paquete. De hecho, rms viene con su propia versión de Coxph, se llama cph.

Usted puede ser que desee probar el siguiente

fit <- cph (Surv(survival, event) ~ statin + rcs(bloodpressure, 3), x=TRUE, y=TRUE, data=data)
# x, y for predict, validate and calibrate;
plot(Predict(fit), data=data)

Usted puede leer más acerca de esto en la conferencia nota de la Profesora Harrell (autor de RMS paquete y un libro con el mismo nombre) http://biostat.mc.vanderbilt.edu/wiki/pub/Main/RmS/rms.pdf Desplácese a la Sección 18.4

Aquí está, sin embargo, otro aspecto a considerar - ¿se han incluido todas las variables relevantes en el modelo? No es necesario tener sólo una variable en el modelo para visualizar su efecto. Un modelo especificado produce estimaciones sesgadas porque no tiene control de todas las posibles variables de confusión. (Véase, por ejemplo: http://en.wikipedia.org/wiki/Omitted-variable_bias)

EDITADO: Para la trama.Predecir el método funcione correctamente, tendrá las siguientes dos líneas antes de la cph línea. (Tengo que confesar que no sabe el significado exacto de este, pero veo que esta en 18.1 de Harrell notas, y ayuda a resolver el mensaje de error para mí) Espero que esto ayude~

dd <- datadist(data)
options(datadist = 'dd')

Answer 2

No acepte esta respuesta, es sólo para la variedad:

En oncología, hay un enfoque diferente, más cerca de riesgos proporcionales. Funcionaría de la siguiente manera: Uno divide la presión arterial escala en pequeñas y los intervalos de las estimaciones de los cocientes de riesgo "local". Usted obtendrá una parcela de la presión arterial vs hazard ratio. Se llama STEPP por Bonetti y Gelber (2004). También hay un paquete de R. Pero hay que tener en cuenta las desventajas de este apporach: necesita bastante grandes tamaños de muestra y los resultados dependen de la (arbitraria) de la longitud de los intervalos. Después de todo, es más útil para explicativo de confirmatoric análisis. También, usted no va a obtener el óptimo de la presión arterial, pero sólo un intervalo. El único pro es que usted no tiene que especificar la forma de la función de riesgo.

(Que también es una desventaja debido a que es tentador pensar menos en el avance de su análisis)

Answer 3

Me gustaría tratar el cox.ph de la familia en la mgcv paquete, que está disponible para los modelos aditivos generalizados implementado por la gam función. Los modelos aditivos generalizados (GAMs) son una especie de regresión lineal múltiple analógico por modelos spline. Si usted tiene un gran conjunto de datos, no te preocupes, porque mgcv es también bastante rápido. Después de montar las flexibles, semi-paramétrico de GAM, puede globo ocular de la curva y, a continuación, empezar a pensar en una paramétrico generalizada modelo no lineal que captura bien su forma.

Answer 4

Para una visualización gráfica de la relación entre pa y eventos, se podría siguientes utilizando ggplot2 y R. Aquí hay un ejemplo usando sus propios datos de ejemplo:

library(ggplot2)
ggplot(data, aes(x=bloodpressure, y=event))+stat_smooth()

Las zonas de sombra indican que el 95% de intervalos de confianza.

Dos curvas pueden ser obtenidos con base en el uso de estatinas:

ggplot(data, aes(x=bloodpressure, y=event, group=factor(statin), color=factor(statin)))+stat_smooth()

Aquí los pacientes que toman estatinas son muy pocos, por lo tanto la curva no es completa.

Código anterior utiliza "loess" como el método. Puesto que usted está esperando a ser una curva, el código puede ser cambiado de la siguiente manera:

ggplot(data, aes(x=bloodpressure, y=event))+stat_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1)

También puede mostrar el suceso real de la tasa en los grupos de presión sanguínea. Dicen que dividir todo en 10 grupos de presión sanguínea. Para muestras de gran tamaño, el número de grupos puede ser aumentado para obtener más exacta de la curva. El eje y en las siguientes curvas de dar proporción de pacientes en los que el grupo que había eventos.

data$grp = cut(data$bloodpressure, 10)
aa = aggregate(event~grp, data=data, mean)
ggplot(aa, aes(x=grp, y=event))+geom_line(aes(group=1))

Errorbars se pueden añadir:

ddt = data.table(data)
aaa = ddt[,list(meanevent=mean(event), se_event=se(event)),by=grp]
ggplot(aaa, aes(x=grp, y=meanevent))+geom_line(aes(group=1)) + geom_errorbar(aes(ymin = meanevent-se_event, ymax = meanevent+se_event), width = 0.2)