[Case Study]Penalized Semiparametric Additive Hazards Model using R

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Authors:	Anders Gorst-Rasmussen, Thomas H. Scheike
Title:	[download] (2000)Coordinate Descent Methods for the Penalized Semiparametric Additive Hazards Model
Reference:	Vol. 47, Issue 9, Apr 2012Submitted 2011-07-05, Accepted 2012-03-06
Type:	Article
Abstract:	For survival data with a large number of explanatory variables, lasso penalized Cox regression is a popular regularization strategy. However, a penalized Cox model may not always provide the best fit to data and can be difficult to estimate in high dimension because of its intrinsic nonlinearity. The semiparametric additive hazards model is a flexible alternative which is a natural survival analogue of the standard linear regression model. Building on this analogy, we develop a cyclic coordinate descent algorithm for fitting the lasso and elastic net penalized additive hazards model. The algorithm requires no nonlinear optimization steps and offers excellent performance and stability. An implementation is available in the R package ahaz. We demonstrate this implementation in a small timing study and in an application to real data.

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：Case study Parametric penalized semipara Additive difficult survival popular provide always

1. ###########
   
2. # PACKAGES
   
3. ###########
   
4. 
   
5. library("timereg")
   
6. library("lars")
   
7. library("ahaz")
   
8. library("xtable")
   
9. 
   
10. ############################
    
11. # FUNCTION DEFINITIONS ETC.
    
12. ############################
    
13. 
    
14. # Generate data from an AFT model
    
15. # Adapted from the R-code accompanying the paper Simon et al. (2011); see http://www.jstatsoft.org/v39/i05/
    
16. gendata <- function(n,p,rho,snr=3)
    
17.   {
    
18.     # Design matrix
    
19.      if(abs(rho) < 1){
    
20.        beta <- sqrt(rho / (1-rho))
    
21.        X0 <- matrix(rnorm(n * p),ncol=p)
    
22.        z <- rnorm(n)
    
23.        X <-beta * matrix(z,nrow=n,ncol=p,byrow=F) + X0
    
24.      }
    
25.      if(abs(rho)==1){X <- matrix(z,nrow=n,ncol=p,byrow=F)}
    
26. 
    
27.      # True survival times
    
28.      b <- ((-1)^(1:p)) * exp(-2 * ( (1:p) - 1) / 20)
    
29.      f <- X %*% b
    
30.      e <- rnorm(n)
    
31.      k <- sqrt(var(f) / (snr * var(e)))
    
32.      T <- exp(f + k * e)
    
33. 
    
34.      # Censoring times
    
35.      b <- ((-1)^(1:p)) * exp(-2 * ( (1:p) - 1) / 20)
    
36.      f <- X %*% b
    
37.      e <- rnorm(n)
    
38.      k <- sqrt(var(f) / (snr * var(e)))
    
39.      C <- exp(k * e)
    
40.      
    
41.      return(list("X"=X,"surv"=Surv(apply(cbind(C,T),1,min),T<C)))
    
42.   }
    
43. 
    
44. # Local ahaz - get D,d but *not* B
    
45. local.ahaz <- function(surv,X)
    
46.   {
    
47.     a <- .C("aha",
    
48.             "X"=X,
    
49.             "time1"=as.double(rep(0,nrow(X))),
    
50.             "time2"=as.double(surv[,1]),
    
51.             "event"=as.integer(surv[,2]),
    
52.             "weights"=as.double(rep(1,nrow(X))),
    
53.             "n"=nrow(X),
    
54.             "p"=ncol(X),
    
55.             "d"=numeric(ncol(X)),
    
56.             "D"=numeric(ncol(X)^2),
    
57.             "B"=numeric(0),
    
58.             "getB"=as.integer(0),
    
59.             "univariate"=as.integer(0),
    
60.             "usethis"= as.integer(rep(1,ncol(X))),
    
61.             "rightcens"=as.integer(1))
    
62.     D <- matrix(a$D,ncol=ncol(X))
    
63.     D <- D + t(D) - diag(diag(D))
    
64.     d <- a$d
    
65.     return(list("D"=D,"d"=d,"n"=nrow(X)))

复制代码

1. ##############
   
2. # TABLE 1 + 2
   
3. ##############
   
4. 
   
5. # NOTE: Calculations for p=10,000 are memory-intensive!  #
   
6. times <- array(0, c(3,3,5,12))
   
7. n<-c(rep(200,4),rep(500,4),rep(1000,4))
   
8. ptmp<-c(100,500,5000,10000)
   
9. p<-rep(ptmp,3)
   
10. rho <- c(0, 0.25, 0.5, 0.9, 0.95)
    
11. max.steps <- 100
    
12. snr <- 3
    
13. 
    
14. for(k in 1:3){                   # Repetitions
    
15.   for(i in 1:12){                  # Loop over n/p combinations
    
16.     cat(paste("Case",i,"- repetition",k,"\n"))
    
17.     for(j in 1:5){                   # Loop over rho
    
18.       data <- gendata(n[i], p[i], rho[j], snr)
    
19.       times[k,1,j,i] <-  system.time(aha <- local.ahaz(data$surv,data$X))[1]
    
20.       times[k,2,j,i] <-  system.time(lar <- surv.lars(aha$D,aha$d,aha$n,max.steps=max.steps))[1] + times[k,1,j,i]
    
21.       lminf <- min(lar$lambda) / max(lar$lambda)
    
22.       times[k,3,j,i]  <- system.time(ccd <- ahazpen(data$surv,data$X,lambda.minf=lminf,standardize=FALSE))[1]
    
23.     }
    
24.   }
    
25. }
    
26. 
    
27. # LaTeX table (CCD versus LARS+pre-calculation)
    
28. avg.tab1 <- apply(times, c(2,3,4), mean)
    
29. i <- 1
    
30. tab <- c()
    
31. for(nn in 1:length(unique(n))){
    
32.   tmp<-c()
    
33.   for(q in 1:length(ptmp)) {
    
34.     tmp <- cbind(tmp, t(avg.tab1[c(3,2),,i]))
    
35.     i <- i + 1
    
36.   }
    
37.   tab <- rbind(tab,tmp)
    
38. }
    
39. xtable(tab)
    
40. 
    
41. # LaTeX table (pre-calculation only)
    
42. avg.tab2 <- apply(times, c(2,4), mean)[1,]
    
43. i <- 1
    
44. tab <- c()
    
45. for(nn in 1:length(unique(n))){
    
46.   tmp <- c()
    
47.   for(q in 1:length(ptmp)) {
    
48.     tmp <- c(tmp,avg.tab2[i])
    
49.     i <- i + 1
    
50.   }
    
51.   tab <- rbind(tab,tmp)
    
52. }
    
53. xtable(tab)

复制代码

1. ##########
   
2. # TABLE 3
   
3. ##########
   
4. 
   
5. times <- array(0, c(3,5,6))
   
6. n <- c(200,200,200,40000,100000,250000)
   
7. p <- c(40000,100000,250000,200,200,200)
   
8. snr <- 3
   
9. rho <- c(0, 0.25, 0.5, 0.9, 0.95)
   
10. for(k in 1:1){                       # Repetitions
    
11.   for(i in 1:6){                       # Loop over n/p combinations
    
12.     cat(paste("Case",i,"- repetition",k,"\n"))
    
13.     for(j in 1:5){                       # Loop over rho
    
14.       
    
15.       data <- gendata(n[i], p[i], rho[j], snr)
    
16.       ccd1 <- ahazpen(data$surv,data$X,dfmax=120,pmax=1000,standardize=FALSE) # Large pmax: many variables considered for rho=.95
    
17.       lminf <- ccd1$lambda[which.min(abs(ccd1$df - 100))] / ccd1$lambda[1]
    
18.       times[k,j,i]  <- system.time(ccd2 <- ahazpen(data$surv,data$X,lambda.minf=lminf,standardize=FALSE))[1]
    
19.     }
    
20.   }
    
21. }

复制代码

1. #################
   
2. # SORLIE EXAMPLE
   
3. #################
   
4. 
   
5. # Load data
   
6. data("sorlie")
   
7. set.seed(10101)
   
8. surv <- Surv(sorlie$time + runif(nrow(sorlie)) * 1e-2, sorlie$status)
   
9. Z <- sorlie[,3:ncol(sorlie)]
   
10. p <- ncol(Z)
    
11. pw.comb <- combn(1:p,2)
    
12. X <- cbind(Z, Z[, pw.comb[1,]] * Z[, pw.comb[2,]])
    
13. 
    
14. # Initial fit
    
15. fit.init <- ahazpen(surv, X, dfmax = 50)
    
16. 
    
17. # Expand lambda grid
    
18. l <- range(fit.init$lambda)
    
19. fit <- ahazpen(surv, X, lambda.minf = l[1] / l[2])
    
20. plot(fit)
    
21. 
    
22. # Cross-validation
    
23. set.seed(10101)
    
24. fit.tune <- tune.ahazpen(surv, X, lambda.minf = l[1] / l[2], tune = "cv")
    
25. plot(fit.tune)
    
26. 
    
27. # Final estimate
    
28. beta <- coef(fit.tune)
    
29. which(as.numeric(beta)!=0)

复制代码