# R/pade.R In survival: Survival Analysis

```# This is a trimmed copy of expm:::expm2.R, optimized for transition matrices
# The row sums of A are 0, with the diagonal the only negative element, and
#    transitions away from the diagonal are normally small.
#  I don't need balancing, don't need Matrix objects
#  The matrices are usually very nice.
#  If I need derivatives, then deriv is an (n,n, k) array, the first
#  slice contains dA/(d beta1), elementwise,  and etc for the second, third etc.
#
n <- nrow(A)
nA <- max(colSums(abs(A)))
I <- diag(n)
if (!missing(deriv)) {
if (!is.array(deriv) || length(dim(deriv))!=3 ||
any(dim(deriv)[1:2] != n)) stop ("invalid derivative array")
nderiv <- dim(deriv)[3]
dmat <- array(0., dim=c(n, n, nderiv))
}
else nderiv <- 0

## If the norm is small enough, use the Pade Approximation (PA) directly
if (nA <= 2.1) {
t <- c(0.015, 0.25, 0.95, 2.1)
## the minimal m for the PA :
l <- which.max(nA <= t)
nterm <- l+1  # for reporting: number of terms used

## Calculate PA
C <- rbind(c(120,60,12,1,0,0,0,0,0,0),
c(30240,15120,3360,420,30,1,0,0,0,0),
c(17297280,8648640,1995840,277200,25200,1512,56,1,0,0),
c(17643225600,8821612800,2075673600,302702400,30270240,
2162160,110880,3960,90,1))
A2 <- A %*% A
P <- I
U <- C[l,2]*I
V <- C[l,1]*I

for (k in 1:l) {
P <- P %*% A2
U <- U + C[l,(2*k)+2]*P
V <- V + C[l,(2*k)+1]*P
}
U <- A %*% U
X <- solve(V-U,V+U)

if (nderiv > 0) {
for (dk in 1:nderiv) {
P <- I
U <- C[l,2]*I
V <- C[l,1]*I
dA <- deriv[,,dk]
dA2 <- A %*% dA + dA %*% A
dU <- dV <- dP <-  0*I
for (k in 1:l) {
dP <- dP %*% A2 + P %*% dA2
P <- P %*% A2
U <- U + C[l,(2*k)+2]*P
dU <- dU + C[l,(2*k)+2]* dP
V <- V + C[l,(2*k)+1]*P
dV <- dV + C[l,(2*k)+1]*dP
}
dU <- dA %*%U + A %*% dU
U <- A %*% U
dmat[,,dk] <- -solve(V-U, dV-dU)%*% X + solve(V-U, dV + dU)
}
}
}
else {
# This section should be rarer than rare
s <- log2(nA/5.4)
B <- A
## Scaling
if (s > 0) {
s <- ceiling(s)
B <- B/(2^s)
if (nderiv>0) deriv <- deriv/(2^s)
}
nterm <- max(0, ceiling(s)) + 6

## Calculate PA
c. <- c(64764752532480000,32382376266240000,7771770303897600,
1187353796428800, 129060195264000,10559470521600, 670442572800,
33522128640, 1323241920, 40840800,960960,16380, 182,1)
B2 <- B %*% B
B4 <- B2 %*% B2
B6 <- B2 %*% B4

U <- B %*% (B6 %*% (c.[14]*B6 + c.[12]*B4 + c.[10]*B2) +
c.[8]*B6 + c.[6]*B4 + c.[4]*B2 + c.[2]*I)
V <- B6 %*% (c.[13]*B6 + c.[11]*B4 + c.[9]*B2) +
c.[7]*B6 + c.[5]*B4 + c.[3]*B2 + c.[1]*I

X <- solve(V-U,V+U)

if (nderiv > 0) {
for (dk in 1:nderiv) {
dB <- deriv[,,dk]
dB2 <- B %*% dB + dB %*% B
dB4 <- B2%*% dB2 + dB2 %*% B2
dB6 <- B2 %*%dB4 + dB2 %*% B4
dU <- dB %*% (B6 %*% (c.[14]*B6 + c.[12]*B4 + c.[10]*B2) +
c.[8]*B6 + c.[6]*B4 + c.[4]*B2 + c.[2]*I) +
B %*% (dB6 %*% (c.[14]*B6 + c.[12]*B4 + c.[10]*B2) +
B6 %*% (c.[14]*dB6 + c.[12]*dB4 + c.[10]*dB2) +
c.[8]*dB6 + c.[6]*dB4 + c.[4]*dB2)
dV <-  dB6 %*% (c.[13]*B6 + c.[11]*B4 + c.[9]*B2) +
B6 %*% (c.[13]*dB6 + c.[11]*dB4 + c.[9]*dB2) +
c.[7]*dB6 + c.[5]*dB4 + c.[3]*dB2

dX <- -solve(V-U, dV-dU)%*% X + solve(V-U, dV + dU)
}
if (s > 0) for (t in 1:s) {
dX <- X %*% dX + dX %*% X
X <- X %*% X
}
dmat[,,dk] <- dX
}
else if (s>0) for (t in 1:s) {
X <- X %*% X
}

}
if (nderiv > 0) list(P=X, dmat=dmat)
else X
}
```

## Try the survival package in your browser

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

survival documentation built on Aug. 24, 2021, 5:06 p.m.