# R/convertCDA.R In smoothSurv: Survival Regression with Smoothed Error Distribution

```###########################################
#### AUTHOR:    Arnost Komarek         ####
####            (2003)                 ####
####                                   ####
#### FILE:      convertCDA.R           ####
####                                   ####
#### FUNCTIONS: c2a                    ####
####            a2c                    ####
####            find.c                 ####
####            give.c                 ####
####            derivative.cc3         ####
###########################################

### ========================================
### c2a: Compute a coefficients from c's
### ========================================
c2a <- function(ccoef, which.zero = which.max(ccoef), toler = 1e-6){
ccoef[ccoef < toler] <- toler
c.zero <- ccoef[which.zero]

acoef <- log(ccoef/c.zero)
return(acoef)
}

### ========================================
### a2c: Function to compute c's from a's
### =========================================
a2c <- function(acoef){
ccoef <- exp(acoef)
sum.exp.a <- sum(ccoef)
if (is.nan(sum.exp.a)) return(NULL)
ccoef <- ccoef/sum.exp.a

return(ccoef)
}

### =======================================================================================
### derivative.expAD: Function to compute derivatives of non-zero exp(a)'s w.r.t. exp(d)'s
### =======================================================================================
##
##  * there are g - 1 non-zero a's
##  * there are g - 3 d's
##
##  * g - 3 a's are equal to d's
##  * 2 a's are a function of the rest
##  * 1 a is equal to zero
##
## INPUT: knots ....... vector of knots
##        sdspline .... standard deviation of a basis G-spline
##        last.three ... vector with indeces of the three a's which are to be computed from the rest
##             a[last.three] = 0
##             a[last.three] = first function(d's)
##             a[last.three] = second function(d's)
##        all ......... do I want the full  matrix or only two columns w.r.t.
##                      to the two a's

## OUTPUT: Omega.....
##               if all == TRUE, matrix (g - 2) x g (there is one zero column)
##                  all == FALSE, matrix (g - 2) x 2
##                         (the first row is always an intercept)
derivative.expAD <- function(knots, sdspline, last.three, all = TRUE){

g <- length(knots)
if (g < 4){
stop("Too short input 'knots' vector ")
}

if (length(last.three) != 3){
stop("Incorrect 'last.three' parameter ")
}

if (sum(last.three %in% (1:g)) != 3){
stop("Incorrect 'last.three' parameter ")
}

if (last.three == last.three || last.three == last.three || last.three == last.three){
stop("Incorrect 'last.three' parameter ")
}

if (sdspline >= 1 || sdspline <= 0){
stop("Incorrect 'sdspline' parameter ")
}

which.zero <- last.three
l1 <- last.three
l2 <- last.three
s02 <- sdspline * sdspline

if (abs(knots[l1]) < 1e-4){
stop("Zero reference knot in derivative.expAD ")
}

kn2.kn1 <- knots[l2] - knots[l1]
jsmm <- 1 - s02 + knots[l1]*knots[l2]

## Knots with removed the two ones corresponding to the two special a's
knotsMin2 <- knots[-c(l1, l2)]

## Index of the zero a in the shorter (by 2) a's sequence
which.zero2 <- ifelse(which.zero < min(l1, l2),
which.zero,
ifelse(which.zero < max(l1, l2),
which.zero - 1,
which.zero - 2))

## Compute the two columns of the resulting matrix corresponding to the two a's
vec2b <- -(1/kn2.kn1) * (knotsMin2 - knots[l1])
vec2c <- (1/jsmm) * (1 - s02 + knots[l1] * knotsMin2)
vec2 <- (vec2b & vec2c)

vec1a <- knots[l2]/knots[l1]
vec1d <- (1/knots[l1]) * knotsMin2
vec1 <- -vec1a * vec2 - vec1d

int1 <- vec1[which.zero2]
int2 <- vec2[which.zero2]

slope1 <- vec1[-which.zero2]
slope2 <- vec2[-which.zero2]

slope <- cbind(slope1, slope2)
intercept <- c(int1, int2)

Omega <- rbind(intercept, slope)

if (all){
leftMat <- rbind(matrix(0, nrow = 1, ncol = g - 3), diag(g - 3))
zeroCol <- matrix(c(1, rep(0, g - 3)), nrow = g - 2, ncol = 1)
k <- 1
kk <- 1
for (j in 1:g){
if (j == which.zero){
Omega <- cbind(Omega, zeroCol)
}
else{
if (j == l1 || j == l2){
Omega <- cbind(Omega, Omega[,k])
k <- k + 1
}
else{
Omega <- cbind(Omega, leftMat[,kk])
kk <- kk + 1
}
}
}
## the first two columns are now obscure -> remove them
Omega <- Omega[,-c(1,2)]
}

return(Omega)
}

### ===================================================================
### find.c: Find mixture proportions that approximate
###         given distribution (dist) by a G-spline mixture with knots
###         and standard deviation sdspline
### ===================================================================
## RETURNS: vector with c coeff.
##          or NULL if problems to find them
find.c <- function(knots, sdspline, dist = "dnorm"){
nsplines <- length(knots)
in.knots <- list(x=as.numeric(knots))
right.side <- do.call(dist, in.knots)
right.side[abs(right.side) < 1e-10] <- 0
mus <- matrix(rep(knots, rep(nsplines, nsplines)), nrow = nsplines)
knotsmat <- matrix(rep(knots, nsplines), nrow = nsplines)
Cmat <- dnorm(knotsmat, mean=mus, sd=sdspline)
Cmat <- qr(Cmat, tol = 1e-07)
if (Cmat\$rank == ncol(Cmat\$qr)){
ccoef <- solve(Cmat, right.side)
tempmin <- min(ccoef[ccoef > 0])
ccoef[ccoef <= 0] <- tempmin           ## this is not theoretically possible but numerically it is
}
else
ccoef <- NULL

return(ccoef)
}

### =================================================================
### give.c: Give a vector of all c's satisfying the three constrains
###         if remaining (g-3) c's are given.
### =================================================================
## INPUT: knots ......... knots
##        sdspline ...... standard deviation of a G-spline
##        c.rest ..... remaining g-3 c coefficients
##        last.three... indeces of the three c coefficients which are a function of the rest ones
## RETURN: a vector of all c's
give.c <- function(knots, sdspline, last.three, c.rest)
{
g <- length(knots)

if (length(c.rest) != g - 3)
stop("Incorrect dimension of the input vector ")

if (length(last.three) != 3) stop("Incorrect 'last.three' parameter ")
if (sum(last.three %in% 1:g) != 3) stop("Incorrect 'last.three' parameter ")
if (length(unique(last.three)) != 3) stop("Incorrect 'last.three' parameter ")

## Matrix to compute last c's from the first g - 3 ones
Omega <- derivative.cc3(knots, sdspline, last.three, all = TRUE)

## compute all c's
Omega0 <- Omega[1,]
Omega1 <- matrix(Omega[2:(g - 2),], nrow = g - 3)
c.all <- (t(Omega1) %*% c.rest) + Omega0

return(as.numeric(c.all))
}

### ==========================================================
### derivative.cc3: Derivatives of all c's w.r.t. (g - 3) c's
### ==========================================================
## INPUT: knots ....... knots
##        sdspline .... standard deviation of a G-spline
##        last.three... indeces of the three c coefficients which are a function of the rest ones
##        all ....... if TRUE, matrix to compute all c's from the first three ones is returned
##                    if FALSE, matrix to compute last three c's from the first three ones is returned
## RETURN: Matrix where the first row is an intercept
##         and remaining (g-3) rows is dc/da
derivative.cc3 <- function(knots, sdspline, last.three, all = TRUE){

g <- length(knots)
if (g < 4) stop("Too short input 'knots' vector ")

if (length(last.three) != 3) stop("Incorrect 'last.three' parameter ")
if (sum(last.three %in% 1:g) != 3) stop("Incorrect 'last.three' parameter ")
if (length(unique(last.three)) != 3) stop("Incorrect 'last.three' parameter ")

last.three <- last.three[order(last.three)]

s02 = sdspline * sdspline

## Indices of the three c coefficients
l0 <- last.three
l1 <- last.three
l2 <- last.three

## Compute first the matrix used to compute last three c's from the first g - 3 ones
Omega = matrix(0, nrow = g - 2, ncol = 3)

## 1st row (intercept)
Omega[1, 1] = (1 - s02 + knots[l2]*knots[l1])/((knots[l2] - knots[l0])*(knots[l1] -knots[l0]))
Omega[1, 2] = -(1 - s02 + knots[l2]*knots[l0])/((knots[l2] - knots[l1])*(knots[l1] -knots[l0]))
Omega[1, 3] = 1 - Omega[1, 1] - Omega[1, 2]

## the rest  (loop over rows)
i <- 2
for (j in 1:g){
if (j == l0 || j == l1 || j == l2) next;
Omega[i, 1] = -((knots[l2] - knots[j])*(knots[l1] - knots[j]))/((knots[l2] - knots[l0])*(knots[l1] -knots[l0]));
Omega[i, 2] = ((knots[l2] - knots[j])*(knots[l0] - knots[j]))/((knots[l2] - knots[l1])*(knots[l1] -knots[l0]));
Omega[i, 3] = -(1 + Omega[i, 1] + Omega[i, 2])
i <- i + 1
}

## Add the components to compute all c's from the first g - 3 ones
if (all){
leftmat <- rbind(rep(0, g - 3), diag(g - 3))
OmegaAll <- NULL
k <- 1; kk <- 1
for (j in 1:g){
if (j == l0 || j == l1 || j == l2){
OmegaAll <- cbind(OmegaAll, Omega[,k])
k <- k + 1
}
else{
OmegaAll <- cbind(OmegaAll, leftmat[, kk])
kk <- kk + 1
}
}
Omega <- OmegaAll
}

return(Omega)
}
```

## Try the smoothSurv package in your browser

Any scripts or data that you put into this service are public.

smoothSurv documentation built on July 8, 2020, 5:46 p.m.