1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
|
# $Id: frailty.controlaic.S 11166 2008-11-24 22:10:34Z therneau $
# Control function to minimize the AIC
# the optional paramater "caic" chooses corrected aic (default=FALSE)
# n is the "effective" sample size
#
frailty.controlaic <- function(parms, iter, old, n, df, loglik) {
if (iter==0) { # initial call
if (is.null(parms$init)) theta <-0.005
else theta <- parms$init[1]
return(list(theta=theta, done=FALSE))
}
# by default, do the corrected AIC
if (length(parms$caic)) correct <- parms$caic
else correct <- FALSE
if (n < df+2) dfc <- (df -n) + (df+1)*df/2 -1 #avoid pathology
else dfc <- -1 + (df+1)/(1- ((df+2)/n))
if (iter==1) { # Second guess in series
history <- c(theta=old$theta, loglik=loglik,
df=df, aic=loglik-df, aicc=loglik - dfc)
if (length(parms$init) <2) theta <-1
else theta <- parms$init[2]
temp <- list(theta=theta, done=FALSE, history=history)
return(temp)
}
history <- rbind(old$history,c(old$theta, loglik, df, loglik-df,
loglik -dfc))
if (is.null(parms$trace)) trace <-FALSE
else trace <- parms$trace
if (iter==2) { #Third guess
theta <- mean(history[,1])
return(list(theta=theta, done=FALSE, history=history))
}
#
# Ok, now we're ready to actually use prior data
# Now, history has iter rows, each row contains the
# value of theta, the Cox PL, the df, aic, and corrected aic
if (correct) aic <- history[,5] #use corrected aic for convergence
else aic <- history[,4]
done <- (abs(1- aic[iter]/aic[iter-1]) < parms$eps)
x <- history[,1]
if (x[iter]== max(aic) && x[iter]==max(x))
newtheta <- 2* max(x)
else newtheta <- frailty.brent(x, aic, lower=parms$lower,
upper=parms$upper)
if (length(parms$trace) && parms$trace) {
print(history)
cat(" new theta=", format(newtheta), "\n\n")
}
list(theta=newtheta, done=done, history=history)
}
|
