Nothing
pentrace <-
function(fit, penalty, penalty.matrix,
method=c('grid', 'optimize'),
which=c('aic.c', 'aic', 'bic'), target.df=NULL,
fitter, pr=FALSE,
tol=1e-7, keep.coef=FALSE, complex.more=TRUE,
verbose=FALSE, maxit=12, subset, noaddzero=FALSE)
{
## Need to check Strata for cph
method <- match.arg(method)
which <- match.arg(which)
tdf <- length(target.df)
if(tdf) method <- 'optimize'
if(! length(X <- fit$x) || ! length(Y <- as.matrix(fit$y)))
stop("you did not specify x=TRUE and y=TRUE in the fit")
fit$x <- fit$y <- NULL
## if(length(pn <- fit$penalty) > 0 && max(unlist(pn)) != 0)
## warning('you did not specify penalty= in fit so that unpenalized model can be a candidate for the best model')
sc.pres <- match("parms", names(fit), 0) > 0
fixed <- NULL
dist <- fit$dist
parms <- fit$parms
clas <- class(fit)[1]
isols <- clas=='ols'
if(!(isols || inherits(fit, 'lrm')))
stop("at present pentrace only works for lrm or ols")
if(missing(fitter))
fitter <- switch(clas,
ols=function(x, y, maxit, ...) lm.pfit(x, y, ...),
lrm=function(x, y, maxit=12, ...)
lrm.fit(x, y, maxit=maxit, ...),
cph=function(x, y, maxit=12, ...) coxphFit(x, y,
strata=Strata, iter.max=maxit,
eps=.0001, method="efron",
toler.chol=tol),
psm=function(x, y, maxit=12,...) survreg.fit2(x, y,
dist=dist, parms=parms, fixed=fixed,
offset=NULL,
init=NULL, maxiter=maxit))
if(!length(fitter))stop("fitter not valid")
Strata <- fit$strata
if(!missing(subset)) {
Y <- Y[subset,, drop=FALSE]
X <- X[subset,, drop=FALSE]
Strata <- Strata[subset, drop=FALSE] # NULL[] is NULL
}
n <- nrow(Y)
atr <- fit$Design
if(missing(penalty.matrix))
penalty.matrix <- Penalty.matrix(atr, X)
obj.best <- -1e10
ns <- num.intercepts(fit)
islist <- is.list(penalty)
if(islist) {
penalty <- expand.grid(penalty)
if(complex.more && ncol(penalty) > 1 && nrow(penalty) > 1) {
ikeep <- NULL
for(i in 1:nrow(penalty)) {
ok <- TRUE
peni <- penalty[i,]
for(j in 2:length(peni))
if(peni[[j]] < peni[[j-1]]) ok <- FALSE
if(ok) ikeep <- c(ikeep, i)
}
penalty <- penalty[ikeep,,drop=FALSE]
}
np <- nrow(penalty)
}
else {
if(method == 'grid' && ! noaddzero) penalty <- c(0, penalty[penalty > 0])
np <- length(penalty)
}
if(method=='optimize') {
stop('method="optimize" not yet implemented in R')
if((islist && nrow(penalty) > 1) || (!islist && length(penalty) > 1))
stop('may not specify multiple potential penalties when method="optimize"')
objective <- function(pen, X, Y, z) {
##Problem with sending so many auxiliary parameters to nlminb -
##nlminb's internal parameters got shifted somehow
n <- z$n; penalty.matrix <- z$penalty.matrix; pennames <- z$pennames
isols <- z$isols; islist <- z$islist; tol <- z$tol; maxit <- z$maxit
ns <- z$ns; fitter <- z$fitter; pr <- z$pr; atr <- z$atr;
tdf <- length(z$target.df)
if(length(pen) > 1) {
pen <- structure(as.list(pen), names=pennames)
penfact <- Penalty.setup(atr, pen)$multiplier
} else penfact <- pen
if(length(penfact)==1 || !islist) pm <- penfact*penalty.matrix
else {
a <- diag(sqrt(penfact))
pm <- a %*% penalty.matrix %*% a
}
f <- fitter(X, Y, penalty.matrix=pm, tol=tol, maxit=maxit)
if(length(f$fail) && f$fail)
stop('fitter failed. Try changing maxit or tol')
if(isols) {
## ols (from lm.pfit) already stored correct LR chisq and effective df
stats <- f$stats
df <- stats['d.f.']
lr <- stats['Model L.R.']
dag <- f$effective.df.diagonal
}
else {
v <- f$var #Later: vcov(f)
f.nopenalty <- fitter(X, Y, initial=f$coef, maxit=1, tol=tol)
if(length(f.nopenalty$fail) && f.nopenalty$fail)
stop('fitter failed. Try changing tol')
info.matrix.unpenalized <-
if(length(f.nopenalty$info.matrix))
f.nopenalty$info.matrix else
solvet(f.nopenalty$var, tol=tol) # -> vcov
dag <- diag(info.matrix.unpenalized %*% v)
df <- if(ns==0)sum(dag) else sum(dag[-(1:ns)])
lr <- f.nopenalty$stats["Model L.R."]
}
obj <- switch(z$which,
aic.c <- lr - 2*df*(1 + (df + 1) / (n - df - 1)),
aic <- lr - 2 * df,
bic <- lr - df * logb(n))
if(tdf) obj <- abs(df - z$target.df)
if(pr) {
w <- if(tdf) df else obj
names(w) <- NULL
pp <- if(islist) unlist(pen) else c(Penalty=pen)
print(c(pp, Objective=w))
}
if(!tdf) obj <- -obj
else
attr(obj,'df') <- df
obj
}
res <- nlminb(unlist(penalty), objective, lower=0, X=X, Y=Y,
z=list(n=n,
penalty.matrix=penalty.matrix, pennames=names(penalty),
isols=isols, islist=islist, tol=tol, maxit=maxit, ns=ns,
fitter=fitter, atr=atr, pr=pr, which=which,
target.df=target.df),
control=list(abs.tol=.00001,
rel.tol=if(tdf)1e-6 else .00001))
return(list(penalty=if(islist)
structure(as.list(res$parameters),names=names(penalty))
else res$parameters,
objective=if(tdf)res$aux$df else -res$objective))
}
df <- aic <- bic <- aic.c <-
if(islist) double(length(penalty[[1]])) else double(length(penalty))
for(i in 1 : np) {
if(islist) {
pen <- penalty[i,]
penfact <- Penalty.setup(atr, pen)$multiplier
} else {
pen <- penalty[i]
penfact <- pen
}
unpenalized <- all(penfact==0)
if(i==1) Coef <- if(keep.coef) matrix(NA,ncol=length(fit$coef),nrow=np)
else NULL
if(unpenalized) f <- fit
else {
if(length(penfact) == 1 || !islist) pm <- penfact * penalty.matrix
else {
a <- diag(sqrt(penfact))
pm <- a %*% penalty.matrix %*% a
}
f <- fitter(X, Y, penalty.matrix=pm, tol=tol, maxit=maxit)
if(length(f$fail) && f$fail)
stop('fitter failed. Try changing maxit or tol')
}
if(keep.coef) Coef[i,] <- f$coef
if(unpenalized || isols) {
## ols (from lm.pfit) already stored correct LR chisq and effective df
stats <- f$stats
df[i] <- stats['d.f.']
lr <- stats['Model L.R.']
dag <- if(unpenalized) rep(1, length(df[i]))
else
f$effective.df.diagonal
}
else {
v <- f$var #Later: vcov(f, regcoef.only=T)
f.nopenalty <- fitter(X, Y, initial=f$coef, maxit=1, tol=tol)
if(length(f.nopenalty$fail) && f.nopenalty$fail)
stop('fitter failed. Try changing tol')
info.matrix.unpenalized <-
if(length(f.nopenalty$info.matrix)) f.nopenalty$info.matrix
else
solvet(f.nopenalty$var, tol=tol) # -> vcov
dag <- diag(info.matrix.unpenalized %*% v)
df[i] <- if(ns == 0)sum(dag) else sum(dag[- (1 : ns)])
lr <- f.nopenalty$stats["Model L.R."]
if(verbose) {
cat('non slopes',ns,'\neffective.df.diagonal:\n')
print(dag)
}
}
aic[i] <- lr - 2 * df[i]
bic[i] <- lr - df[i] * logb(n)
aic.c[i] <- lr - 2 * df[i] * (1 + (df[i] + 1) / (n - df[i] - 1))
obj <- switch(which, aic.c=aic.c[i], aic=aic[i], bic=bic[i])
if(obj > obj.best) {
pen.best <- pen
df.best <- df[i]
obj.best <- obj
f.best <- f
var.adj.best <- if(unpenalized || isols) f$var
else
v %*% info.matrix.unpenalized %*% v
diag.best <- dag
}
if(pr) {
d <- if(islist) as.data.frame(pen, row.names='') else
data.frame(penalty=pen, row.names='')
d$df <- df[i]
d$aic <- aic[i]
d$bic <- bic[i]
d$aic.c <- aic.c[i]
cat('\n'); print(d)
}
}
mat <- if(islist) as.data.frame(penalty)
else
data.frame(penalty=penalty)
mat$df <- df
mat$aic <- aic
mat$bic <- bic
mat$aic.c <- aic.c
structure(list(penalty=pen.best, df=df.best, objective=obj.best,
fit=f.best, var.adj=var.adj.best, diag=diag.best,
results.all=mat, Coefficients=Coef), class="pentrace")
}
plot.pentrace <- function(x, method=c('points', 'image'),
which=c('effective.df', 'aic', 'aic.c', 'bic'),
pch=2, add=FALSE, ylim, ...)
{
method <- match.arg(method)
x <- x$results.all
penalty <- x[[1]]
effective.df <- x$df
aic <- x$aic
bic <- x$bic
aic.c <- x$aic.c
if(length(x) == 5) { ## only one variable given to penalty=
if('effective.df' %in% which) {
if(add) lines(penalty, effective.df) else
plot(penalty, effective.df, xlab="Penalty", ylab="Effective d.f.",
type="l", ...)
if(length(which) == 1) return(invisible())
}
if(!add) plot(penalty, aic,
ylim=if(missing(ylim)) range(c(aic, bic)) else ylim,
xlab="Penalty",
ylab=expression(paste("Information Criterion (", chi^2,
" scale)")),
type=if('aic' %in% which)"l" else "n", lty=3, ...)
else
if('aic' %in% which) lines(penalty, aic, lty=3, ...)
if('bic' %in% which) lines(penalty, bic, lty=2, ...)
if('aic.c' %in% which) lines(penalty, aic.c, lty=1, ...)
if(!add && length(setdiff(which, 'effective.df')) > 1)
title(sub=paste(if('aic.c' %in% which) "Solid: AIC_c",
if('aic' %in% which) "Dotted: AIC",
if('bic' %in% which) "Dashed: BIC",sep=' '),
adj=0,cex=.75)
return(invisible())
}
## At least two penalty factors
if(add) stop('add=TRUE not implemented for >=2 penalty factors')
X1 <- x[[1]]
X2 <- x[[2]]
nam <- names(x)
x1 <- sort(unique(X1))
x2 <- sort(unique(X2))
n1 <- length(x1)
n2 <- length(x2)
aic.r <- rank(aic); aic.r <- aic.r/max(aic.r)
if(method=='points') {
plot(0, 0, xlim=c(1,n1), ylim=c(1,n2), xlab=nam[1], ylab=nam[2],
type='n', axes=FALSE, ...)
mgp.axis(1, at=1:n1, labels=format(x1))
mgp.axis(2, at=1:n2, labels=format(x2))
ix <- match(X1, x1)
iy <- match(X2, x2)
for(i in 1:length(aic)) points(ix[i], iy[i], pch=pch, cex=(.1+aic.r[i])*3)
return(invisible(aic.r))
}
z <- matrix(NA,nrow=n1,ncol=n2)
for(i in 1:n1)
for(j in 1:n2) z[i,j] <- aic.r[X1==x1[i] & X2==x2[j]]
image(x1, x2, z, xlab=nam[1], ylab=nam[2], zlim=range(aic.r), ...)
invisible(aic.r)
}
print.pentrace <- function(x, ...)
{
cat('\nBest penalty:\n\n')
pen <- if(is.list(x$penalty)) as.data.frame(x$penalty,row.names='')
else
data.frame(penalty=x$penalty, row.names='')
pen$df <- x$df
pen$aic <- x$aic
print(pen)
cat('\n')
if(is.data.frame(x$results.all)) print(x$results.all, row.names=FALSE) else
print(as.data.frame(x$results.all,), row.names=FALSE)
# row.names=rep('',length(x$results.all[[1]]))))
invisible()
}
effective.df <- function(fit, object)
{
atr <- fit$Design
dag <- if(missing(object)) fit$effective.df.diagonal
else
object$diag
if(length(dag)==0) stop('object not given or fit was not penalized')
ia.or.nonlin <- param.order(atr, 2)
nonlin <- param.order(atr, 3)
ia <- param.order(atr, 4)
ia.nonlin <- param.order(atr, 5)
ns <- num.intercepts(fit)
if(ns > 0) dag <- dag[-(1:ns)]
z <- rbind(c(length(dag), sum(dag)),
c(sum(!ia.or.nonlin), sum(dag[!ia.or.nonlin])),
c(sum(ia.or.nonlin), sum(dag[ia.or.nonlin])),
c(sum(nonlin), sum(dag[nonlin])),
c(sum(ia), sum(dag[ia])),
c(sum(ia.nonlin), sum(dag[ia.nonlin])))
dimnames(z) <- list(c('All','Simple Terms','Interaction or Nonlinear',
'Nonlinear', 'Interaction','Nonlinear Interaction'),
c('Original','Penalized'))
cat('\nOriginal and Effective Degrees of Freedom\n\n')
print(round(z,2))
invisible(z)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.