Librerías que se van a utilizar.

library(kableExtra)
library(dplyr)
library(survival)
library(splines)
library(Boruta)
library(ggplot2)

Sea T una v.a continua. Demuestre que la vida media de T está dada por:

\[ E(T) = \int_{0}^{\infty} tf(t;\theta)dt = \int_{0}^{+\infty} [1 - F(t;\theta)] dt\] \[\text{Tenemos: } E(t) = \int_0^\infty tf(t;\theta)dt ~~~~\text{ recordedemos que: } f(t;\theta) = -S'(t;\theta)\] \[\text{Luego: } E(T) = \int_{0}^{\infty} t \cdot -S'(t;\theta)dt\] \[E(T) = - \int_{0}^{\infty} t \cdot S'(t;\theta)dt\] \[\text{Haciendo integración por partes: }\] \[E(T) = \underbrace{-tS(t; \theta) \displaystyle \left. \right|_{0}^{\infty}}_{\text{se cancela}} + \int_{0}^{\infty} S(t;\theta)dt\] \[E(T) = \int_{0}^{\infty} S(t;\theta)dt~~~~\text{ recordedemos que: } S(t;\theta) = 1 - F(t;\theta)\] \[\text{Entonces: } E(T) = \int_{0}^{\infty} 1 - F(t;\theta)dt\]

Hay varias formas de parametrizar una distribución de Weibull. La función survreg del R se basa en una familia de localización y escala, que es una parametrización diferente a la función rweibull (rweibull usa la fórmula vista en la clase 9), y a menudo lleva a confusión. En realidad, la relación es como sigue:

# survreg’s scale = 1/(rweibull shape)
# survreg’s intercept = log(rweibull scale)

Genere 200 valores de una Weibull con shape=7.2393 y scale=12.3559. Halle las estimaciones de máxima verosimilitud para \(\eta\) y \(\beta\) (ver clase 9) a partir de los datos generados; ¿qué se puede comentar de estas estimaciones respecto a las que se usaron para simular los valores? Con estos valores estimados para la Weibull, obtenga E(T) usando la formula del numeral anterior. Puede hacerlo numéricamente, por ejemplo:

#f<-function(x){exp(-(x/eta)^beta) }
#integrate(f, 0, Inf)

y compare su resultado con el promedio aritmético de los valores simulados.

Vamos a simular los datos de una Weibull con \(\eta = 12.3559\) y \(\beta = 7.2393\).

#Generar 200 observaciones de una weibull
wei <- rweibull(200, shape = 7.2393, scale = 12.3559)

Ahora se procederá a estimar los parámetros de escala y forma por medio de verosimilitud, haciendo uso de la función survreg.

survreg(Surv(wei) ~ 1, dist="weibull")
Call:
survreg(formula = Surv(wei) ~ 1, dist = "weibull")

Coefficients:
(Intercept) 
    2.53386 

Scale= 0.1355561 

Loglik(model)= -407.8   Loglik(intercept only)= -407.8
n= 200

Se observa que la estimación para el parámetro de escala es 0.1559346 y para el parámetro de forma es 2.514713, recordemos que survreg utiliza una parametrización diferente a rweibull, por lo que para poder comparar los valores estimados con los que se fijo para la simulación debemos realizar la transformación correspondiente:

# survreg’s scale = 1/(rweibull shape)
# survreg’s intercept = log(rweibull scale)

\[\text{rweibull shape } = \frac{1}{\text{survreg’s scale}} = \frac{1}{0.1559346} = 6.412945\]

\[\text{rweibull scale } = e^{\text{survreg’s intercept}} = e^{2.514713} = 12.36306\]

Recordemos que los valores que fijamos para la simulación son \(\eta = 12.3559\) y \(\beta = 7.2393\) y los estimados son \(\hat{\eta} = 12.36306\) y \(\hat{\beta} = 6.412945\), se observa que las diferencias entre ambos valores no son muy grandes, por otro lado, la estimación que se hace es sin censura, por lo que se espera que sean similares ambos valores, entonces podemos decir que survreg realiza una buena estimación de los parámetros de escala y forma.

Ahora se procederá a calcular la esperanza con la formula del literal a):

f <- function(x){exp(-(x/12.3559)^7.2393)}
integrate(f, 0, Inf)
11.57823 with absolute error < 3.7e-06

La esperanza de la weibull con parámetros \(\eta = 12.3559\) y \(\beta = 7.2393\) es 11.57823, comparemos con la media de los valores simulados.

mean(wei)
[1] 11.81423

Se observa que ambos valores son similares, tanto la esperanza teórica como la estimada, por lo tanto, la simulación que se realiza con rweibull es muy buena.

Sea T que tiene una distribución exponencial de dos parámetros la cual se denota \(T ~ exp(\theta, \gamma)\). Demuestre que:

la función definida por:

\[ f(t; \theta, \gamma) = \frac{1}{\theta} e^{- \frac{t-\gamma}{\theta}}, t \geq \gamma\]

es en efecto una densidad de probabilidad.

Para ser una densidad de probabilidad \(f(t; \theta, \gamma)\) debe cumplir las siguientes propiedades:

1) \[f(t; \theta, \gamma) \geq 0\]

Esto es cierto ya que toda función exponencial siempre toma valores entre \([0,+\infty)\) y con \(\theta > 0\) se puede evidenciar que en este caso particular \(f(t; \theta, \gamma)> 0\).

2) \[\int_{-\infty}^\infty f(t; \theta, \gamma) dt = 1\]

Lo anterior se puede escribir como:

\[\int_{\gamma}^\infty f(t; \theta, \gamma) dt= \int_{\gamma}^\infty \frac{1}{\theta} e^{- \frac{t-\gamma}{\theta}} dt = \frac{1}{\theta} (-e^{- \frac{t-\gamma}{\theta}})(\theta) \mid_{\gamma}^{h= \infty} = - e^{- \frac{t-\gamma}{\theta}} \mid_{\gamma}^{h= \infty} = 0 - (- e^{- \frac{h-\gamma}{\theta}}) \] Donde \(lim_{h \rightarrow 0} - e^{- \frac{h-\gamma}{\theta}} = - e^0 = -1\) entonces:

\[\int_{\gamma}^\infty f(t; \theta, \gamma) dt = 0 - (-1) = 1\] 3) Luego, como antes se mostró que la integral existe para todo el espaico parametral entonces la probabilidad de que t tome un valor en el intervalo \([a,b]\) es el área bajo la curva de la función de densidad en ese intervalo o lo que es lo mismo, la integral definida en dicho intervalo.

\[\int_{a}^b f(t; \theta, \gamma) dt = F(b) - F(a)\]

Luwgo, puesto que \(f(t; \theta, \gamma)\) satisface esas tres propiedades anteriores, se puede por tanto aseguar que \(f(t; \theta, \gamma)\) en efecto es una densidad de probabilidad.

y que efectivamente su vida media está dada por:

\[E(T) = \gamma + \theta\]

Se tiene que:

\[E(T)= \int_{\gamma}^{+\infty} tf(t;\theta, \gamma)dt = \int_{\gamma}^\infty \frac{t}{\theta} e^{- \frac{t-\gamma}{\theta}} dt\] Luego, sea \(u = \frac{t-\gamma}{\theta}\) entonces \(t = \theta u + \gamma\), de lo cual cuando \(t \longrightarrow \gamma\) por consiguiente \(u \longrightarrow 0\), además cuando \(t \longrightarrow \infty\) por consiguiente \(u \longrightarrow \infty\)

después, derivando a ambos lados se obtiene lo siguiente:

\[\frac{du}{dt} = \frac{1}{\theta}\] Luego la integral con respecto a u es como sigue:

\[\int_{0}^\infty (\theta u + \gamma) e^{-u} du = \theta \int_{0}^\infty u e^{-u} du + \gamma \int_{0}^\infty e^{-u} du \\ = \underbrace {\theta \int_{0}^\infty u e^{-u} du}_{\text{Gamma(2,1)}} + \gamma e^{-u} \mid_{0}^{\infty} \\ = \gamma * 1 + \underbrace {\theta \int_{0}^\infty \frac{u^{2-1}}{\Gamma(2)} e^{-u} du}_{\text{ Densidad Gamma(2,1)}} \\ = \gamma + \theta \] Por tanto, se concluye efectivamente que:

\[E(T) = \gamma + \theta\]

Considere los datos anexos acerca de pacientes con Cirrosis Biliar Primaria (PBC siglas en inglés) (PBC.txt). La descripción de los datos se encuentra en el archivo: Descripción datos PBC.txt).

Cirrosis Biliar Primaria

De un conjunto de datos casi idéntico que se encuentra en el apéndice D de Fleming y Harrington. Donde las diferencias con este conjunto de datos son: la edad que está en días, el status está codificado con 3 niveles y las variables sexo y status no faltan para las observaciones 313-418.

Los datos provienen del ensayo de Mayo Clinic sobre cirrosis biliar primaria (PBC, por sus siglas en inglés) del hígado realizado entre 1974 y 1984. Un total de 424 pacientes con PBC, derivados a Mayo Clinic durante ese intervalo de diez años, cumplieron con los criterios de elegibilidad para el estudio aleatorizado controlado con placebo. ensayo del fármaco D-penicilamina. Los primeros 312 casos en el conjunto de datos participaron en el ensayo aleatorio y contienen datos en gran parte completos. Los 112 casos adicionales no participaron en el ensayo clínico, pero aceptaron que se registraran las medidas básicas y que se les hiciera un seguimiento para sobrevivir. Seis de esos casos se perdieron durante el seguimiento poco después del diagnóstico, por lo que los datos aquí son de 106 casos adicionales, así como de los 312 participantes asignados al azar. Los elementos de datos que faltan se indican con un punto.

Preprocesamiento de datos:

Primero se procede a realizar la lectura de la base de datos.

# lectura de la base de datos
PBC <- read.table("PBC.txt", quote="\"", comment.char="")

Para facilitar el uso de la base de datos se renombra las variables.

# renombrando variables
colnames(PBC) <- c("ID","dias", "status", "medicamento","edad", "sexo","ascitis", "hepatomegalia","varices","edema", "bilirrubina","colesterol","albumina", "cobre_orina","fosfatasa", "SGOT","trigliceridos", "plaquetas","protombina","histologico")
head(PBC)

Descripción de las variables

Las variables que contemplan en la base de datos PBC son:

• ID: Número único de identificación del paciente.
• dias: número de días entre el registro y la fecha anterior de muerte, trasplante o análisis del estudio en julio de 1986.
• status: Tres niveles que son: 0 = vivo, 1 = trasplante de hígado, 2 = muerto.
• medicamento: Dos niveles 1= D-penicilamina, 2=placebo.
• edad: Medida en días.
• sexo: 0= Masculino, 1= Femenino.
• ascitis: Presencia de ascitis: 0=no 1=sí.
• hepatomegalia: Presencia de hepatomegalia 0=no 1=si.
• varices: Presencia de varices 0=no 1=si.
• edema: Presencia de edema 0=sin edema y sin terapia diurética para el edema; 0,5 = edema presente sin diuréticos o edema resuelto con diuréticos; 1 = edema a pesar de la terapia con diuréticos.
• bilirrubina: Bilirrubina sérica en mg/dl.
• colesterol: Colesterol sérico en mg/dl.
• albumina: Albúmina en gm/dl.
• cobre_orina: Cobre en orina en ug/día (ug = un microgramo= 1 millonésimo gramo).
• fosfatasa: Fosfatasa alcalina en U/litro (U = ¼ de pulgada).
• SGOT: Valores de transaminasa en la sangre en U/ml (U = ¼ de pulgada).
• triglicéridos: Triglicéridos en mg/dl.
• plaquetas: Plaquetas por ml cúbico / 1000.
• protombina: Tiempo de protrombina en segundos.
• histologico: Etapa histológica de la enfermedad.

Se adecuará la variable status, para que solo quede con dos niveles, 1 = censura, 0 = vivo.

PBC$status <- ifelse(PBC$status %in% c(0, 1), 0, no = 1)

Como se menciono anteriormente los datos faltantes se indican con un punto, entonces para efectos prácticos se transformaran a un objeto del tipo NA.

PBC[PBC == "."] <- NA

Se organizan las variables para que las que son categóricas sean de tipo factor en R.

#Convertir las variables necesarias
#PBC$status <- as.factor(PBC$status)
PBC$medicamento <- as.factor(PBC$medicamento)
PBC$sexo <- as.factor(PBC$sexo)
PBC$ascitis <- as.factor(PBC$ascitis)
PBC$hepatomegalia <- as.factor(PBC$hepatomegalia)
PBC$varices <- as.factor(PBC$varices)
PBC$edema <- as.factor(PBC$edema)
PBC$colesterol <- as.integer(PBC$colesterol)
PBC$cobre_orina <- as.integer(PBC$cobre_orina)
PBC$fosfatasa <- as.numeric(PBC$fosfatasa)
PBC$SGOT <- as.numeric(PBC$SGOT)
PBC$trigliceridos <- as.integer(PBC$trigliceridos)
PBC$plaquetas <- as.integer(PBC$plaquetas)
PBC$protombina <- as.numeric(PBC$protombina)
PBC$histologico <- as.factor(PBC$histologico)

Se realizará una tabla para contar los NA.

#Contar los NA
nas <- apply(X = is.na(PBC), MARGIN = 2, FUN = sum)
kable(nas)
x
ID 0
dias 0
status 0
medicamento 106
edad 0
sexo 0
ascitis 106
hepatomegalia 106
varices 106
edema 0
bilirrubina 0
colesterol 134
albumina 0
cobre_orina 108
fosfatasa 106
SGOT 106
trigliceridos 136
plaquetas 11
protombina 2
histologico 6

De la tabla se observa que hay una buena cantidad de datos faltantes, pero se sabe que hay 112 casos adicionales que no participaron en el ensayo clínico, pero aceptaron que se registraran las medidas básicas, lo que causo una gran cantidad de datos faltantes en otras variables, sin embargo, se construirán modelos con los datos faltantes ya que no es adecuado suprimir información epidemiológica.

Graficos descriptivos

Se realiza boxplots para todas las variables categóricas:

par(pty = 's')
sexo <-PBC$sexo
dias = PBC$dias

#sexo
p1 <- ggplot(PBC, aes(x=sexo,y=dias)) + geom_boxplot(aes(fill=(sexo))) + theme_bw()

#ascitis
ascitis<- PBC$ascitis
p2 <- ggplot(PBC, aes(x=ascitis,y=dias)) + geom_boxplot(aes(fill=(ascitis))) + theme_bw()
# hepatomegalia
hepatomegalia<- PBC$hepatomegalia
p3 <- ggplot(PBC, aes(x=hepatomegalia,y=dias)) + geom_boxplot(aes(fill=(hepatomegalia))) + theme_bw()

varices<- PBC$varices
p4 <- ggplot(PBC, aes(x=varices,y=dias)) + geom_boxplot(aes(fill=(varices))) + theme_bw()
#edema
edema<- PBC$edema
p5 <- ggplot(PBC, aes(x=edema,y=dias)) + geom_boxplot(aes(fill=(edema))) + theme_bw()
# histologico
histologico<- PBC$histologico
p6 <- ggplot(PBC, aes(x=histologico,y=dias)) + geom_boxplot(aes(fill=(histologico))) + theme_bw()

#medicamento
medicamento<- PBC$medicamento
p7 <- ggplot(PBC, aes(x=medicamento,y=dias)) + geom_boxplot(aes(fill=(medicamento))) + theme_bw()


library(gridExtra)
grid.arrange(p1, p2, p3, p4, p5, p6, p7, ncol=2, nrow = 4)

De manera descriptiva, se observa que las variables edema, ascitis y histologico afectan en el comportamiento de la variable días.

A manera descriptiva se realizará una matriz de correlación sin los datos faltantes.

library(corrplot)
base <- na.omit(PBC)
continuas <- select_if(base, is.numeric)[-1]

# create correlation matrix
corr.data <- cor(continuas)
# correlation plot of quantitative variables
corrplot(corr.data, method = 'ellipse', order='AOE', type = 'upper')

Se observa que las variables que tienen mayor correlación con la variable dias es cobre_orina y bilirubina.

Usando estos datos seleccione y ajuste un modelo Weibull.

Resumen modelo completo

attach(PBC)

#Resumen de la varible de censura segun el tiempo de falla

res <- by(PBC$dias,PBC$status,summary)

# MODELO WEIBULL CON TODAS LAS VARIABLES (17)

resultW <- survreg(Surv(dias,status)~ medicamento+edad+sexo+ascitis+hepatomegalia+varices+
                  edema+bilirrubina+colesterol+albumina+cobre_orina+fosfatasa+SGOT+
                  trigliceridos+plaquetas+protombina+histologico,  
                dist="weibull",data=PBC)

# Resumen del modelo
summary(resultW)

Call:
survreg(formula = Surv(dias, status) ~ medicamento + edad + sexo + 
    ascitis + hepatomegalia + varices + edema + bilirrubina + 
    colesterol + albumina + cobre_orina + fosfatasa + SGOT + 
    trigliceridos + plaquetas + protombina + histologico, data = PBC, 
    dist = "weibull")
                   Value Std. Error     z       p
(Intercept)     1.15e+01   1.47e+00  7.80 6.4e-15
medicamento2    1.04e-01   1.32e-01  0.79  0.4303
edad           -4.92e-05   1.92e-05 -2.56  0.0106
sexo1           2.13e-01   1.90e-01  1.12  0.2628
ascitis1       -4.40e-02   2.34e-01 -0.19  0.8507
hepatomegalia1 -3.54e-02   1.53e-01 -0.23  0.8175
varices1       -6.22e-03   1.51e-01 -0.04  0.9672
edema0.5       -1.58e-01   2.01e-01 -0.78  0.4326
edema1         -7.85e-01   2.44e-01 -3.22  0.0013
bilirrubina    -4.59e-02   1.50e-02 -3.05  0.0023
colesterol     -2.79e-04   2.73e-04 -1.02  0.3077
albumina        3.97e-01   1.82e-01  2.18  0.0291
cobre_orina    -1.51e-03   7.09e-04 -2.13  0.0328
fosfatasa      -5.12e-06   2.45e-05 -0.21  0.8346
SGOT           -2.40e-03   1.20e-03 -2.00  0.0452
trigliceridos   1.83e-04   8.62e-04  0.21  0.8321
plaquetas      -4.70e-04   7.15e-04 -0.66  0.5113
protombina     -1.81e-01   7.21e-02 -2.51  0.0122
histologico2   -8.87e-01   6.61e-01 -1.34  0.1796
histologico3   -1.03e+00   6.44e-01 -1.60  0.1088
histologico4   -1.31e+00   6.53e-01 -2.00  0.0455
Log(scale)     -5.04e-01   7.69e-02 -6.55 5.8e-11

Scale= 0.604 

Weibull distribution
Loglik(model)= -966.1   Loglik(intercept only)= -1053.2
    Chisq= 174.24 on 20 degrees of freedom, p= 1.3e-26 
Number of Newton-Raphson Iterations: 7 
n=276 (142 observations deleted due to missingness)

Call:
survreg(formula = Surv(dias, status) ~ edad + edema + bilirrubina + 
    albumina + cobre_orina + SGOT + protombina + histologico, 
    data = base, dist = "weibull")
                 Value Std. Error     z       p
(Intercept)   1.14e+01   1.36e+00  8.34 < 2e-16
edad         -5.24e-05   1.71e-05 -3.07  0.0022
edema0.5     -1.16e-01   1.84e-01 -0.63  0.5297
edema1       -6.21e-01   2.07e-01 -3.00  0.0027
bilirrubina  -5.10e-02   1.12e-02 -4.55 5.3e-06
albumina      4.18e-01   1.64e-01  2.54  0.0109
cobre_orina  -1.80e-03   5.93e-04 -3.04  0.0024
SGOT         -2.50e-03   1.11e-03 -2.26  0.0241
protombina   -1.72e-01   6.93e-02 -2.48  0.0132
histologico2 -8.59e-01   6.64e-01 -1.29  0.1959
histologico3 -1.03e+00   6.46e-01 -1.60  0.1105
histologico4 -1.27e+00   6.43e-01 -1.97  0.0491
Log(scale)   -5.01e-01   7.56e-02 -6.63 3.3e-11

Scale= 0.606 

Weibull distribution
Loglik(model)= -968   Loglik(intercept only)= -1053.2
    Chisq= 170.45 on 11 degrees of freedom, p= 9.5e-31 
Number of Newton-Raphson Iterations: 7 
n= 276

resumen modelo reducido

# Ajustando el modelo weibull con esas variables:

modWeibull <- survreg(Surv(dias, status) ~ edad + edema + bilirrubina + albumina + 
    cobre_orina + SGOT + protombina + histologico, dist="weibull",data=PBC)

summary(modWeibull)

Call:
survreg(formula = Surv(dias, status) ~ edad + edema + bilirrubina + 
    albumina + cobre_orina + SGOT + protombina + histologico, 
    data = PBC, dist = "weibull")
                 Value Std. Error     z       p
(Intercept)   1.17e+01   1.31e+00  8.98 < 2e-16
edad         -5.37e-05   1.59e-05 -3.39 0.00070
edema0.5     -9.60e-02   1.75e-01 -0.55 0.58319
edema1       -6.11e-01   1.86e-01 -3.30 0.00098
bilirrubina  -5.07e-02   1.08e-02 -4.68 2.9e-06
albumina      4.63e-01   1.54e-01  3.01 0.00260
cobre_orina  -1.84e-03   5.58e-04 -3.29 0.00099
SGOT         -2.60e-03   1.03e-03 -2.53 0.01135
protombina   -2.04e-01   6.47e-02 -3.16 0.00160
histologico2 -1.04e+00   6.63e-01 -1.57 0.11689
histologico3 -1.17e+00   6.47e-01 -1.80 0.07123
histologico4 -1.39e+00   6.44e-01 -2.15 0.03143
Log(scale)   -4.97e-01   7.13e-02 -6.98 2.9e-12

Scale= 0.608 

Weibull distribution
Loglik(model)= -1078.3   Loglik(intercept only)= -1179.7
    Chisq= 202.83 on 11 degrees of freedom, p= 1.9e-37 
Number of Newton-Raphson Iterations: 7 
n=310 (108 observations deleted due to missingness)

Usando estos datos seleccione y ajuste un modelo de Cox.

Resumen modelo completo

#Modelo completo 

resultadoC <- coxph(Surv(dias,status)~ medicamento+edad+sexo+ascitis+hepatomegalia+
                      varices+edema+bilirrubina+colesterol+albumina+cobre_orina
                    +fosfatasa+SGOT+trigliceridos+plaquetas+protombina+histologico,  
                    data=PBC)
summary(resultadoC)
Call:
coxph(formula = Surv(dias, status) ~ medicamento + edad + sexo + 
    ascitis + hepatomegalia + varices + edema + bilirrubina + 
    colesterol + albumina + cobre_orina + fosfatasa + SGOT + 
    trigliceridos + plaquetas + protombina + histologico, data = PBC)

  n= 276, number of events= 111 
   (142 observations deleted due to missingness)

                     coef  exp(coef)   se(coef)      z Pr(>|z|)   
medicamento2   -1.765e-01  8.382e-01  2.189e-01 -0.806  0.42003   
edad            7.910e-05  1.000e+00  3.186e-05  2.483  0.01304 * 
sexo1          -3.757e-01  6.868e-01  3.111e-01 -1.208  0.22720   
ascitis1        5.899e-04  1.001e+00  3.951e-01  0.001  0.99881   
hepatomegalia1  5.648e-02  1.058e+00  2.565e-01  0.220  0.82575   
varices1        7.011e-02  1.073e+00  2.491e-01  0.281  0.77839   
edema0.5        2.436e-01  1.276e+00  3.354e-01  0.726  0.46764   
edema1          1.146e+00  3.146e+00  4.166e-01  2.751  0.00595 **
bilirrubina     8.014e-02  1.083e+00  2.628e-02  3.050  0.00229 **
colesterol      4.731e-04  1.000e+00  4.534e-04  1.043  0.29675   
albumina       -7.496e-01  4.725e-01  3.104e-01 -2.415  0.01572 * 
cobre_orina     2.442e-03  1.002e+00  1.171e-03  2.087  0.03693 * 
fosfatasa       5.228e-07  1.000e+00  4.066e-05  0.013  0.98974   
SGOT            3.710e-03  1.004e+00  1.992e-03  1.863  0.06247 . 
trigliceridos  -5.177e-04  9.995e-01  1.436e-03 -0.361  0.71846   
plaquetas       8.410e-04  1.001e+00  1.189e-03  0.707  0.47929   
protombina      2.763e-01  1.318e+00  1.167e-01  2.367  0.01793 * 
histologico2    1.406e+00  4.081e+00  1.080e+00  1.302  0.19288   
histologico3    1.683e+00  5.380e+00  1.052e+00  1.600  0.10968   
histologico4    2.119e+00  8.325e+00  1.067e+00  1.985  0.04711 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

               exp(coef) exp(-coef) lower .95 upper .95
medicamento2      0.8382     1.1930    0.5458    1.2872
edad              1.0001     0.9999    1.0000    1.0001
sexo1             0.6868     1.4560    0.3732    1.2637
ascitis1          1.0006     0.9994    0.4613    2.1704
hepatomegalia1    1.0581     0.9451    0.6400    1.7494
varices1          1.0726     0.9323    0.6582    1.7479
edema0.5          1.2759     0.7838    0.6611    2.4623
edema1            3.1456     0.3179    1.3901    7.1178
bilirrubina       1.0834     0.9230    1.0291    1.1407
colesterol        1.0005     0.9995    0.9996    1.0014
albumina          0.4725     2.1162    0.2572    0.8682
cobre_orina       1.0024     0.9976    1.0001    1.0047
fosfatasa         1.0000     1.0000    0.9999    1.0001
SGOT              1.0037     0.9963    0.9998    1.0076
trigliceridos     0.9995     1.0005    0.9967    1.0023
plaquetas         1.0008     0.9992    0.9985    1.0032
protombina        1.3183     0.7586    1.0487    1.6572
histologico2      4.0811     0.2450    0.4914   33.8956
histologico3      5.3797     0.1859    0.6845   42.2793
histologico4      8.3253     0.1201    1.0274   67.4610

Concordance= 0.848  (se = 0.018 )
Likelihood ratio test= 169.7  on 20 df,   p=<2e-16
Wald test            = 175  on 20 df,   p=<2e-16
Score (logrank) test = 295.6  on 20 df,   p=<2e-16
Call:
coxph(formula = Surv(dias, status) ~ edad + edema + bilirrubina + 
    albumina + cobre_orina + SGOT + protombina + histologico, 
    data = base)

  n= 276, number of events= 111 

                   coef  exp(coef)   se(coef)      z Pr(>|z|)    
edad          8.575e-05  1.000e+00  2.817e-05  3.043  0.00234 ** 
edema0.5      1.598e-01  1.173e+00  3.055e-01  0.523  0.60090    
edema1        9.122e-01  2.490e+00  3.545e-01  2.573  0.01009 *  
bilirrubina   8.693e-02  1.091e+00  1.981e-02  4.388 1.14e-05 ***
albumina     -7.387e-01  4.777e-01  2.792e-01 -2.645  0.00816 ** 
cobre_orina   2.787e-03  1.003e+00  9.912e-04  2.811  0.00493 ** 
SGOT          3.956e-03  1.004e+00  1.841e-03  2.148  0.03168 *  
protombina    2.642e-01  1.302e+00  1.123e-01  2.353  0.01862 *  
histologico2  1.360e+00  3.895e+00  1.081e+00  1.258  0.20843    
histologico3  1.682e+00  5.378e+00  1.048e+00  1.605  0.10839    
histologico4  2.063e+00  7.867e+00  1.043e+00  1.977  0.04801 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

             exp(coef) exp(-coef) lower .95 upper .95
edad            1.0001     0.9999    1.0000    1.0001
edema0.5        1.1733     0.8523    0.6447    2.1352
edema1          2.4897     0.4017    1.2427    4.9881
bilirrubina     1.0908     0.9167    1.0493    1.1340
albumina        0.4777     2.0932    0.2764    0.8258
cobre_orina     1.0028     0.9972    1.0008    1.0047
SGOT            1.0040     0.9961    1.0003    1.0076
protombina      1.3024     0.7678    1.0451    1.6230
histologico2    3.8947     0.2568    0.4682   32.3981
histologico3    5.3782     0.1859    0.6897   41.9364
histologico4    7.8672     0.1271    1.0182   60.7865

Concordance= 0.845  (se = 0.019 )
Likelihood ratio test= 165.7  on 11 df,   p=<2e-16
Wald test            = 176.9  on 11 df,   p=<2e-16
Score (logrank) test = 278.8  on 11 df,   p=<2e-16

Resumen del modelo reducido

# Ajustando el modelo COX con esas variables:

modCox <- coxph(Surv(dias,status)~ edad + edema + bilirrubina + albumina + 
    cobre_orina + SGOT + protombina + histologico,  data=PBC)

summary(modCox)
Call:
coxph(formula = Surv(dias, status) ~ edad + edema + bilirrubina + 
    albumina + cobre_orina + SGOT + protombina + histologico, 
    data = PBC)

  n= 310, number of events= 124 
   (108 observations deleted due to missingness)

                   coef  exp(coef)   se(coef)      z Pr(>|z|)    
edad          8.964e-05  1.000e+00  2.610e-05  3.434 0.000595 ***
edema0.5      1.446e-01  1.156e+00  2.879e-01  0.502 0.615454    
edema1        9.208e-01  2.511e+00  3.165e-01  2.909 0.003626 ** 
bilirrubina   8.641e-02  1.090e+00  1.909e-02  4.527 5.98e-06 ***
albumina     -7.970e-01  4.507e-01  2.599e-01 -3.066 0.002167 ** 
cobre_orina   2.721e-03  1.003e+00  9.304e-04  2.924 0.003454 ** 
SGOT          4.229e-03  1.004e+00  1.681e-03  2.516 0.011860 *  
protombina    3.158e-01  1.371e+00  1.038e-01  3.044 0.002336 ** 
histologico2  1.655e+00  5.232e+00  1.076e+00  1.538 0.124085    
histologico3  1.899e+00  6.681e+00  1.047e+00  1.813 0.069779 .  
histologico4  2.243e+00  9.423e+00  1.041e+00  2.154 0.031252 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

             exp(coef) exp(-coef) lower .95 upper .95
edad            1.0001     0.9999    1.0000    1.0001
edema0.5        1.1556     0.8653    0.6572    2.0320
edema1          2.5114     0.3982    1.3504    4.6704
bilirrubina     1.0903     0.9172    1.0502    1.1318
albumina        0.4507     2.2189    0.2708    0.7501
cobre_orina     1.0027     0.9973    1.0009    1.0046
SGOT            1.0042     0.9958    1.0009    1.0076
protombina      1.3714     0.7292    1.1190    1.6807
histologico2    5.2322     0.1911    0.6349   43.1162
histologico3    6.6806     0.1497    0.8576   52.0382
histologico4    9.4235     0.1061    1.2238   72.5650

Concordance= 0.851  (se = 0.018 )
Likelihood ratio test= 197.2  on 11 df,   p=<2e-16
Wald test            = 205.7  on 11 df,   p=<2e-16
Score (logrank) test = 328.3  on 11 df,   p=<2e-16

Compare ambos modelos en términos de los parámetros estimados comunes.

Una pregunta de interés es si la supervivencia en los grupos sex, presence of ascites y presence of spiders está relacionada con las covariables. ¿Qué se puede decir en este caso? Grafique las respectivas curvas de supervivencia para el modelo Weibull y el de Cox. Comente.