swh:1:snp:ff0951ca787d0b7f47dc2335f47fed43820a6324
Tip revision: 8aca325bb4949c0a75f3e27fabfa944a7a8a503d authored by Venkatraman E. Seshan on 18 January 2012, 00:00:00 UTC
version 0.9.9
version 0.9.9
Tip revision: 8aca325
coxphCPE.R
##
## coded by E. S. Venkatraman based on the R code by Gonen 02/08/2007
##
coxphCPE <- function(phfit) {
if (class(phfit) != "coxph") stop("phfit shoud be coxph class object")
n <- phfit$n
betahat <- phfit$coefficients
p <- length(phfit$coefficients)
vbetahat <- phfit$var
xbeta <- phfit$linear.predictors
if (getRversion() <= "2.9.1") {
xMat <- as.matrix(model.matrix(phfit)[,-1])
} else {
xMat <- as.matrix(model.matrix(phfit))
}
bw <- 0.5*sd(xbeta)*(n^(-1/3))
zzz <- .Fortran("cpesub",
as.integer(n),
as.integer(p),
as.double(xMat),
as.double(xbeta),
as.double(bw),
CPE=double(1),
CPEsmooth=double(1),
varDeriv=double(p),
uRowSum=double(n),
uSSQ=double(1),
PACKAGE="clinfun")
CPE <- 2*zzz$CPE/(n*(n-1))
CPEsmooth <- 2*zzz$CPEsmooth/(n*(n-1))
varTerm1 <- 4*(sum((zzz$uRowSum+rep(0.5,n)-n*CPEsmooth)^2) - (zzz$uSSQ + n/4 - n*CPEsmooth - n*(n-2)*CPEsmooth^2))/(n*(n-1))^2
varDeriv <- 2*zzz$varDeriv/(n*(n-1))
varTerm2 <- t(varDeriv)%*%vbetahat%*%varDeriv
varCPE <- varTerm1 + varTerm2
out <- c(CPE, CPEsmooth, sqrt(varCPE))
names(out) <- c("CPE", "smooth.CPE", "se.CPE")
out
}