Nothing
# $Id: coxpenal.df.S 11166 2008-11-24 22:10:34Z therneau $
#
# degrees of freedom computation, based on Bob Gray's paper
#
# hmat = right hand slice of cholesky of H
# hinv = right hand slice of cholesky of H-inverse
# fdiag= diagonal of D-inverse
# assign.list: terms information
# ptype= 1 or 3 if a sparse term exists, 2 or 3 if a non-sparse exists
# nvar = # of non-sparse terms
# pen1 = the penalty matrix (diagonal thereof) for the sparse terms
# pen2 = the penalty matrix for the non-sparse terms
# sparse = indicates which term is the sparse one
coxpenal.df <- function(hmat, hinv, fdiag, assign.list, ptype, nvar,
pen1, pen2, sparse) {
if (ptype ==1 & nvar==0) { #only sparse terms
hdiag <- 1/fdiag
list(fvar2=(hdiag-pen1)*fdiag^2, df=sum((hdiag-pen1)*fdiag),
fvar = fdiag, trH=sum(fdiag))
}
else if (ptype==2) { # only dense ones
hmat.full <- t(hmat) %*% (ifelse(fdiag==0, 0,1/fdiag)* hmat)
hinv.full <- hinv %*% (fdiag* t(hinv))
if (length(pen2)==length(hmat.full)) imat <- hmat.full - pen2
else imat <- hmat.full - diag(pen2)
var <- hinv.full %*% imat %*% hinv.full
if (length(assign.list)==1)
list(var2=var, df=sum(imat * hinv.full),
trH=sum(diag(hinv.full)), var=hinv.full)
else {
df <- trH <- NULL
d2 <- diag(hinv.full)
for (i in assign.list) {
temp <- coxph.wtest(hinv.full[i,i], var[i,i])$solve
if (is.matrix(temp)) df <- c(df, sum(diag(temp)))
else df <- c(df, sum(temp))
trH<- c(trH, sum(d2[i]))
}
list(var2=var, df=df, trH=trH, var = hinv.full)
}
}
else {
# sparse terms + other vars
nf <- length(fdiag) - nvar
nr1 <- 1:nf
nr2 <- (nf+1):(nf+nvar)
d1 <- fdiag[nr1]
d2 <- fdiag[nr2]
temp <- t(hinv[nr1,])
temp2<- t(hinv[nr2,,drop=FALSE])
A.diag <- d1 + c(rep(1,nvar) %*% (temp^2*d2))
B <- hinv[nr1,] %*% (d2 * temp2)
C <- hinv[nr2,] %*% (d2 * temp2) #see notation in paper
var2 <- C - t(B) %*% (pen1 * B)
if (ptype==3) {
#additional work when we have penalties on both the sparse term
# and on non-sparse terms
hmat.22 <- t(hmat) %*%(ifelse(fdiag==0, 0,1/fdiag)* hmat)
temp <- C - coxph.wtest(hmat.22, diag(nvar))$solve
if (nvar==1) {
var2 <- var2 - C*pen2*C # C will be 1 by 1
temp2 <- c(temp*pen2)
}
else if (length(pen2) == nvar) {
var2 <- var2 - C %*% (pen2 * C) #diagonal penalty
temp2 <- sum(diag(temp) * pen2)
}
else {
var2 <- var2 - C %*% matrix(pen2,nvar) %*% C
temp2 <- sum(diag(temp * pen2))
}
}
else temp2 <- 0 #temp2 contains trace[B'A^{-1}B P2], this line: P2=0
df <- trH <- NULL
cdiag <- diag(C)
for (i in 1:length(assign.list)) {
if (sparse==i){
df <- c(df, nf - (sum(A.diag * pen1) + temp2))
trH <- c(trH, sum(A.diag))
}
else {
j <- assign.list[[i]]
temp <- coxph.wtest(C[j,j], var2[j,j])$solve
if (is.matrix(temp)) df <- c(df, sum(diag(temp)))
else df <- c(df, sum(temp))
trH <- c(trH, sum(cdiag[j]))
}
}
list(var=C, df=df, trH=trH, fvar=A.diag, var2=var2)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.