En aplicaciones biomédicas las técnicas no paramétricas (p. ej. estimador KM) y semiparamétricas (p. ej. modelo de riesgos proporcionales de Cox) suelen ser las más empleadas por su flexibilidad para ajustarse a formas arbitrarias de la función de riesgo. Sin embargo, los modelos paramétricos siguen siendo útiles cuando los datos de supervivencia se ajustan aproximadamente a una familia conocida, ya que:
Anteriormente se introdujeron varias familias paramétricas (exponencial, Weibull, gamma, …). En este capítulo se hace hincapié en la distribución exponencial y la Weibull, las más utilizadas en práctica, y se comentan brevemente otros modelos.
La distribución exponencial es el caso más simple y se caracteriza por tener una tasa de riesgo constante:
\[ h(t)=\lambda, \]
lo que le otorga la propiedad de pérdidad de memoria: el riesgo instantáneo es igual en cualquier punto del tiempo al riesgo inicial.
Para construir la verosimilitud con datos censurados a la derecha, basta multiplicar
Por su sencillez, la exponencial es muy útil en cálculos de potencia y tamaño de muestra. No obstante, al carecer de flexibilidad en la forma de \(h(t)\) (siempre constante), a menudo se prefiere la Weibull —que añade un parámetro de forma— para modelar supervivencias reales.
Nota: cuando \(T\sim\mathrm{Exp}(\lambda)\), entonces \(T^\alpha\sim\mathrm{Weibull}(\alpha,\lambda)\), lo que muestra que la exponencial es un caso particular de la Weibull.
La distribución Weibull, con parámetros de forma \(\alpha\) y escala \(\lambda\), tiene
\[
h(t)=\alpha\,\lambda^{\alpha}\,t^{\alpha-1},
\qquad
S(t)=\exp\bigl[-(\lambda t)^{\alpha}\bigr].
\]
Usando en su lugar
\[
\sigma=\frac1\alpha,
\quad
\mu=-\log\lambda,
\]
se reescriben
\[
h(t)=\frac1\sigma\,e^{-\mu/\sigma}\,t^{1/\sigma-1},
\qquad
S(t)=\exp\bigl[-e^{-\mu/\sigma}\,t^{1/\sigma}\bigr].
\]
Si aplicamos la transformación complemento‐log‐log
\[
g(u)=\log\bigl[-\log(u)\bigr],
\]
a \(S_i=S(t_i)\) obtenemos la relación lineal
\[
\log\bigl[-\log S_i\bigr]
=\alpha\log\lambda+\alpha\log t_i
=-\frac\mu\sigma+\frac1\sigma\log t_i.
\]
Por tanto, para valorar si los datos siguen una Weibull, se procede así:
y_i <- log(-log( survfit...$surv ))
x_i <- log( survfit...$time )Ajustar la regresión lineal
fit <- lm( y_i ~ x_i )cuyo coeficiente de pendiente \(m\) estima \(1/\sigma\) y la ordenada en el origen \(b\) estima \(-\mu/\sigma\).
Ejemplo con gastricXelox:
library(survival)
library(asaur)
timeM <- gastricXelox$timeWeeks * 7/30.25
delta <- gastricXelox$delta
library(survival)
km <- survfit(Surv(timeM, delta) ~ 1)
survEst <- km$surv
survTime <- km$time
logLogSurv <- log(-log(survEst))
logTime <- log(survTime)
plot(logLogSurv ~ logTime)
lm1 <- lm(logLogSurv ~ logTime)
abline(lm1)
El gráfico anterior muestra desviaciones de la recta, sugiriendo que Weibull no es adecuado.
Ejemplo con pharmacoSmoking:
attach(pharmacoSmoking)
ttr[ttr == 0] <- 0.5 # corregir ceros
km2 <- survfit(Surv(ttr, relapse) ~ 1)
survEst2 <- km2$surv
survTime2 <- km2$time
logLogSE <- log(-log(survEst2))
logST <- log(survTime2)
fit2 <- lm(logLogSE ~ logST)
coef(fit2)(Intercept) logST
-2.003177 0.438493 # (Intercept) logST
# -2.0032 0.4385
plot(logLogSE ~ logST)
abline(fit2)
Sea una muestra de tamaño \(n\) con tiempos \(t_i\) y indicadores \(\delta_i\) (1=evento, 0=censura). La función de verosimilitud para el modelo Weibull, con tasa de eventos \(d=\sum\delta_i\), es
\[ l(\lambda,\alpha) =\sum_{i=1}^n \Bigl[\delta_i\log h(t_i)+\log S(t_i)\Bigr] =d\log\alpha+d\,\alpha\log\lambda +(\alpha-1)\sum_{i=1}^n\delta_i\log t_i -\lambda^\alpha\sum_{i=1}^n t_i^\alpha, \]
donde \(h(t)=\alpha\lambda^\alpha t^{\alpha-1}\) y \(S(t)=\exp[-(\lambda t)^\alpha]\). Esto debido a que:
\[L(\lambda,\alpha)=\prod_{i=1}^nf(t_i;\lambda,\alpha)^{\delta_i}S(t_i)^{1-\delta_i}=\prod_{i=1}^nh(t_i;\lambda,\alpha)^{\delta_i}S(t_i)\]
Para facilitar la optimización, reparametrizamos en términos de
\[
\sigma=\frac1\alpha,
\quad
\mu=-\log\lambda,
\]
y codificamos la log‐verosimilitud en R como:
logLikWeib <- function(par, tt, status) {
mu <- par[1]
sigma <- par[2]
lambda <- exp(-mu)
alpha <- 1/sigma
d <- sum(status)
sum1 <- sum(status * log(tt))
sum2 <- sum(tt^alpha)
term1 <- d*log(alpha) + alpha*d*log(lambda)
term2 <- (alpha - 1)*sum1
term3 <- (lambda^alpha)*sum2
term1 + term2 - term3
}Usamos optim para maximizarla (fnscale=-1) partiendo de \((\hat\mu,\hat\sigma)\) de la regresión lineal anterior:
# tt = vector de tiempos; status = vector de 0/1
result <- optim(par = c(2.5, 2.5),
fn = logLikWeib,
method = "L-BFGS-B",
lower = c(0.001, 0.01),
upper = c(5, 5),
control=list(fnscale = -1),
tt = ttr, status = relapse)
result$par[1] 4.656329 2.041061# mu.hat = 4.6563, sigma.hat = 2.0411Una forma más sencilla es usar survreg con distribución "weibull":
library(survival)
result.sr <- survreg(Surv(ttr, relapse) ~ 1, dist="weibull")
summary(result.sr)
Call:
survreg(formula = Surv(ttr, relapse) ~ 1, dist = "weibull")
Value Std. Error z p
(Intercept) 4.6563 0.2170 21.46 < 2e-16
Log(scale) 0.7135 0.0919 7.76 8.3e-15
Scale= 2.04
Weibull distribution
Loglik(model)= -466.1 Loglik(intercept only)= -466.1
Number of Newton-Raphson Iterations: 5
n= 125 Se aprovecha la propiedad de que, para un valor fijo de \(\alpha\): \[T^*=T^\alpha\sim\mathrm{Exp}(\lambda^\alpha)\] de modo que el MLE de \(\lambda\) es
\[ \hat\lambda(\alpha) =\Bigl(\frac{d}{\sum_i t_i^\alpha}\Bigr)^{1/\alpha}, \quad d=\sum_i\delta_i. \]
Sustituyendo \(\lambda=\hat\lambda(\alpha)\) en la log‐verosimilitud completa obtenemos la verosimilitud perfilada:
\[ l^*(\alpha)=l\bigl(\hat\lambda(\alpha),\,\alpha\bigr), \]
que solo depende de \(\alpha\). En la práctica, trabajamos con \(\sigma=1/\alpha\) y reescribimos la función en R así:
logLikWeibProf <- function(par, tt, status) {
# par = sigma
sigma <- par
alpha <- 1/sigma
d <- sum(status) # número de eventos
sumt <- sum(status * log(tt))
sumta <- sum(tt^alpha)
# MLE de lambda dado alpha:
lambda <- (d/sumta)^(1/alpha)
# términos de la log‐verosimilitud:
term1 <- d*log(alpha) + alpha*d*log(lambda)
term2 <- (alpha - 1)*sumt
term3 <- (lambda^alpha)*sumta
term1 + term2 - term3
}Luego maximizamos sobre \(\sigma\) con optim (fnscale = –1 para buscar el máximo):
resultProf <- optim(par=2.280,
fn=logLikWeibProf,
method="L-BFGS-B",
lower=0.01, upper=5,
control=list(fnscale=-1),
tt=ttr, status=relapse)
sigma.hat <- resultProf$par
# sigma.hat ≈ 2.041063Obtenemos entonces
alpha.hat <- 1/sigma.hat
sumta <- sum(ttr^alpha.hat)
lambda.hat <- (sum(relapse)/sumta)^(1/alpha.hat)
mu.hat <- -log(lambda.hat)Finalmente, para visualizar el perfil de \(l^\*(\sigma)\):
sigma.list <- seq(1, 5, by=0.01)
logLik.list <- sapply(sigma.list,
logLikWeibProf,
tt=ttr, status=relapse)
plot(sigma.list, logLik.list, type="l",
xlab="σ", ylab="Log‐verosimilitud perfilada")
abline(v=sigma.hat, col="gray")
Este gráfico muestra el máximo en \(\hat\sigma\), confirmando los valores obtenidos por optimización.
Podemos ajustar una Weibull que pase exactamente por dos puntos \((t_{1},s_{1})\) y \((t_{2},s_{2})\) del estimador de Kaplan–Meier resolviendo el sistema lineal que surge de la transformación \[y_i=\log\bigl[-\log(s_i)\bigr]=\alpha\log\lambda+\alpha\log(t_i),\quad i=1,2.\] De ahí se obtiene la solución cerrada
\[ \tilde\alpha=\frac{y_1-y_2}{\log(t_1)-\log(t_2)}, \qquad \tilde\lambda =\exp\!\Bigl(\frac{y_2\,\log(t_1)-y_1\,\log(t_2)}{y_1-y_2}\Bigr). \]
Ejemplo (grupo “patchOnly” en pharmacoSmoking, días 28 y 84):
# 1) Estimador KM en t=28,84:
result.surv <- survfit(Surv(ttr, relapse) ~ 1,
subset=grp=="patchOnly")
result.summ <- summary(result.surv, time=c(28,84))
t.vec <- result.summ$time # 28, 84
s.vec <- result.summ$surv # 0.4375, 0.265625
data.frame(t.vec, s.vec) t.vec s.vec
1 28 0.437500
2 84 0.265625# 2) Ajuste Weibull "ad hoc" con Hmisc::Weibull2
library(Hmisc)
Attaching package: 'Hmisc'The following objects are masked from 'package:base':
format.pval, unitspharmWeib <- Weibull2(t.vec, s.vec)
t.vals <- 1:200
s.vals <- pharmWeib(t.vals)# 3) Ajuste MLE completo vía survreg
model.pharm.weib.basic <- survreg(Surv(ttr, relapse) ~ 1,
dist="weibull",
subset=grp=="patchOnly")
mu.hat <- model.pharm.weib.basic$coef # intercept = -log λ
sigma.hat <- model.pharm.weib.basic$scale # = 1/α
lambda.hat <- exp(-mu.hat)
alpha.hat <- 1/sigma.hat
# Función de supervivencia Weibull MLE
s.mle.vals <- 1 - pweibull(t.vals,
shape=alpha.hat,
scale=1/lambda.hat)# 4) Comparación gráfica
plot(result.surv, conf.int=FALSE,
xlab="Días hasta recaída",
ylab="Prob. supervivencia")
lines(s.mle.vals ~ t.vals, col="blue") # Weibull MLE
lines(s.vals ~ t.vals, col="red") # Weibull por dos puntos
points(t.vec, s.vec, col="red") # puntos de ajuste
En la figura anterior se ve que la curva roja (Weibull por dos puntos) y la azul (MLE) son bastante similares, y ambas aproximan bien al escalón de Kaplan–Meier.
En el modelo Accelerated Failure Time (AFT) se asume que el tiempo de supervivencia bajo tratamiento, \(T_{1}\), es un múltiplo del que habría tenido bajo control, \(T_{0}\):
\[ T_{1}=e^{\gamma}T_{0},\quad S_{1}(t)=S_{0}(e^{-\gamma}t). \]
Para una distribución Weibull con parámetros \((\mu,\sigma)\) (donde recordemos \(\lambda=e^{-\mu}\) y \(\alpha=1/\sigma\)), se tiene
\[ h_{0}(t)=\frac{1}{\sigma}e^{-\mu/\sigma}\,t^{1/\sigma-1},\quad h_{1}(t)=e^{-\gamma}h_{0}(e^{-\gamma}t) =e^{-\gamma/\sigma}h_{0}(t). \]
Así, en Weibull el modelo AFT es equivalente al modelo de hazard proporcional con
\[ e^{\beta}=e^{-\gamma/\sigma} \quad\Longrightarrow\quad \beta=-\frac{\gamma}{\sigma}. \]
Ajuste con survreg (Weibull AFT)
# Ajuste AFT Weibull
result.survreg.grp <- survreg(
Surv(ttr, relapse) ~ grp,
dist="weibull",
)
summary(result.survreg.grp)
Call:
survreg(formula = Surv(ttr, relapse) ~ grp, dist = "weibull")
Value Std. Error z p
(Intercept) 5.2859 0.3320 15.92 <2e-16
grppatchOnly -1.2514 0.4348 -2.88 0.004
Log(scale) 0.6888 0.0911 7.56 4e-14
Scale= 1.99
Weibull distribution
Loglik(model)= -461.8 Loglik(intercept only)= -466.1
Chisq= 8.63 on 1 degrees of freedom, p= 0.0033
Number of Newton-Raphson Iterations: 5
n= 125 Comparación con Cox (PH semiparamétrico)
# Ajuste Cox
result.coxph.grp <- coxph(Surv(ttr, relapse) ~ grp,
data=pharmacoSmoking)
summary(result.coxph.grp)Call:
coxph(formula = Surv(ttr, relapse) ~ grp, data = pharmacoSmoking)
n= 125, number of events= 89
coef exp(coef) se(coef) z Pr(>|z|)
grppatchOnly 0.6050 1.8313 0.2161 2.8 0.00511 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
grppatchOnly 1.831 0.5461 1.199 2.797
Concordance= 0.581 (se = 0.027 )
Likelihood ratio test= 7.99 on 1 df, p=0.005
Wald test = 7.84 on 1 df, p=0.005
Score (logrank) test = 8.07 on 1 df, p=0.004Curvas de supervivencia paramétricas vs Cox:
mu0.hat <- result.survreg.grp$coef[1] # intercept = -log λ₀
sigma.hat <- result.survreg.grp$scale # = 1/α
alpha.hat <- 1/sigma.hat
lambda0.hat <- exp(-mu0.hat)tt.vec <- 0:182
surv0.vec <- exp(- (lambda0.hat * tt.vec)^alpha.hat)gamma.hat <- result.survreg.grp$coef[2]
# hazard ratio = exp(-γ/σ)
hr <- exp(-gamma.hat/sigma.hat)
surv1.vec <- surv0.vec^hrcoxph.surv.est <- survfit(result.coxph.grp,
newdata=data.frame(grp=c("combination","patchOnly")))plot(coxph.surv.est, col=c("red","black"),
xlab="Días hasta recaída", ylab="Survival prob.")
lines(surv0.vec ~ tt.vec, col="red") # Weibull baseline
lines(surv1.vec ~ tt.vec, col="black") # Weibull patchOnly
En la Figura resultante se visualizan en step las curvas de Cox y en líneas suaves las estimaciones Weibull, mostrando su buena concordancia.
En el enfoque de regresión AFT Weibull, en lugar de trabajar directamente con funciones de riesgo, modelamos el log-tiempo de supervivencia como un modelo de localización-escala:
\[ \log(t)=\mu+\gamma\,z+\sigma\,\epsilon^{*}, \]
donde:
Bajo este modelo:
Esta formulación sugiere además que, si en lugar de \(\epsilon^{*}=\log\varepsilon\) (valor extremo) se elige otro error:
se obtienen otras familias paramétricas de supervivencia.
En los modelos Weibull AFT con varias covariables, se ajustan estadísticos de supervivencia de la misma forma que en Cox, pero usando survreg(..., dist="weibull"). Por ejemplo, en los datos pharmacoSmoking:
# Cox proporcional de referencia
modelAll2.coxph <- coxph(Surv(ttr, relapse) ~ grp + age + employment)
summary(modelAll2.coxph)Call:
coxph(formula = Surv(ttr, relapse) ~ grp + age + employment)
n= 125, number of events= 89
coef exp(coef) se(coef) z Pr(>|z|)
grppatchOnly 0.60788 1.83654 0.21837 2.784 0.00537 **
age -0.03529 0.96533 0.01075 -3.282 0.00103 **
employmentother 0.70348 2.02077 0.26929 2.612 0.00899 **
employmentpt 0.65369 1.92262 0.32732 1.997 0.04581 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
grppatchOnly 1.8365 0.5445 1.1971 2.8176
age 0.9653 1.0359 0.9452 0.9859
employmentother 2.0208 0.4949 1.1920 3.4256
employmentpt 1.9226 0.5201 1.0122 3.6518
Concordance= 0.638 (se = 0.03 )
Likelihood ratio test= 22.03 on 4 df, p=2e-04
Wald test = 21.91 on 4 df, p=2e-04
Score (logrank) test = 22.48 on 4 df, p=2e-04A continuación, el ajuste Weibull AFT con las mismas covariables:
model.pharm.weib <- survreg(
Surv(ttr, relapse) ~ grp + age + employment,
dist="weibull"
)
summary(model.pharm.weib)
Call:
survreg(formula = Surv(ttr, relapse) ~ grp + age + employment,
dist = "weibull")
Value Std. Error z p
(Intercept) 2.4024 0.9653 2.49 0.0128
grppatchOnly -1.1902 0.4133 -2.88 0.0040
age 0.0697 0.0203 3.43 0.0006
employmentother -1.3890 0.5029 -2.76 0.0057
employmentpt -1.3143 0.6132 -2.14 0.0321
Log(scale) 0.6313 0.0900 7.02 2.3e-12
Scale= 1.88
Weibull distribution
Loglik(model)= -454.1 Loglik(intercept only)= -466.1
Chisq= 23.96 on 4 degrees of freedom, p= 8.2e-05
Number of Newton-Raphson Iterations: 5
n= 125 Diferencias clave respecto a Cox:
grppatchOnly, age, … son parámetros AFT (\(\gamma_j\)): un valor negativo indica acortamiento de tiempos; el factor de aceleración es \(e^{\gamma_j}\).Bajo Weibull la correspondencia con riesgos proporcionales es
\[ \beta_j \;=\; -\frac{\gamma_j}{\sigma}, \]
donde \(\sigma=e^{\text{Log(scale)}}\). En R:
# Extraer gamma (AFT) y convertir a beta (PH)
weib.coef.all <- model.pharm.weib$coef
weib.gamma <- weib.coef.all[2:5] # grppatchOnly, age, ...
sigma.hat <- model.pharm.weib$scale
weib.beta.ph <- -weib.gamma / sigma.hat
# Coeficientes Cox
coxph.beta <- modelAll2.coxph$coef
# Comparar
data.frame(
predictor = names(coxph.beta),
beta_PH_weibull = weib.beta.ph,
beta_PH_cox = coxph.beta
) predictor beta_PH_weibull beta_PH_cox
grppatchOnly grppatchOnly 0.63301278 0.60788405
age age -0.03708786 -0.03528934
employmentother employmentother 0.73878031 0.70347664
employmentpt employmentpt 0.69903157 0.65369019Vemos que los \(\beta_j\) convertidos del AFT Weibull coinciden muy bien (diferencias < 7 %) con los del modelo de Cox.
En el contexto de un modelo Weibull con múltiples covariables, podemos reutilizar las herramientas de selección de modelo y diagnóstico por residuos vistas para Cox. Por ejemplo, ajustamos un modelo inicial con todas las covariables y aplicamos eliminación para atrás usando AIC:
modelAll.pharm.weib <- survreg(
Surv(ttr, relapse) ~ grp + gender + race + employment +
yearsSmoking + levelSmoking +
age + priorAttempts + longestNoSmoke,
dist="weibull"
)
model.step.pharm.weib <- step(modelAll.pharm.weib)Start: AIC=930.93
Surv(ttr, relapse) ~ grp + gender + race + employment + yearsSmoking +
levelSmoking + age + priorAttempts + longestNoSmoke
Df AIC
- race 3 928.34
- levelSmoking 1 928.96
- gender 1 929.08
- priorAttempts 1 929.14
- longestNoSmoke 1 929.49
- yearsSmoking 1 929.72
<none> 930.93
- age 1 936.20
- employment 2 936.37
- grp 1 936.94
Step: AIC=928.34
Surv(ttr, relapse) ~ grp + gender + employment + yearsSmoking +
levelSmoking + age + priorAttempts + longestNoSmoke
Df AIC
- gender 1 926.37
- levelSmoking 1 926.47
- priorAttempts 1 926.51
- yearsSmoking 1 926.76
- longestNoSmoke 1 927.59
<none> 928.34
- age 1 932.13
- employment 2 932.53
- grp 1 934.56
Step: AIC=926.37
Surv(ttr, relapse) ~ grp + employment + yearsSmoking + levelSmoking +
age + priorAttempts + longestNoSmoke
Df AIC
- levelSmoking 1 924.51
- priorAttempts 1 924.54
- yearsSmoking 1 924.80
- longestNoSmoke 1 925.61
<none> 926.37
- age 1 930.36
- employment 2 931.04
- grp 1 932.57
Step: AIC=924.51
Surv(ttr, relapse) ~ grp + employment + yearsSmoking + age +
priorAttempts + longestNoSmoke
Df AIC
- priorAttempts 1 922.69
- yearsSmoking 1 922.83
- longestNoSmoke 1 923.67
<none> 924.51
- age 1 928.38
- employment 2 929.21
- grp 1 930.71
Step: AIC=922.69
Surv(ttr, relapse) ~ grp + employment + yearsSmoking + age +
longestNoSmoke
Df AIC
- yearsSmoking 1 921.05
- longestNoSmoke 1 921.85
<none> 922.69
- age 1 926.64
- employment 2 927.28
- grp 1 928.75
Step: AIC=921.05
Surv(ttr, relapse) ~ grp + employment + age + longestNoSmoke
Df AIC
- longestNoSmoke 1 920.30
<none> 921.05
- employment 2 925.49
- grp 1 927.55
- age 1 929.48
Step: AIC=920.3
Surv(ttr, relapse) ~ grp + employment + age
Df AIC
<none> 920.30
- employment 2 925.38
- grp 1 926.86
- age 1 930.14La selección final coincide con el modelo de la sección anterior, conservando sólo grp, age y employment.
Para el análisis de residuos de devianza, usamos:
smoothSEcurve <- function(yy, xx) {
# use after a call to "plot"
# fit a lowess curve and 95% confidence interval curve
# make list of x values
xx.list <- min(xx) + ((0:100)/100)*(max(xx) - min(xx))
# Then fit loess function through the points (xx, yy)
# at the listed values
yy.xx <- predict(loess(yy ~ xx), se=T,
newdata=data.frame(xx=xx.list))
lines(yy.xx$fit ~ xx.list, lwd=2)
lines(yy.xx$fit -
qt(0.975, yy.xx$df)*yy.xx$se.fit ~ xx.list, lty=2)
lines(yy.xx$fit +
qt(0.975, yy.xx$df)*yy.xx$se.fit ~ xx.list, lty=2)
}
resid.deviance <- residuals(model.pharm.weib, type="deviance")
par(mfrow=c(2,2))
plot(resid.deviance ~ age, data=pharmacoSmoking)
smoothSEcurve(resid.deviance, age)
title("Deviance residuals\nversus age")
plot(resid.deviance ~ grp, data=pharmacoSmoking)
title("Deviance residuals\nversus treatment group")
plot(resid.deviance ~ employment, data=pharmacoSmoking)
title("Deviance residuals\nversus employment")
En la figura anterior, los residuos de desviancia en función de grp y employment están centrados alrededor de cero, y para age el patrón es consistente con un efecto lineal dado el ancho de los intervalos de confianza al 95 %.
Para identificar observaciones influyentes en el coeficiente de age, calculamos los dfbeta:
resid.dfbeta <- residuals(model.pharm.weib, type="dfbeta")
n.obs <- length(pharmacoSmoking$ttr)
index.obs <- 1:n.obs
plot(resid.dfbeta[, 3] ~ index.obs, type="h",
xlab="Observación", ylab="Cambio en coeficiente para age",
ylim=c(-0.0065, 0.004))
abline(h=0)
En la figura anterior aparecen nuevamente los pacientes 46 y 68 (y además el 114) como los que más alteran la estimación de \(\beta_{\text{age}}\), al igual que en el modelo de Cox.
En esta sección se muestran dos extensiones del modelo de tiempo acelerado (AFT) usando distribuciones distintas para \(\varepsilon\):
Si \(\varepsilon\sim N(0,1)\) entonces \(T\sim\mathrm{Log\text{-}Normal}\) y ajustamos en R con
model.pharm.lognormal <- survreg(
Surv(ttr, relapse) ~ grp + age + employment,
dist="lognormal"
)
summary(model.pharm.lognormal)
Call:
survreg(formula = Surv(ttr, relapse) ~ grp + age + employment,
dist = "lognormal")
Value Std. Error z p
(Intercept) 1.6579 1.0084 1.64 0.1002
grppatchOnly -1.2623 0.4523 -2.79 0.0053
age 0.0648 0.0203 3.20 0.0014
employmentother -1.1711 0.5316 -2.20 0.0276
employmentpt -0.9543 0.7198 -1.33 0.1849
Log(scale) 0.8754 0.0796 10.99 <2e-16
Scale= 2.4
Log Normal distribution
Loglik(model)= -451.5 Loglik(intercept only)= -461.7
Chisq= 20.4 on 4 degrees of freedom, p= 0.00042
Number of Newton-Raphson Iterations: 3
n= 125 Aquí los coeficientes son en la escala AFT: un efecto negativo significa menor tiempo esperado.
Si \(\varepsilon\sim\mathrm{Logistic}(0,1)\,,\quad S(u)=\frac1{1+e^u},\) entonces \(T\sim\mathrm{Log\text{-}Logistic}\) y ajustamos con
model.pharm.loglogistic <- survreg(
Surv(ttr, relapse) ~ grp + age + employment,
dist="loglogistic"
)
summary(model.pharm.loglogistic)
Call:
survreg(formula = Surv(ttr, relapse) ~ grp + age + employment,
dist = "loglogistic")
Value Std. Error z p
(Intercept) 1.9150 0.9708 1.97 0.0485
grppatchOnly -1.3260 0.4588 -2.89 0.0038
age 0.0617 0.0196 3.15 0.0017
employmentother -1.2605 0.5392 -2.34 0.0194
employmentpt -1.0991 0.7050 -1.56 0.1190
Log(scale) 0.3565 0.0884 4.03 5.5e-05
Scale= 1.43
Log logistic distribution
Loglik(model)= -453.4 Loglik(intercept only)= -463.6
Chisq= 20.47 on 4 degrees of freedom, p= 4e-04
Number of Newton-Raphson Iterations: 4
n= 125 En este caso la proporción de odds de supervivencia cumple
\[ \frac{S_1(t)}{1 - S_1(t)} \;=\; e^{z\beta}\;\frac{S_0(t)}{1 - S_0(t)}\,, \]
y el modelo AFT log-logístico es equivalente a un modelo de odds proporcionales, del mismo modo que el Weibull a PH.
Nota: todos los ajustes survreg(..., dist=...) devuelven coeficientes en la escala AFT.