Estos son los datos:
año <- 4:22
peso <- c(11.78, 18.43, 25.21, 30.78, 33.03, 35.66, 36.96, 37.97, 38.04, 39.2,
36.5, 37.21, 39.97, 38.45, 33.65, 34.71, 37.75, 32.81, 37.99)
plot(año, peso, pch = 16, col = "purple")
\[ p=\phi (t)=\frac{\alpha}{1+\beta e^{-kt}}+\varepsilon \]
Donde \( p= peso \) y \( t=año \) mientras que \( \left\{ \alpha,\beta,k \right\} \) son los parámetros a estimar.
\[ SCE=\sum_{i=1}^n\varepsilon_i^2=\sum_{i=1}^n \left (p_i- \phi(t_i) \right )^2=\sum_{i=1}^n\left (p_i - \frac{\alpha}{1+\beta e^{-kt_i}}\right ) ^2 \]
Donde:
\( \varepsilon_i= \) El error en la i-ésima observación.
\( p_i= \) El peso de la i-ésima observación
\( t_i= \) El año en la i-ésima observación
En este caso lo haremos en R, sin embargo se requiere introducir a la función valores iniciales de la estimación. Utilizaremos las recomendaciones de John Fox & Sanford Weisberg. Sin embargo el modelo en este texto es un poco diferente:
\[ p\approx \frac{\theta_1}{1+\exp[-(\theta_2+\theta_3 t)]} \]
Entonces podemos utilizar este nuevo modelo, utilizando las siguientes igualdades:
y a partir de aquí, podemos obtener:
\[ \ln\left[\frac {p/\theta_1}{1-y/\theta_1}\right]\approx \theta_2+\theta_3 t \]
Ahora es necesario únicamente una estimación de \( \theta_1 \), para el cual siguiendo el texto de Fox y Weisberg, podemos utilizar un valor mayor al máximo de \( p \), esto dado que el lado izquierdo de la ecuación es una razón de momios, i.e. una función logística, el \( \max(p)=39.97 \) y por lo tanto utilizaremos como punto de partida \( \hat{\theta}_1=40 \), utilizando este valor aplicamos una regresión lineal con los valores transformados utilizando la función logística:
library(boot)
lm(formula = logit(peso/40) ~ año)
##
## Call:
## lm(formula = logit(peso/40) ~ año)
##
## Coefficients:
## (Intercept) año
## 0.303 0.150
De tal manera que los puntos de partida para las estimaciones son:
Con estos puntos de partida, optimizamos la ecuación de sumas de cuadrados del error, pero primero es necesario definir la función a optimizar:
scuad <- function(estimadores) {
teta1 <- estimadores[1]
teta2 <- estimadores[2]
teta3 <- estimadores[3]
sum((peso - (teta1/(1 + exp(-(teta2 + (teta3 * año))))))^2)
}
Con la función definida en términos de los \( \theta \)'s, podemos optimizarla con la función optim() a partir de los puntos de inicio ya definidos:
inicio <- c(40, coef(lm(formula = logit(peso/40) ~ año))[1], coef(lm(formula = logit(peso/40) ~
año))[2])
names(inicio) <- c("teta1", "teta2", "teta3")
optim(par = inicio, fn = scuad, method = "BFGS")
## $par
## teta1 teta2 teta3
## 37.0806 -3.9163 0.7816
##
## $value
## [1] 56.15
##
## $counts
## function gradient
## 44 13
##
## $convergence
## [1] 0
##
## $message
## NULL
De tal manera que los valores que minimizan la suma de cuadrados del error son:
En este caso en lugar de utilizar SAS utilizamos R, si incluso la redacción del documento lo estamos haciendo en R no hay razón para cambiar, utilizaremos la función nls() que es para ajustar modelos no lineales por mínimos cuadrados, y utilizaremos los mismos puntos de partida que para la optimización de las suma de cuadrados del error:
mod_log <- nls(peso ~ teta1/(1 + exp(-(teta2 + (teta3 * año)))), start = inicio)
summary(mod_log)
##
## Formula: peso ~ teta1/(1 + exp(-(teta2 + (teta3 * año))))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## teta1 37.081 0.533 69.56 < 2e-16 ***
## teta2 -3.916 0.583 -6.72 5.0e-06 ***
## teta3 0.782 0.109 7.17 2.2e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.87 on 16 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 1.72e-06
Notemos que los estimadores de los coeficientes son iguales a los obtenidos utilizando la optimización en el punto anterior, sin embargo difieren con respecto al texto de Azme et. al., donde utilizando la transformación previamente mencionada tenemos que los estimadores son:
Es claro que el modelo estimado es diferente, veamos gráficamente que implica esta diferencia:
alfa0 <- 37.0806
alfa1 <- 4.8149
alfa2 <- 0.7817
pred2 <- alfa0/(1 + alfa1 * exp(-alfa2 * seq(from = 0, to = 24, by = 0.3)))
plot(año, peso, pch = 16, col = "purple")
lines(año, predict(mod_log), lty = 2, col = 2)
lines(seq(from = 0, to = 24, by = 0.3), pred2, lty = 2, col = 3)
legend("bottomright", bty = "n", lty = 2, col = 2:3, legend = c("Modelo de mínimos cuadrados",
"Modelo de Azme, et.al."))
A simple vista resulta evidente que los estimadores encontrados en este ejercicio tienden a minimizar de mejor manera los errores del modelo.
Recordemos que la R cuadrada está dada por:
\[ R^2=1-\frac{SCE}{SCT} \]
A partir de esta ecuación podemos obtener el \( R^2 \) utilizando la función diseñada previa de suma de cuadrados(scuad()):
\[ R^2_{lsq}=0.9435 \]
También lo podemos comparar con el modelo con los coeficientes del texto:
\[ R^2_{Azme}=0.1977 \]
Esto va acorde con lo visto gráficamente sin embargo, al revisar el artículo resulta que la SCE reportado es el mismo que el obtenido por nuestro modelo, lo cual hace pensar que más bien se trata de un problema de redacción.
La \( R^2 \) nos indica que el modelo explica el 94% de la dispersión de los datos, sin embargo la interpretación de los coeficientes en un poco más complicada, sobre todo al ser un modelo no linealizable, básicamente el interés residen en \( \theta_3 \) que es el que multiplica al tiempo, en este sentido por cada año que pase el aumento en el peso varía dependiendo del valor del año. Es decir que durante los primeros 10 años el peso muestra un crecimiento acelerado, que al llegar aproximadamente a los diez años se estabiliza en un valor cercano al 35.
Para este problema trabajaremos con datos sobre moscas de la fruta, el ejercicio es una adaptación del realizado por Hanley & Shapiro. Sin embargo las instrucciones no son muy claras, y me parece que hay un problema en cuanto se refiere a información con es proporcionada, y de hecho es lo que se lee, que no existe información al respecto:
Compliance' of the males in the two experimental groups was documented as follows: On two days per week throughout the life of each experimental male, the females that had been supplied as virgins to that male were kept and examined for fertile eggs. The insemination rate declined from approximately 7 females/day at age one week to just under 2/day at age eight weeks in the males supplied with eight virgin females per day, and from just under 1/day at age one week to approximately 0.6/day at age eight weeks in the males supplied with one virgin female per day. These 'compliance' data were not supplied for individual males, but the authors say that “There were no significant differences between the individual males within each experimental group.”
El hecho es que no se cuenta con la tas de inseminación de las hembras, los datos con los que se cuentan son:
| Variable | Descripción |
|---|---|
| ID | Serial No. (1-25) within each group of 25 (the order in which data points were abstracted) |
| PARTNERS | Number of companions (0, 1 or 8) |
| TYPE | Type of companion |
| * 0: newly pregnant female | |
| * 1: virgin female | |
| * 9: not applicable (when PARTNERS=0) | |
| LONGEVITY | Lifespan, in days |
| THORAX | Length of thorax, in mm (x.xx) |
| SLEEP | Percentage of each day spent sleeping |
Y como menciona el documento de Hanley & Shapiro el objetivo del estudio es observar el impacto que tienen las variables, particularmente las reproductivas, en el tiempo de vida de las moscas:
…an experiment designed to test if increased reproduction reduces longevity for male fruitflies. (Such a cost has already been established for females.)
De tal manera que a partir de esta idea aplicaremos un modelo que incluya interacciones entre las variables, pero veamos primero gráficamente de que se trata:
moscas <- read.csv("moscas.csv")
moscas$type <- as.factor(moscas$type)
par(mfrow = c(2, 2))
plot(moscas$parners, moscas$longevity, pch = 17, ylab = "Días", xlab = "partners",
col = 2)
plot(moscas$type, moscas$longevity, pch = 17, ylab = "Días", xlab = "type",
col = 3)
plot(moscas$thorax, moscas$longevity, pch = 17, ylab = "Días", xlab = "thorax",
col = 4)
plot(moscas$sleep, moscas$longevity, pch = 17, ylab = "Días", xlab = "sleep",
col = 6)
A partir de estas gráficas el único que parece claramente relacionado es el tamaño del tórax, pero justamente la idea es aplicar un modelo de regresión lineal múltiple para observar estas situaciones, sin embargo dado el número de variables, que interesan las interacciones entre ellas, y que la muestra no es tan grande, iniciaremos con un modelo sin interacciones, y a partir de este modificaremos para llegar a un modelo más manejable, esto a riesgo de inflar el error tipo 1, así como caer en un sobre ajuste, pero como de moscas yo se lo mismo que de osos bailarines me parece que es la mejor opción:
mod_lin_1 <- lm(longevity ~ parners + type + thorax + sleep, data = moscas)
drop1(mod_lin_1, test = "F")
## Single term deletions
##
## Model:
## longevity ~ parners + type + thorax + sleep
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 14292 604
## parners 1 833 15125 609 6.93 0.0096 **
## type 2 7119 21411 651 29.64 3.6e-11 ***
## thorax 1 13735 28027 687 114.36 < 2e-16 ***
## sleep 1 111 14403 603 0.93 0.3376
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Utilizamos para decidir el ajuste del modelo un análisis de varianza a partir de una suma de cuadrados tipo III, donde no importa el orden de inclusión de las variables en el modelo, con la tabla ANOVA anterior concluimos que la variable de sueño no resulta significativa, y ajustamos un nuevo modelo incluyendo las interacciones de las variables que resultaron significativas, sin embargo hay que tener en consideración que existe una relación lógica determinística entre la clasificación de las variables, (véase la descripción), de tal manera que partners y type no pueden interactuar dado que en al menos un valor posible de partners¨ determina el valor de type:
mod_lin_2 <- lm(longevity ~ parners * thorax + type * thorax, data = moscas)
anova(mod_lin_2)
## Analysis of Variance Table
##
## Response: longevity
## Df Sum Sq Mean Sq F value Pr(>F)
## parners 1 3513 3513 28.56 4.5e-07 ***
## thorax 1 13271 13271 107.88 < 2e-16 ***
## type 2 7065 3533 28.72 7.1e-11 ***
## parners:thorax 1 0 0 0.00 0.97
## thorax:type 2 10 5 0.04 0.96
## Residuals 117 14393 123
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
En este caso utilizamos una suma de cuadrados tipo I que nos indica que las interacciones no resultan significativas según el orden de inclusión de la variable en el modelo, probaremos mediante F's parciales (suma de cuadrados tipo I), si las variables resultan significativas a pesar del orden de inclusión en el modelo:
anova(lm(longevity ~ thorax + type + parners, data = moscas))
## Analysis of Variance Table
##
## Response: longevity
## Df Sum Sq Mean Sq F value Pr(>F)
## thorax 1 15497 15497 129.11 < 2e-16 ***
## type 2 7545 3773 31.43 1.1e-11 ***
## parners 1 808 808 6.73 0.011 *
## Residuals 120 14403 120
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(longevity ~ type + parners + thorax, data = moscas))
## Analysis of Variance Table
##
## Response: longevity
## Df Sum Sq Mean Sq F value Pr(>F)
## type 2 7845 3923 32.7 4.7e-12 ***
## parners 1 2372 2372 19.8 2.0e-05 ***
## thorax 1 13633 13633 113.6 < 2e-16 ***
## Residuals 120 14403 120
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod_lin_3 <- lm(longevity ~ thorax + type + parners, data = moscas)
A partir de esto concluimos con el modelo que incluye las tres variables: Partners, Type y Thorax. Ahora probaremos los supuestos de este modelo
Si bien apoyándonos en el tamaño de la muestra y el teorema central del límite podemos asumir normalidad en la media de la longevidad, vale la pena probar este supuesto, hay dos maneras para realizarlo, gráficamente y mediante una prueba estadística:
qqnorm(mod_lin_3$residuals)
qqline(mod_lin_3$residuals, col = 2)
shapiro.test(mod_lin_3$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod_lin_3$residuals
## W = 0.9899, p-value = 0.4984
Tanto la gráfica de normalidad, como la prueba de Shapiro Wilks, indican que no es posible rechazar que los residuos se distribuyen según una distribución normal, esto con un 95% de confianza. De tal manera que el supuesto de normalidad se cumple.
Básicamente un outlier es un punto que se encuentra lejos del resto, veamos gráficamente los residuales estudentizados,adicionalmente aplicamos la prueba de Bonferroni para residuos estudentizados:
library(car)
## Warning: package 'car' was built under R version 3.0.3
##
## Attaching package: 'car'
##
## The following object is masked from 'package:boot':
##
## logit
plot(rstudent(mod_lin_3), xlab = "", ylab = "Residuales estudentizados")
abline(h = 0, lty = 2, col = 2)
outlierTest(mod_lin_3)
##
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferonni p
## 67 -2.853 0.005114 0.6392
Desde que observamos gráficamente vemos que los residuales estudentizados se mantienen dentro de los límites de \( |3| \), adicionalmente sabemos que el residual más grande es de 2.5267, a partir de esto podemos prácticamente concluir que no existen outliers; sin embargo la prueba Bonferroni nos indica que no existe ningún dato que podamos clasificar como outlier.
Para este caso graficaremos la distancia de Cook:
plot(cooks.distance(mod_lin_3), xlab = "", ylab = "Distancia de Cook")
Se muestra que en general los valores son pequeños, adicionalmente el valor más alto es de 0.054 por lo cual concluimos que no existen problemas de puntos de influencia.
Ahora probaremos la falta de linealidad, mediante gráficas de residuales contra las variables continuas en el modelo:
par(mfrow = c(1, 2))
plot(moscas$parners, rstudent(mod_lin_3), ylab = "Residuos Estudentizados",
xlab = "Partner", col = 11, pch = 3)
abline(h = 0, col = 2, lty = 2)
plot(moscas$thorax, rstudent(mod_lin_3), ylab = "Residuos Estudentizados", xlab = "Partner",
col = 12, pch = 3)
abline(h = 0, col = 2, lty = 2)
No observamos ninguna tendencia en el comportamiento de los residuales que nos hagan pensar en un comportamiento no lineal, claro que dado el diseño del experimento queda la pregunta si no sería mejor utilizar la variable de número de parejas como una variable categórica, sin embargo insistiré sobre la parte que de moscas conozco que son negras y vuelan. Concluimos que no existen problemas de falta de linealidad en el modelo.
La colinealidad se puede observar considerando la correlación entre las variables explicativas en el modelo, en este caso únicamente tenemos dos variables numéricas cuyo \( \rho_{partners,thorax}=-0.1933 \), que nos indica que la correlación entre estas dos variables es poca, y nos da a pensar que no existe colinealidad; sin embargo existen otras maneras de verificar colinealidad o falta de independencia, por ejemplo el indicador vif que es el indicador de la inflación de la varianza, el cual utilizaremos
sqrt(vif(mod_lin_3))
## GVIF Df GVIF^(1/(2*Df))
## thorax 1.019 1.000 1.010
## type 1.154 1.414 1.036
## parners 1.170 1.000 1.082
Al obtener valores cercanos a 1 concluimos que no existen problemas de colinealidad.
Ahora probaremos la igualdad de varianzas, nuevamente lo podemos realizar mediante una gráfica de residuales, contra valores ajustados, adicionalmente se realiza una prueba de scores sobre no constancia del la varianza del error (Breusch-Pagan test):
plot(mod_lin_3$fitted, rstudent(mod_lin_3), ylab = "Residuos estudentizados",
xlab = "Valores ajustados de la regresión")
abline(h = 0, col = 2, lty = 2)
ncvTest(mod_lin_3)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 7.213 Df = 1 p = 0.007237
Concluimos que existen problemas de igualdad de varianzas, podemos suponer a partir de las gráficas previas que de hecho se debe a que conforme el tamaño del tórax aumenta las desviaciones del modelo también lo hacen.
Dado que no encontramos colinealidad en el modelo pero sí heterocedsaticidad, optamos para solucionar este problema por la versión particular de la suma de cuadrados generalizados (GLS) que es la suma de cuadrados ponderados (WLS) vista en clase, dado que consideramos que el cambio de varianza se debe a la variable del tamaño del tórax, utilizaremos esta variable para determinar el peso. En la siguiente gráfica se observa la relación entre el tamaño del error cuadrado (desviación) y el tamaño del tórax:
library(grDevices)
plot(moscas$thorax, (rstudent(mod_lin_3))^2, xlab = "Thorax", ylab = expression(epsilon^2))
abline(lm((rstudent(mod_lin_3))^2 ~ moscas$thorax), col = 2)
A partir de esta interpretación de los errores utilizaremos la variable thorax como los pesos de la regresión, obteniendo el siguiente modelo:
mod_lin_p <- lm(longevity ~ parners + type + thorax, data = moscas, weights = thorax)
summary(mod_lin_p)
##
## Call:
## lm(formula = longevity ~ parners + type + thorax, data = moscas,
## weights = thorax)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -27.942 -5.884 -0.551 5.885 24.780
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -44.027 11.610 -3.79 0.00024 ***
## parners -0.827 0.322 -2.57 0.01149 *
## type1 -16.723 2.226 -7.51 1.1e-11 ***
## type9 -7.029 3.055 -2.30 0.02315 *
## thorax 137.092 13.595 10.08 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.1 on 120 degrees of freedom
## Multiple R-squared: 0.604, Adjusted R-squared: 0.591
## F-statistic: 45.8 on 4 and 120 DF, p-value: <2e-16
El modelo explica aproximadamente el 60% de la variabilidad de los datos, en cuanto a los coeficientes de estos, por cada pareja (partner) que la mosca macho tenga su longevidad en días se reduce en 0.8 días; mientras que se reduce en 16 días promedio cuando el tipo de compañeras es de hembras vírgenes, frente a hembras recientemente preñadas; y solo se reduce en 7 días cuando las parejas no existen, siendo pues el caso donde comparte con hembras vírgenes cuando la longevidad es mas afectada, y cuando se trata de hembras recientemente preñadas el caso donde la longevidad es mayor. El tamaño del tórax, es decir la fisonomía de la mosca es también significativa, por cada milímetro que aumente el tórax de la mosca su vida promedio aumenta en 137 días, es decir que por cada, décima de milímetro la vida promedio aumenta en casi 14 días.
En conclusión, los datos indican que las moscas robustas viven más, y que entre menos actividad reproductiva tengan también viven más; algo así como que las hembras les chupan la vida, pero eso no es políticamente correcto. Además no existe interacción entre estos aspectos, y el tiempo de sueño no es relevante para la vida esperada de una mosca…nunca he visto una mosca dormir…o al menos no he sabido.
En este último ejercicio vamos a utilizar técnicas de análisis de supervivencia a un caso de recuperación de heroinómanos,tomados del artículo de Luis Palmer Pol, donde los datos son:
| Variable | Descripción |
|---|---|
| tiempo | El tiempo que un heroinómano ha estado incluido en el programa, medido en semanas |
| estado | Situación del sujeto en su última observación |
| 0= Sano | |
| 1= Recaída | |
| grupo | Diferencia a los sujetos según el tipo de desintoxicación seguida |
| 0=Psico: Con soporte psico-social | |
| 1=Farma: Con medicación | |
| edad | La edad del sujeto en años |
Es evidente que la variable de interés, y que adicionalmente determina la censura por la derecha es el estado, es decir que en este caso el recaer en la heroína es equivalente a morir…algo así como la canción de The Velvet Underground.
Se requieren varias librerías
library(MASS)
library(graphics)
library(lattice)
library(foreign)
library(survival)
library(splines)
library(KMsurv)
require(grDevices)
library(stats)
library(mvtnorm)
library(modeltools)
library(stats4)
library(coin)
tiempo <- c(3, 5, 5, 8, 8, 13, 16, 19, 20, 20, 21, 25)
estado <- c(1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0)
grupo <- c(1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0)
edad <- c(45, 20, 39, 51, 30, 35, 44, 28, 35, 35, 38, 24)
Iniciaremos con un análisis de curvas de supervivencia de Kaplan-Meier, para lo cual es necesario formar los objetos Survival:
t.surv <- Surv(tiempo, estado)
km <- survfit(t.surv ~ grupo + edad, type = "kaplan-meier", conf.int = 0)
km1 <- survfit(t.surv ~ grupo, type = "kaplan-meier", conf.int = 1)
plot(km1, col = c("Purple", "Red"), ylab = "Probabilidad de supervivencia S(t)",
xlab = "Semanas", main = "Curvas de Kaplan-Meier, para grupo")
legend("bottomleft", bty = "n", col = c("Purple", "Red"), lty = 1, legend = c("Con soporte psico-social",
"Con medicanto"))
Si bien tenemos dos variables de interés: el grupo y la edad, dada la naturaleza de las curvas de KM, únicamente introducimos la variable categórica de grupo; gráficamente se muestra como el tener soporte psicosocial aumenta la probabilidad de supervivencia frente al tratamiento médico; sin embargo es prudente utilizar pruebas de hipótesis para verificar este supuesto, así como la relación entre el tiempo de supervivencia y la edad de las personas:
survdiff(t.surv ~ grupo, rho = 0)
## Call:
## survdiff(formula = t.surv ~ grupo, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## grupo=0 6 3 4.5 0.5 2.61
## grupo=1 6 3 1.5 1.5 2.61
##
## Chisq= 2.6 on 1 degrees of freedom, p= 0.106
survdiff(t.surv ~ edad, rho = 0)
## Call:
## survdiff(formula = t.surv ~ edad, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## edad=20 1 0 0.0833 0.0833 0.0909
## edad=24 1 0 1.2222 1.2222 1.8093
## edad=28 1 0 0.4722 0.4722 0.5746
## edad=30 1 1 0.3056 1.5783 1.9350
## edad=35 3 1 1.7500 0.3214 0.5364
## edad=38 1 1 1.2222 0.0404 0.0598
## edad=39 1 0 0.0833 0.0833 0.0909
## edad=44 1 1 0.4722 0.5899 0.7177
## edad=45 1 1 0.0833 10.0833 11.0000
## edad=51 1 1 0.3056 1.5783 1.9350
##
## Chisq= 19.3 on 9 degrees of freedom, p= 0.0231
A partir de la prueba log-rank, concluimos que no existe evidencia para rechazar que las curvas de supervivencia de ambos grupos (tratamientos) sean iguales, esto con un 95% de confianza. Es decir que las curvas no son diferentes, i.e. no hay evidencia de que los tratamientos difieran en cuanto al tiempo de supervivencia, es decir el tiempo en que tienen una recaída. Al considerar la edad, observamos que a partir de la prueba podemos rechazar con un 95% de confianza, sin embargo hay que considerar que en este caso lo que nos dice la prueba es que existe al menos un para de edades para las cuales su curva de supervivencia es diferente, no nos indica cuales son las que son diferentes. Podríamos buscar aplicar pruebas pareadas, sin embargo dado el tamaño de la muestra y la naturaleza de la variable edad esto no resulta conveniente. En su lugar es más conveniente utilizar una regresión semiparamétrica de Cox, un modelo de riesgos proporcionales:
mod_ph <- coxph(t.surv ~ grupo + edad)
summary(mod_ph)
## Call:
## coxph(formula = t.surv ~ grupo + edad)
##
## n= 12, number of events= 6
##
## coef exp(coef) se(coef) z Pr(>|z|)
## grupo 0.7006 2.0149 1.4335 0.49 0.63
## edad 0.1157 1.1226 0.0804 1.44 0.15
##
## exp(coef) exp(-coef) lower .95 upper .95
## grupo 2.01 0.496 0.121 33.45
## edad 1.12 0.891 0.959 1.31
##
## Concordance= 0.75 (se = 0.158 )
## Rsquare= 0.367 (max possible= 0.83 )
## Likelihood ratio test= 5.49 on 2 df, p=0.0643
## Wald test = 4.65 on 2 df, p=0.0978
## Score (logrank) test = 5.99 on 2 df, p=0.0501
cox.zph(mod_ph)
## rho chisq p
## grupo 0.00526 5.38e-05 0.994
## edad 0.30918 2.48e-01 0.619
## GLOBAL NA 3.76e-01 0.829
La función summary() aplicada a un objeto del tipo modelo de riesgos proporcionales, además de regresar los estimadores de los coeficientes nos presenta una prueba de Wald para cada coeficiente, así como tres pruebas simultaneas sobre la significancia del modelo:Razón de verosimilitud, Wald y Score-logrank. Ninguna de estas pruebas presenta un p-valor menor a 0.05 por lo cual concluimos que el modelo con las variables edad y grupo resulta significativo. Adicionalmente se probó el supuesto de riesgos proporcionales, que para ambas variables se cumple.
Sin embargo en el artículo únicamente se prueba la edad en el modelo de Cox, y ciertamente lo que sucede con este tipo de modelos es que se prueba si son significativas, dado que las otras variables se encuentran dentro del modelo por lo cual, probaremos de manera separa cada variable.
coxph(t.surv ~ grupo)
## Call:
## coxph(formula = t.surv ~ grupo)
##
##
## coef exp(coef) se(coef) z p
## grupo 1.67 5.3 1.16 1.43 0.15
##
## Likelihood ratio test=2.45 on 1 df, p=0.117 n= 12, number of events= 6
coxph(t.surv ~ edad)
## Call:
## coxph(formula = t.surv ~ edad)
##
##
## coef exp(coef) se(coef) z p
## edad 0.142 1.15 0.0683 2.08 0.038
##
## Likelihood ratio test=5.25 on 1 df, p=0.0219 n= 12, number of events= 6
Ciertamente al incluir únicamente la variable de edad en el modelo, este es significativo con un 95% de confianza. La interpretación del modelo y su coeficientes consiste en que por cada año de edad que el paciente tenga más con respecto a otro, la tasa de riesgo es 1.15 veces la del segundo, es decir 0.15 veces mayor. A manera más sencilla entre más grande sea la persona corre un mayor riesgo de recaer en la heroína. Faltaría preguntarse sobre otra variable que quizá sea más explicativa, que sería el tiempo que tiene la persona consumiendo heroína, y que claramente guarda relación con la de la edad de la persona, pero en términos explicativos puede aportar más.