Comparando un modelo lineal con un modelo spline

Question

Trato de averiguar cómo describir mi variable continua. Por desgracia, yo no entiendo todas las estadísticas. Realmente agradecería que me ustedes podría ayudarme a salir de aquí. Para ilustrar mejor mi problema, me escribió el siguiente ejemplo.

library(rms)
library(survival)

data(pbc)
d <- pbc
rm(pbc, pbcseq)
d$status <- ifelse(d$status != 0, 1, 0)

dd = datadist(d)
options(datadist='dd')

# linear model
f1 <- cph(Surv(time, status) ~  albumin, data=d)
p1 <- Predict(f1, fun=exp)
(a1 <- anova(f1))
Function(f1)
plot(p1, anova=a1, pval=TRUE, ylab="Hazard Ratio")

# rcs model
f2 <- cph(Surv(time, status) ~  rcs(albumin, 4), data=d)
p2 <- Predict(f2, fun=exp)
(a2 <- anova(f2))
Function(f2)
plot(p2, anova=a2, pval=TRUE, ylab="Hazard Ratio")

# minimal CI width
p1$diff <- p1$upper-p1$lower
    min(p1$diff) # = 0.002321521
p1[which(p1$diff==min(p1$diff)),]$albumin # = 3.494002
    describe(d$albumin) # mean = 3.497

p2$diff <- p2$upper-p2$lower
    min(p2$diff) # = 0.2039817
p2[which(p2$diff==min(p2$diff)),]$albumin # = 3.502447
    describe(d$albumin) # mean = 3.497

# both models in a single figure
p <- rbind(linear.model=p1, rcs.model=p2)
library(ggplot2)
df <- data.frame(albumin=p$albumin, yhat=p$yhat, lower=p$lower, upper=p$upper, predictor=p$.set.)
(g <- ggplot(data=df, aes(x=albumin, y=yhat, group=predictor, color=predictor)) + geom_line(size=1))
(g <- g + geom_ribbon(data=df, aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0))
(g <- g + theme_bw())
(g <- g + xlab("Albumin"))
(g <- g + ylab("Hazard Ratio"))
(g <- g + theme(axis.line = element_line(color='black', size=1)))
(g <- g + theme(axis.ticks = element_line(color='black', size=1)))
(g <- g + theme( plot.background = element_blank() ))
(g <- g + theme( panel.grid.minor = element_blank() ))
(g <- g + theme( panel.border = element_blank() ))

1. Por qué se muestra la trama de un modelo lineal (p1) no es una línea recta?
2. ¿Cómo puedo trazar los modelos de f1 y f2 en la misma figura?
3. ¿Cómo puedo comparar los modelos de f1 y f2 para investigar los modelos que se ajusta a los datos mejor? ... como anova() para coxph en el paquete de supervivencia
4. ¿Por qué es la mínima CI ancho cerca de la media de la albúmina más pronunciado en el lineal (f1) modelo?
5. ¿Qué significa el valor de P en las parcelas decir? ¿Cómo tengo que interpretar la salida de anova(...)

Actualización #1 Tras la respuesta de Harrell, he actualizado el código anterior muestra cómo combinar spline parcelas de dos predictores en una sola figura. Una última pregunta: ¿Cómo puedo comparar los dos rms modelos como anova(m1, m2) de supervivencia del paquete como se muestra a continuación?

> m1 <- coxph(Surv(time, status) ~ albumin, data=d)
> m2 <- coxph(Surv(time, status) ~ pspline(albumin), data=d)
> anova(m1, m2) # compare models
Analysis of Deviance Table
 Cox model: response is  Surv(time, status)
 Model 1: ~ albumin
 Model 2: ~ pspline(albumin)
   loglik  Chisq Df P(>|Chi|)
1 -975.61                    
2 -973.26 4.6983 11    0.9449
> summary(m1)
Call:
coxph(formula = Surv(time, status) ~ albumin, data = d)

  n= 418, number of events= 186 

           coef exp(coef) se(coef)      z Pr(>|z|)    
albumin -1.4695    0.2300   0.1714 -8.574   <2e-16 ***
---
Signif. codes:  0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1

        exp(coef) exp(-coef) lower .95 upper .95
albumin      0.23      4.347    0.1644    0.3219

Concordance= 0.688  (se = 0.023 )
Rsquare= 0.147   (max possible= 0.992 )
Likelihood ratio test= 66.6  on 1 df,   p=3.331e-16
Wald test            = 73.51  on 1 df,   p=0
Score (logrank) test = 72.38  on 1 df,   p=0

ACTUALIZACIÓN #2 Creo que me acaba de responder a mi "la última pregunta" por mí mismo (ver más abajo). Espero que esto no mostrar correcta accidentalmente. Me gustaría pensar que puedo comparar los modelos de cph y coxph de que manera, yo no puedo? Es la forma de calcular los grados de libertad df correcta?

> # using coxph from survival
> m1 <- coxph(Surv(time, status) ~  albumin, data=d)
> m2 <- coxph(Surv(time, status) ~  albumin + age, data=d)
> # loglik = a vector of length 2 containing the log-likelihood with the initial values and with the final values of the coefficients.
> m1$loglik[2]
    [1] -975.6126
    > m2$loglik[2]
[1] -973.2272
> (df <- abs(length(m1$coefficients) - length(m2$coefficients)))
[1] 1
> (LR <- 2 * (m2$loglik[2] - m1$loglik[2]))
[1] 4.770787
> pchisq(LR, df, lower=FALSE)
[1] 0.02894659
> anova(m2, m1)
Analysis of Deviance Table
 Cox model: response is  Surv(time, status)
 Model 1: ~ albumin + age
 Model 2: ~ albumin
   loglik  Chisq Df P(>|Chi|)  
1 -973.23                      
2 -975.61 4.7708  1   0.02895 *
---
Signif. codes:  0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1
> m1 <- cph(Surv(time, status) ~  albumin, data=d)
> m2 <- cph(Surv(time, status) ~  albumin + age, data=d)
> # loglik = a vector of length 2 containing the log-likelihood with the initial values and with the final values of the coefficients.
> m1$loglik[2]
    [1] -975.6126
    > m2$loglik[2]
[1] -973.2272
> (df <- abs(length(m1$coefficients) - length(m2$coefficients)))
[1] 1
> (LR <- 2 * (m2$loglik[2] - m1$loglik[2]))
[1] 4.770787
> pchisq(LR, df, lower=FALSE)
[1] 0.02894659

ACTUALIZACIÓN #3 He cambiado el ejemplo siguiente el tipo de respuesta de DWin de la siguiente manera. De esta manera los grados de libertad debe ser calculado correctamente:

library(Hmisc)
library(rms)
library(ggplot2)
library(gridExtra)

data(pbc)
d <- pbc
rm(pbc, pbcseq)
d$status <- ifelse(d$status != 0, 1, 0)

### log likelihood test using a coxph model
m1 <- coxph(Surv(time, status) ~  albumin, data=d)
m2 <- coxph(Surv(time, status) ~  albumin + age, data=d)
# loglik = a vector of length 2 containing the log-likelihood with the initial values and with the final values of the coefficients.
m1$loglik[2]
    m2$loglik[2]
(df <- abs(sum(anova(m1)$Df, na.rm=TRUE) - sum(anova(m2)$Df, na.rm=TRUE)))
(LR <- 2 * (m2$loglik[2] - m1$loglik[2])) # the most parsimonious models have to be first
pchisq(LR, df, lower=FALSE)
anova(m2, m1)

### log likelihood test using a cph model
dd = datadist(d)
options(datadist='dd')
m3 <- cph(Surv(time, status) ~  albumin, data=d)
m4 <- cph(Surv(time, status) ~  albumin + age, data=d)
# loglik = a vector of length 2 containing the log-likelihood with the initial values and with the final values of the coefficients.
m3$loglik[2]
    m4$loglik[2]
(df <- abs(print(anova(m3)[, "d.f."])[['TOTAL']] - print(anova(m4)[, "d.f."])[['TOTAL']]))
(LR <- 2 * (m4$loglik[2] - m3$loglik[2])) # the most parsimonious models have to be first
pchisq(LR, df, lower=FALSE)

Answer 1

No tengo tiempo para responder a todas sus preguntas. Para responder a dos de ellos: el $P$-valor para el test de todos los coeficientes correspondientes a un determinado predictor. En otras palabras es la influencia global de todos los términos, ya sea lineal o no lineal o interacciones con ellos. Esto es una prueba de asociación, donde la hipótesis nula de no asociación (planitud). Este particular el cálculo se utiliza la Wald $\chi^2$ prueba, la cual es una generalización de la $z$-prueba.

Usted obtendrá el mejor producto con ggplot(Predict(...)) en lugar de plot(). Para cualquiera de los dos puede combinar los resultados de Predict() el uso de la rbind función. Tipo ?rbind.Predict ayuda.

Answer 2

El "p-valor" podría ser llamado un "muy poco probable estadísticas". Se supone que debe ser la medida en que los datos no está de acuerdo con el valor Null o "modelo alternativo". TYu son en realidad la realización de una "prueba de razón de verosimilitud" cuando se comparan dos (anidada) modelos. A pesar de que usted ha elegido utilizar la pspline función que toma un número bastante elevado de nudos, probablemente no importa mucho en este caso. En un verdadero análisis de datos podría tener sentido para limitar el número de nudos. Me parece que más de 4 nudos a una spline plazo parece añadir muy poco en el significado adicional. Tanto la diferencia en la desviación o ( -2*(LL1-LL0) ) contribuir a un t-estadístico y el correspondiente p-valor. El más nudos, cuanto mayor es el obstáculo más que la prueba de razón de verosimilitud debe saltar para alcanzar convencional niveles de "importancia".

Estoy teniendo un poco de dificultad con el anova sobre el pspline modelo en comparación con la no-spline modelo. Si sólo tiene una diferencia de -975.61 - -973.26) no sé por qué el estadístico de chi-cuadrado es 4.6983. Parece como si debería ser un poco más de 2. En cualquier caso, la df=11 va a impedir la conclusión de tat la spline modelo es significativamente mejor, ya que el valor crítico de chisq para una diferencia de 10 df es de alrededor de 20.

Respuesta a P en el comentario.: El df de la anova-objeto de la coxph modelos es

anova(m1,m2)$Df

El df de la f2 objeto puede ser visto con:

a2[, "d.f."]

(Para examinar el R objetos puede print de ellos o la utilización str )