# R/frailty.controldf.R In survival: Survival Analysis

#### Defines functions frailty.controldf

```# \$Id: frailty.controldf.S 11373 2009-10-28 17:12:59Z therneau \$
# A function to calibrate the df
#    very empirical
# Find the closest 3 points that span the target value
#   We know the function is monotone, so fit the function
#     dy = a * (dx)^p   to the 3 points, where dx and dy are the distance
#     from the leftmost of the three points.
#   This method can fail near a boundary, so use step halving if things don't
#     go well
# On input, parms\$df = target degrees of freedom
#           parms\$dfs, parms\$thetas = known values (usually 0,0)
#           parms\$guess = first guess
#
frailty.controldf <- function(parms, iter, old, df) {
if (iter==0) {
theta <- parms\$guess
return(list(theta=theta, done=FALSE,
history=cbind(thetas=parms\$thetas, dfs=parms\$dfs)))
}

eps <- parms\$eps
if (length(eps)==0) eps <- .1

thetas <- c(old\$history[,1], old\$theta)
dfs    <- c(old\$history[,2], df)
nx <- length(thetas)
if (nx==2) {
#linear guess based on first two
# but try extra hard to bracket the root
theta <- thetas[1] + (thetas[2]-thetas[1])*(parms\$df - dfs[1])/
(dfs[2] - dfs[1])
if (parms\$df > df) theta <- theta * 1.5
return(list(theta=theta, done=FALSE,
history=cbind(thetas=thetas, dfs=dfs), half=0))
}
else{
# Now, thetas= our guesses at theta
#  dfs = the degrees of freedom for each guess
done <- (iter>1 &&
(abs(dfs[nx]-parms\$df) < eps))

# look for a new minimum
x <- thetas
y <- dfs
target <- parms\$df

# How am I doing
if ( abs( (y[nx]-target)/(y[nx-1]-target)) > .6) doing.well <- FALSE
else doing.well <- TRUE

ord <- order(x)
if ((x[1]-x[2])*(y[1]-y[2]) >0)  y <- y[ord]  #monotone up
else  { #monotone down
y <- -1* y[ord]
target <- -target
}
x <- x[ord]

if (all(y>target)) b1 <- 1     #points 1:3 are the closest then
else if (all(y<target)) b1 <- nx-2
else {
b1 <- max((1:nx)[y <= target]) #this point below target, next above
if (!doing.well && (is.null(old\$half) ||  old\$half<2)) {
#try bisection
if (length(parms\$trace) && parms\$trace){
print(cbind(thetas=thetas, dfs=dfs))
cat("  bisect:new theta=" , format( mean(x[b1+0:1])),
"\n\n")
}
return(list(theta= mean(x[b1+0:1]),done=done,
history=cbind(thetas=thetas, dfs=dfs),
half=max(old\$half, 0) +1))
}
# use either b1,b1+1,b1+2 or  b1-1, b1, b1+1, whichever is better
#  better = midpoint of interval close to the target

if ((b1+1)==nx ||
(b1>1 &&  ((target -y[b1]) < (y[b1+1] -target))))
b1 <- b1-1
}

#now have the best 3 points
# fit them with a power curve anchored at the leftmost one
b2 <- b1 + 1:2
xx <- log(x[b2] - x[b1])
yy <- log(y[b2] - y[b1])
power <- diff(yy)/diff(xx)
a <- yy[1] - power*xx[1]
newx <- (log(target -y[b1]) - a)/power
if (length(parms\$trace) && parms\$trace){
print(cbind(thetas=thetas, dfs=dfs))
cat("  new theta=" , format(x[b1] + exp(newx)), "\n\n")
}
list(theta=x[b1] + exp(newx), done=done,
history=cbind(thetas=thetas, dfs=dfs), half=0)
}
}
```

## 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.