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

#### Defines functions frailty.controlgauss

```# \$Id: frailty.controlgauss.S 11166 2008-11-24 22:10:34Z therneau \$
#
# The control function for REML on a gaussian
#
frailty.controlgauss <- function(opt, iter, old, fcoef, trH, loglik){
if (iter==0) {
# initial call
#  Because of how the iteration works, 0 is not a useful trial value
if (!is.null(opt\$theta)) theta <- opt\$theta  #fixed theta case
else {
if (is.null(opt\$init)) theta <- 1
else theta <- opt\$init
}
list(theta=theta)
}

else {
if (is.null(opt\$trace)) trace <-FALSE
else trace <- opt\$trace

nfrail <- length(fcoef)
fsum   <- sum(fcoef^2)
theta <- old\$theta
resid <- fsum/(nfrail - trH/theta) - theta

# save history of the iteration, and get the next theta
if (iter==1) {
history <- c(theta=theta, resid=resid, fsum=fsum, trace=trH)
if (is.null(opt\$init )) {
if (resid>0)  theta <- theta*3
else          theta <- theta/3  	        }
else theta <- opt\$init
list(theta=theta, done=FALSE, history=history)
}
else {
history <- rbind(old\$history,
as.vector(c(theta, resid, fsum, trH)))
if (iter ==2) {
if (all(history[,2] > 0))     theta <- history[2,1]*2
else if (all(history[,2] <0)) theta <- history[2,1]/2
else    		theta <- mean(history[1:2,1])
if (trace) {
print(history)
cat("    new theta=", theta, "\n\n")
}
list(theta=theta, done=FALSE, history=history)
}
else {
done <- (abs(history[iter,2]) < opt\$eps)
ord <- order(history[,1])
tempy <- history[ord,2]  #x & y from left to right
tempx <- history[ord,1]
# make sure we have one positve and one negative y value
#  y must be positive near 0, and negative for large x
if (all(tempy>0))  newtheta <- 2*max(tempx)
else if (all(tempy<0)) newtheta <- .5 * min(tempx)
else{
#find the latest point, and one on each side of 0
b1 <- (1:iter)[ord==iter]
if (b1==1) b1 <-2
else if (b1==iter) b1 <- iter-1

# Brent's formula, straight from Numerical Recipies
# Why, you may ask, don't we use the uniroot() function
#  which is built into S, and implements Brent's method?
# Because all we want is the next guess for x.  The interal
#  loop of coxph is calling us, not the other way around.
R <- tempy[b1]/ tempy[b1+1]
S <- tempy[b1]/ tempy[b1-1]
U <- R/S
P <- S* (U*(R-U)*(tempx[b1+1]-tempx[b1]) -
(1-R)*(tempx[b1]-tempx[b1-1]))
Q <- (U-1)*(R-1)*(S-1)
newtheta <- tempx[b1] + P/Q
# if the new guess is outside the brackets, do a binomial
#   search step
if (newtheta > tempx[b1+1]) newtheta <- mean(tempx[b1+0:1])
if( newtheta < tempx[b1-1]) newtheta <- mean(tempx[b1-0:1])
}
if (trace) {
print(history)
cat("    new theta=", format(newtheta), "\n\n")
}
list(theta=newtheta, done=done, history=history)
}
}
}
}

```

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