# SCCS @(#)frailty.controlgam.s 1.1 10/28/98 # # The control function for a single Gamma frailty term. # frailty.controlgam <- function(opt, iter, old, group, status, loglik){ if (iter==0) { # initial call if (!is.null(opt$theta)) theta <- opt$theta #fixed theta case else { if (is.null(opt$init)) theta <- 0 #no initial value -- use 0 else theta <- opt$init[1] } list(theta=theta) } else { if (is.null(opt$trace)) trace <-FALSE else trace <- opt$trace theta <- old$theta #compute correction to the loglik if (theta==0) correct <- 0 else { if (is.matrix(group)) group <- c(group %*% 1:ncol(group)) d <- tapply(status,group,sum) correct <- frailty.gammacon(d, 1/theta) } if (!is.null(opt$theta)) # fixed theta case list(theta=theta, done=TRUE, c.loglik=loglik + correct) else { # save history of the iteration, and get the next theta if (iter==1) history <- c(theta=theta, loglik=loglik, c.loglik=loglik + correct) else history <- rbind(old$history, as.vector(c(theta, loglik, loglik + correct))) if (iter==1) { if (is.null(opt$init )) theta <-1 else theta <- opt$init[2] list(theta=theta, done=FALSE, history=history, c.loglik= loglik+correct) } else if (iter ==2) { if (history[2,3] < (history[1,3] +1)) theta <- mean(history[1:2,1]) else theta <- 2*history[2,1] if (trace) { print(history) cat(" new theta=", theta, "\n\n") } list(theta=theta, done=FALSE, history=history, c.loglik= loglik+correct) } else { #Now, history has iter rows, each row contains the value # of theta, the Cox PL, and the full LL done <- (abs(1- history[iter,3]/history[iter-1,3]) < opt$eps) x <- history[,1] y <- history[,3] if (y[iter]== max(y) && x[iter]==max(x)) newtheta <- 2* max(x) else newtheta <- frailty.brent(sqrt(x), y, lower=0)^2 if (trace) { print(history) cat(" new theta=", format(newtheta), "\n\n") } list(theta=newtheta, done=done, history=history, c.loglik = loglik + correct) } } } }