R cox.ph 附加 Cox 比例風險模型

說明

cox.ph 係列實現了 Cox 比例風險模型，其中包含 Peto 對關係的校正、可選分層以及通過懲罰部分似然最大化進行的估計，與 gam 一起使用。在模型公式中，事件時間是響應。在分層下，響應有兩列：時間和層的數字索引。 weights 向量提供審查信息(0 表示審查，1 表示事件)。 cox.ph 處理每個受試者都有一個事件/審查時間和一行協變量值的情況。當每個受試者有多個時間相關協變量時，請參閱cox.pht。

請參閱下麵的條件邏輯回歸示例。

用法

cox.ph(link="identity")

參數

link	目前(可能永遠)僅支持"identity"。

細節

與 gam 一起使用，使 Cox 比例風險模型適合生存數據。模型公式的左側將包含事件/審查時間，右側將包含線性預測器規範。審查信息由 weights 參數提供給 gam ，其中 1 表示事件，0 表示審查。

分層是可能的，允許不同層存在不同的基線危險。在這種情況下，響應有兩列：第一列是事件/審查時間，第二列是數字層索引。請參閱下麵的示例。

使用 type="response" 從擬合模型對象(使用 predict 方法)進行的預測將在幸存者函數尺度上進行預測。這需要提供評估時間以及協變量(參見示例)。另請參閱下麵的示例代碼，以直接提取累積基線危險/生存率。模型對象中存儲的 fitted.values 是每個受試者在事件/審查時間的生存函數估計。

可以提取 deviance 、 martingale 、 score 或 schoenfeld 殘差。參見 Klein amd Moeschberger (2003) 的說明。分數殘差以與模型矩陣相同維度的矩陣形式返回，具有 "terms" 屬性，該屬性是一個列表，指示哪些模型矩陣列屬於哪些模型項。分數殘差被縮放。對於參數項，這是由相關模型係數的標準偏差決定的。對於平滑項，該項的分數殘差子矩陣後乘以該項係數的協方差矩陣的轉置 Cholesky 因子。這是一種使係數近似獨立的變換，根據需要使分數殘差與事件時間的關係圖可解釋以檢查比例風險假設(參見 Klein amd Moeschberger，2003，p376)。懲罰會導致分數殘差發生漂移，該漂移也被刪除，以允許殘差近似解釋為未懲罰的分數殘差。 Schoenfeld 和分數殘差按分層計算。有關詳細信息，請參閱通過繪製分數殘差進行簡單 PH 假設檢查的示例，以及 Klein amd Moeschberger(2003 年，第 11.4 節)。請注意，術語之間的高度相關性可能會破壞這些檢查。

模型係數的估計是通過最大化受平滑懲罰懲罰的 log-partial 似然來實現的。參見例如Hastie 和 Tibshirani，1990 年，第 8.3 節。對於所使用的部分似然(使用 Peto 近似關係)，但請注意，部分似然的優化並不遵循 Hastie 和 Tibshirani。有關殘差、累積基線風險、生存函數和相關標準誤差的估計，請參閱 Klein amd Moeschberger (2003)(生存標準誤差表達式有一個拚寫錯誤)。

Cox PH 模型報告的偏差百分比解釋基於偏差殘差的平方和(作為模型偏差)以及協變量效應設置為零時的偏差殘差的平方和(作為零偏差)。兩者使用相同的基線危險估計。

該係列有效地處理每個受試者具有一個事件/審查時間和一行協變量值的情況。對於每個受試者有多個時變協變量測量的研究，則應使用 bam(...,discrete=TRUE) 將等效泊鬆模型擬合到合適的偽數據。請參閱cox.pht。

值

繼承自類 general.family 的對象。

例子

library(mgcv)
library(survival) ## for data
col1 <- colon[colon$etype==1,] ## concentrate on single event
col1$differ <- as.factor(col1$differ)
col1$sex <- as.factor(col1$sex)

b <- gam(time~s(age,by=sex)+sex+s(nodes)+perfor+rx+obstruct+adhere,
         family=cox.ph(),data=col1,weights=status)

summary(b) 

plot(b,pages=1,all.terms=TRUE) ## plot effects

plot(b$linear.predictors,residuals(b))

## plot survival function for patient j...

np <- 300;j <- 6
newd <- data.frame(time=seq(0,3000,length=np))
dname <- names(col1)
for (n in dname) newd[[n]] <- rep(col1[[n]][j],np)
newd$time <- seq(0,3000,length=np)
fv <- predict(b,newdata=newd,type="response",se=TRUE)
plot(newd$time,fv$fit,type="l",ylim=c(0,1),xlab="time",ylab="survival")
lines(newd$time,fv$fit+2*fv$se.fit,col=2)
lines(newd$time,fv$fit-2*fv$se.fit,col=2)

## crude plot of baseline survival...

plot(b$family$data$tr,exp(-b$family$data$h),type="l",ylim=c(0,1),
     xlab="time",ylab="survival")
lines(b$family$data$tr,exp(-b$family$data$h + 2*b$family$data$q^.5),col=2)
lines(b$family$data$tr,exp(-b$family$data$h - 2*b$family$data$q^.5),col=2)
lines(b$family$data$tr,exp(-b$family$data$km),lty=2) ## Kaplan Meier

## Checking the proportional hazards assumption via scaled score plots as
## in Klein and Moeschberger Section 11.4 p374-376... 

ph.resid <- function(b,stratum=1) {
## convenience function to plot scaled score residuals against time,
## by term. Reference lines at 5% exceedance prob for Brownian bridge
## (see KS test statistic distribution).
  rs <- residuals(b,"score");term <- attr(rs,"term")
  if (is.matrix(b$y)) {
    ii <- b$y[,2] == stratum;b$y <- b$y[ii,1];rs <- rs[ii,]
  }
  oy <- order(b$y)
  for (i in 1:length(term)) {
    ii <- term[[i]]; m <- length(ii)
    plot(b$y[oy],rs[oy,ii[1]],ylim=c(-3,3),type="l",ylab="score residuals",
         xlab="time",main=names(term)[i])
    if (m>1) for (k in 2:m) lines(b$y[oy],rs[oy,ii[k]],col=k);
    abline(-1.3581,0,lty=2);abline(1.3581,0,lty=2) 
  }  
}
par(mfrow=c(2,2))
ph.resid(b)

## stratification example, with 2 randomly allocated strata
## so that results should be similar to previous....
col1$strata <- sample(1:2,nrow(col1),replace=TRUE) 
bs <- gam(cbind(time,strata)~s(age,by=sex)+sex+s(nodes)+perfor+rx+obstruct
          +adhere,family=cox.ph(),data=col1,weights=status)
plot(bs,pages=1,all.terms=TRUE) ## plot effects

## baseline survival plots by strata...

for (i in 1:2) { ## loop over strata
## create index selecting elements of stored hazard info for stratum i...
ind <- which(bs$family$data$tr.strat == i)
if (i==1) plot(bs$family$data$tr[ind],exp(-bs$family$data$h[ind]),type="l",
      ylim=c(0,1),xlab="time",ylab="survival",lwd=2,col=i) else
      lines(bs$family$data$tr[ind],exp(-bs$family$data$h[ind]),lwd=2,col=i)
lines(bs$family$data$tr[ind],exp(-bs$family$data$h[ind] +
      2*bs$family$data$q[ind]^.5),lty=2,col=i) ## upper ci
lines(bs$family$data$tr[ind],exp(-bs$family$data$h[ind] -
      2*bs$family$data$q[ind]^.5),lty=2,col=i) ## lower ci
lines(bs$family$data$tr[ind],exp(-bs$family$data$km[ind]),col=i) ## KM
}


## Simple simulated known truth example...
ph.weibull.sim <- function(eta,gamma=1,h0=.01,t1=100) { 
  lambda <- h0*exp(eta) 
  n <- length(eta)
  U <- runif(n)
  t <- (-log(U)/lambda)^(1/gamma)
  d <- as.numeric(t <= t1)
  t[!d] <- t1
  list(t=t,d=d)
}
n <- 500;set.seed(2)
x0 <- runif(n, 0, 1);x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1);x3 <- runif(n, 0, 1)
f0 <- function(x) 2 * sin(pi * x)
f1 <- function(x) exp(2 * x)
f2 <- function(x) 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
f3 <- function(x) 0*x
f <- f0(x0) + f1(x1)  + f2(x2)
g <- (f-mean(f))/5
surv <- ph.weibull.sim(g)
surv$x0 <- x0;surv$x1 <- x1;surv$x2 <- x2;surv$x3 <- x3

b <- gam(t~s(x0)+s(x1)+s(x2,k=15)+s(x3),family=cox.ph,weights=d,data=surv)

plot(b,pages=1)

## Another one, including a violation of proportional hazards for
## effect of x2...

set.seed(2)
h <- exp((f0(x0)+f1(x1)+f2(x2)-10)/5)
t <- rexp(n,h);d <- as.numeric(t<20)

## first with no violation of PH in the simulation...
b <- gam(t~s(x0)+s(x1)+s(x2)+s(x3),family=cox.ph,weights=d)
plot(b,pages=1)
ph.resid(b) ## fine

## Now violate PH for x2 in the simulation...
ii <- t>1.5
h1 <- exp((f0(x0)+f1(x1)+3*f2(x2)-10)/5)
t[ii] <- 1.5 + rexp(sum(ii),h1[ii]);d <- as.numeric(t<20)

b <- gam(t~s(x0)+s(x1)+s(x2)+s(x3),family=cox.ph,weights=d)
plot(b,pages=1)
ph.resid(b) ## The checking plot picks up the problem in s(x2) 


## conditional logistic regression models are often estimated using the 
## cox proportional hazards partial likelihood with a strata for each
## case-control group. A dummy vector of times is created (all equal). 
## The following compares to 'clogit' for a simple case. Note that
## the gam log likelihood is not exact if there is more than one case
## per stratum, corresponding to clogit's approximate method.
library(survival);library(mgcv)
infert$dumt <- rep(1,nrow(infert))
mg <- gam(cbind(dumt,stratum) ~ spontaneous + induced, data=infert,
          family=cox.ph,weights=case)
ms <- clogit(case ~ spontaneous + induced + strata(stratum), data=infert,
             method="approximate")
summary(mg)$p.table[1:2,]; ms

參考

Hastie and Tibshirani (1990) Generalized Additive Models, Chapman and Hall.

Klein, J.P and Moeschberger, M.L. (2003) Survival Analysis: Techniques for Censored and Truncated Data (2nd ed.) Springer.

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and model selection for general smooth models. Journal of the American Statistical Association 111, 1548-1575 doi:10.1080/01621459.2016.1180986

也可以看看

cox.pht , cnorm