API tools faq
paste
Advertisement
mgordon

coxvc with a bug-fix for the plot & added feature

mgordon
Dec 1st, 2011
1,056
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 11.13 KB | None | 0 0
raw download clone embed print report
  1. .packageName <- "coxvc"
  2. calc.h0 <- function(obj) {
  3.     if (is(obj,"coxph")) {
  4.         ha <- coxph.detail(obj)$haz
  5.         ha*exp(-sum(obj$means*obj$coef))}
  6.     else {
  7.         ha <- obj$haz
  8.         ha * exp(-sum(obj$means * obj$theta))}}
  9. coxvc <- function (formula = formula(data), Ft, rank, iter.max = 30, data = sys.parent())
  10. {
  11.     call <- match.call()
  12.     m <- match.call(expand = FALSE)
  13.     temp <- c("", "formula", "data")
  14.     m <- m[match(temp, names(m), nomatch = 0)]
  15.     Terms <- terms(formula, data = data)
  16.     m$formula <- Terms
  17.     m[[1]] <- as.name("model.frame")
  18.     m <- eval(m, sys.parent())
  19.     if (NROW(m) == 0)
  20.         stop("No (non-missing) observations")
  21.     Y <- model.extract(m, "response")
  22.     if (!inherits(Y, "Surv"))
  23.         stop("Response must be a survival object")
  24.     X <- model.matrix(Terms, m)[, -1, drop = FALSE]
  25.     time <- Y[, 1]
  26.     ord <- order(time)
  27.     Ftime <- Ft[ord, ]
  28.     n <- nrow(X)
  29.     d <- Y[, 2]
  30.     p <- ncol(X)
  31.     q <- ncol(Ftime)
  32.     r <- rank
  33.     if (sum(Ft[, 1]) != n)
  34.         stop("The matrix of time functions should include a constant in the first column")
  35.     nm <- attr(Terms, "term.labels")
  36.     tmp <- paste("f", 1:(q - 1), "(t)", sep = "")
  37.    
  38.     X <- apply(X, 2, function(x) {
  39.                 x - mean(x)
  40.             })
  41.    
  42.     means<-  matrix(apply(matrix(apply(X,2,function(x)x*Ft),ncol=p*q),2,mean),ncol=q,byrow=TRUE)
  43.    
  44.     if (r > p | r > q)
  45.         stop("Rank r cannot exceed the number of covariates or time functions",
  46.                 "\n", "Full rank model for this data is: R=", min(p,
  47.                         q))
  48.     eventtimes <- (1:n)[d == 1]
  49.     k <- length(eventtimes)
  50.     Ftime <- Ftime[eventtimes, ]
  51.     if (r == min(p, q)) {
  52.         theta <- matrix(rep(0, p * q), ncol = q)
  53.         fvc <- full.fit(eventtimes, X, Ftime, theta, iter.max,
  54.                 n, p, q)
  55.         nm <- c(nm, paste(rep(nm, q - 1), rep(tmp, each = p),
  56.                         sep = ":"))
  57.         fit <- list(call = call, theta = fvc$th, ster = fvc$se, rank=r,
  58.                 logL = fvc$lik, iter = fvc$iter, coef.names = nm, hazard = fvc$hi, means=means,time=Y[,1],death=Y[, 2],Ftime=Ft)
  59.         class(fit) <- "coxvc"
  60.         fit
  61.     }
  62.     else {
  63.         b <- matrix(rep(0, r * p), ncol = r)
  64.         gama <- matrix(runif(q * r)/10, ncol = r)
  65.         if (r == 1) {
  66.             gama[1, 1] <- 1
  67.         }
  68.         else {
  69.             gama[, 1] <- 1
  70.         }
  71.         frr <- rr.fit(eventtimes, X, Ftime, b, gama, iter.max,
  72.                 n, p, q, r)
  73.         nm <- c(nm, paste(rep(nm, q - 1), rep(tmp, each = p),
  74.                         sep = ":"))
  75.         fit <- list(call = call, b = frr$b, gama = frr$gama,
  76.                 theta = frr$th, ster = frr$se, logL = frr$lik, iter = frr$iter,
  77.                 rank = r, coef.names = nm, hazard=frr$hi,means=means,time=Y[,1],death=Y[, 2],Ftime=Ft)
  78.         class(fit) <- "coxrr"
  79.         fit
  80.     }
  81. }
  82. expand.haz <- function(haz,status,fun=c("baseline","cumulative")){
  83.     c <- status
  84.     c[status==0]  <- 1
  85.     c[status==1] <- 0
  86.     h2 <- NULL
  87.     it <- 1
  88.     k <- 1
  89.     while(it < (length(status)+1)){
  90.         if (c[it]==0){h2[it] <- haz[it-k+1]}
  91.         if (c[it]==1 && fun=="cumulative"){
  92.             h2[it] <-  h2[it-1]
  93.             k <- k+1}
  94.         if (c[it]==1 && fun=="baseline"){
  95.             h2[it] <-  0
  96.             k <- k+1}
  97.         it <- it+1}
  98.     h2}      
  99. full.fit <-
  100.         function (eventtimes, X, Ftime, theta, iter.max, n, p, q)
  101. {
  102.     delta <- 5
  103.     iter <- 0
  104.     k <- length(eventtimes)
  105.     while (delta > 1e-06 && iter < iter.max) {
  106.         th <- theta
  107.         iter <- iter + 1
  108.         tmp <- sapply(eventtimes, sumevents, X, Ftime, theta,
  109.                 n, p, eventtimes)
  110.         hi <- unlist(tmp[1, ])
  111.         lik <- sum(log(hi)) + sum(X[eventtimes, ] * t(theta %*%
  112.                                 t(Ftime)))
  113.         scores <- matrix(unlist(tmp[2, ]), ncol = k) %*% rep.int(1,
  114.                 k) + matrix(t(X[eventtimes, ]) %*% Ftime, ncol = 1)
  115.         hess <- -(matrix(matrix(unlist(tmp[3, ]), ncol = k) %*%
  116.                                     rep(1, k), ncol = p * q, byrow = FALSE))
  117.         theta <- matrix(matrix(theta, ncol = 1) - solve(hess,
  118.                         matrix(t(scores), ncol = 1)), ncol = q)
  119.         delta <- max(abs(theta - th))
  120.         if (iter > iter.max)
  121.             stop("Likelihood did not converge after ", iter,
  122.                     " iterations")
  123.     }
  124.     sterr <- round(matrix(sqrt(abs(diag(solve(hess)))), nrow = p),
  125.             3)
  126.     outlist <- list(th = round(theta, 3), se = sterr, lik = lik,
  127.             iter = iter, hi = hi)
  128.     outlist
  129. }
  130.  
  131. plotcoxvc <- function(object, fun=c("survival","effects"),
  132.         ylab="",xlab="",lwd=1,col=1,lty=1,ylim=NULL,
  133.         legend_y=NULL, legend_x=NULL, terms=NULL, ...){
  134.     if(!inherits(object, "coxvc")  && !inherits(object,"coxrr"))
  135.         stop("First arg must be the result of coxvc")
  136.  
  137.     # It the time variables are scrambled generating a very irratic line
  138.     getMatrixForPlot <- function (time, death, y_axis){
  139.         x_axis <- time[death==1]
  140.         base_2d_matrix <- cbind(x_axis, y_axis)
  141.         time_ordered_matrix <- base_2d_matrix[order(base_2d_matrix[ ,1]), ]
  142.         return(time_ordered_matrix)
  143.     }
  144.    
  145.     if (fun=="survival"){
  146.         haz <- object$haz * exp(-sum(object$means * object$theta))
  147.         surv <- exp(-cumsum(haz*(exp( sum(object$means * object$theta)))))    
  148.         if (length(ylim)!=2){
  149.             ylim <- c(0,1)
  150.         }
  151.         plot(getMatrixForPlot(object$time, object$death, surv),
  152.                 ylim=ylim, type='s',
  153.                 ylab=ylab, xlab=xlab,
  154.                 lwd=lwd, col=col, lty=lty)
  155.     }
  156.     if (fun=="effects"){
  157.         th <- object$theta
  158.         p <- nrow(th)
  159.         nms <- object$coef.names[1:nrow(th)]
  160.         Ft <- object$Ftime[object$death==1,]
  161.         ef <- th%*%t(Ft)
  162.         if (length(terms) > 0){
  163.             ef <- ef[terms, ]
  164.             p <- length(terms)
  165.             nms <- nms[terms]
  166.         }
  167.  
  168.         if (length(ylim)!=2){
  169.             ylim <- c(floor(min(ef)),ceiling(max(ef)))
  170.         }
  171.         plot(getMatrixForPlot(object$time, object$death, ef[1,]),
  172.                 ylim=ylim, type='l',
  173.                 lwd=lwd, lty=lty, col=col,
  174.                 ylab=ylab, xlab=xlab)    
  175.         i <- 1
  176.         while (i < p) {
  177.             i <- i+1
  178.             lines(getMatrixForPlot(object$time, object$death, ef[i,]),lty=i,col=i,lwd=lwd)
  179.         }
  180.        
  181.         if (length(legend_y) == 0){
  182.             legend_y <- ceiling(max(ylim))
  183.         }
  184.        
  185.         if (length(legend_x) == 0){
  186.             legend_x <- min(object$time[object$death==1])
  187.         }
  188.        
  189.         legend(x=legend_x, y=legend_y,
  190.                 legend=nms, lty=1:p, col=1:p)
  191.     }
  192.    
  193. }
  194.  
  195. print.coxrr <- function(x,...) {
  196.     cat("call:"); cat("\n")
  197.     print(x$call);cat("\n")
  198.     coef <- as.vector(x$theta)
  199.     se <- as.vector(x$ster)
  200.     tmp <- cbind(coef,exp(coef),se)
  201.     dimnames(tmp) <- list(x$coef.names,c("coef","exp(coef)","se(coef)"))
  202.     prmatrix(tmp)
  203.     cat("\n")
  204.     cat("log-likelihood= ",x$logL,", ","Rank=",x$rank,"\n")
  205.     cat("algorithm converged in", x$iter ,"iterations","\n")
  206.     cat("\n")
  207.     cat("Beta :"); cat("\n"); prmatrix(round(x$b,4))   ;cat("\n")
  208.     cat("Gamma:"); cat("\n"); prmatrix(round(x$gama,4));cat("\n") }
  209.  
  210. print.coxvc <- function (x, digits =max(options()$digits-4,3),...)
  211. {
  212.     cat("call:")
  213.     cat("\n")
  214.     print(x$call)
  215.     cat("\n")
  216.     coef <- as.vector(x$theta)
  217.     se <- as.vector(x$ster)
  218.     tmp <- cbind( round(coef,4), round(exp(coef),3), round(se,3), round(coef/se,4),   signif(1 - pchisq((coef/se)^2, 1), digits-1 ))
  219.     dimnames(tmp) <- list(x$coef.names, c("coef", "exp(coef)",
  220.                     "se(coef)", "z", "p"))
  221.     prmatrix(tmp)
  222.     cat("\n")
  223.     cat("log-likelihood= ", x$logL)
  224.     cat("\n")
  225.     cat("\n")
  226.     cat("algorithm converged in", x$iter, "iterations")
  227. }      
  228. rr.fit <-
  229.         function (eventtimes, X, Ftime, b, gama, iter.max, n, p, q, r)
  230. {
  231.     delta <- 5
  232.     iter <- 0
  233.     k <- length(eventtimes)
  234.     while (delta > 1e-06 && iter < iter.max) {
  235.         theta <- b %*% t(gama)
  236.         iter <- iter + 1
  237.         tmp <- sapply(eventtimes, sumevents, X, Ftime, theta,
  238.                 n, p, eventtimes)
  239.         hi <- unlist(tmp[1, ])
  240.         lik <- sum(log(hi)) + sum(X[eventtimes, ] * t(theta %*%
  241.                                 t(Ftime)))
  242.         scores <- matrix(unlist(tmp[2, ]), ncol = k) %*% rep.int(1,
  243.                 k) + matrix(t(X[eventtimes, ]) %*% Ftime, ncol = 1)
  244.         hess <- -(matrix(matrix(unlist(tmp[3, ]), ncol = k) %*%
  245.                                     rep(1, k), ncol = p * q, byrow = FALSE))
  246.         vecu <- as.vector(scores)
  247.         parb <- as.vector(b)
  248.         db <- kronecker(gama, diag(p))
  249.         d1 <- t(db) %*% vecu
  250.         d2 <- -t(db) %*% hess %*% db
  251.         parb <- parb + solve(d2) %*% d1
  252.         b <- matrix(parb, ncol = r)
  253.         theta <- b %*% t(gama)
  254.         tmp <- sapply(eventtimes, sumevents, X, Ftime, theta,
  255.                 n, p, eventtimes)
  256.         hi <- unlist(tmp[1, ])
  257.         lik2 <- sum(log(hi)) + sum(X[eventtimes, ] * t(theta %*%
  258.                                 t(Ftime)))
  259.         scores <- matrix(unlist(tmp[2, ]), ncol = k) %*% rep.int(1,
  260.                 k) + matrix(t(X[eventtimes, ]) %*% Ftime, ncol = 1)
  261.         hess <- -(matrix(matrix(unlist(tmp[3, ]), ncol = k) %*%
  262.                                     rep(1, k), ncol = p * q, byrow = FALSE))
  263.         vecu <- scores
  264.         parg <- as.vector(t(gama))
  265.         dg <- kronecker(diag(q), b)
  266.         d1 <- t(dg) %*% vecu
  267.         d2 <- -t(dg) %*% hess %*% dg
  268.         parg <- parg + solve(d2) %*% d1
  269.         gama <- matrix(parg, ncol = r, byrow = TRUE)
  270.         delta <- lik2 - lik
  271.     }
  272.     parb <- as.vector(b)
  273.     parg <- as.vector(t(gama))
  274.     parbg <- c(parb, parg)
  275.     db <- kronecker(gama, diag(p))
  276.     dg <- kronecker(diag(q), b)
  277.     dbg <- cbind(db, dg)
  278.     d2 <- t(dbg) %*% hess %*% dbg
  279.     covtheta <- dbg %*% ginv(d2, tol = 10^-15) %*% t(dbg)
  280.     sterr <- sqrt(abs(diag(covtheta)))
  281.     sterr <- round(matrix(sterr, nrow = p), 3)
  282.     outlist <- list(b = b, gama = gama, th = round(theta, 3),
  283.             se = sterr, lik = lik2, iter = iter, hi =hi)
  284. }
  285.  
  286.  
  287.  
  288. sumevents <- function(i,X,Ftime,theta,n,p,eventtimes){
  289. #function to compute risk sets, weighted means,part of scores and information
  290.     xj <- matrix(X[i:n,],nrow=(n-i+1))
  291.     Fi <- Ftime[eventtimes==i,]
  292.     rr <-exp(xj %*% theta %*% Fi) #relative risk
  293.     hi <- 1/sum(rr)  #baseline hazard
  294.     xmean <- hi*t(xj)%*%rr        
  295.     xdif <- xj -matrix(rep(xmean,n-i+1),ncol=p, byrow=TRUE)
  296.     XX <- t(xdif)%*% (xdif*matrix(rep(rr,p),ncol=p,byrow=FALSE))
  297.     FF <- Fi %*% t(Fi)
  298.     list (hi,                      # hazard hi
  299.             -xmean%*% t(Fi),            # part of term i of score
  300.             hi*kronecker(FF,XX))}
  301.  
  302. summary.coxrr <- function(object,...) {
  303.     cat("call:"); cat("\n")
  304.     print(object$call);cat("\n")
  305.     cat("Beta :"); cat("\n"); prmatrix(round(object$b,4))   ;cat("\n")
  306.     cat("Gamma:"); cat("\n"); prmatrix(round(object$gama,4));cat("\n")
  307. }
  308.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!