Nothing
## TO DO: roll back into bbmle?
## allow multiple 'transects'?
## (i.e. if two sets of parameters are given ...)
## * robustification
## print method
## allow manual override of ranges
## allow log-scaling?
## preserve order of parameters in 1D plots
## substitute values of parameters into full parameter vector
mkpar <- function(params,p,i) {
params[i] <- p
params
}
## get reasonable range for slice
## document what is done here
## implement upper bound
## more robust approaches;
## try not to assume anything about signs of parameters
## inherit bounds from fitted value
get_trange <- function(pars, ## baseline parameter values
i, ## focal parameter
fun, ## objective function
lower=-Inf, ## lower bound
upper=Inf, ## upper bound
cutoff=10, ## increase above min z-value
maxit=200, ## max number of iterations
steptype=c("mult","addprop"),
step=0.1) {
## step possibilities: multiplicative
## additive (absolute scale) [not yet implemented]
addabs <- NULL ## fix false positive test
steptype <- match.arg(steptype)
v <- v0 <- fun(pars)
lowval <- pars[i]
it <- 1
if (steptype=="addprop") step <- step*pars[i]
while (it<maxit && lowval>lower && v<(v0+cutoff)) {
lowval <- switch(steptype,
addabs,
addpropn=lowval-step,
mult=lowval*(1-step))
v <- fun(mkpar(pars,lowval,i))
it <- it+1
}
lowdev <- v
lowit <- it
upval <- pars[i]
it <- 1
v <- v0 <- fun(pars)
if (upval==0) upval <- 1e-4
while (it<maxit && v<(v0+cutoff)) {
upval <- switch(steptype,
addabs,
addpropn=lowval+step,
mult=lowval*(1+step))
v <- fun(mkpar(pars,upval,i))
## cat(it,upper,v,"\n")
it <- it+1
}
updev <- v
upit <- it
c(low_val=lowval,up_valr=upval,low_dev=lowdev,up_dev=updev,
low_it=lowit,up_it=upit)
}
get_all_trange <- function(params,fun,lower,upper,cutoff=10,...) {
arglist <- c(list(pars=params,fun=fun,cutoff=cutoff),list(...))
tranges <- t(mapply(FUN=get_trange,
seq(length(params)),
lower,
upper,
MoreArgs=arglist,
SIMPLIFY=TRUE))
rownames(tranges) <- names(params)
tranges
}
## generic function (S3)
slice <- function (x, dim=1, ...) {
UseMethod("slice")
}
slice.mle2 <- function(x, ...) {
ff <- x@minuslogl
## vectorize objective function: return minus log L
## (switching to deviance screws things up???)
ff2 <- function(p) {
do.call(ff,as.list(p))
}
slice0(coef(x),ff2, ...)
}
slice0 <- function(params,fun,dim=1,params2=NULL,...) {
if (dim==1) {
if (is.null(params2)) {
slice1D(params,fun,...)
} else {
slicetrans(params,params2,fun,...)
}
} else {
if (!is.null(params2)) stop("can't do transect in 2D")
slice2D(params,fun,...)
}
}
plot.slice <- function(x,...) {
switch(x$dim,xyplot(x,...),splom(x,...))
}
slicetrans <- function(params, params2, fun, extend=0.1, nt=401,
lower=-Inf, upper=Inf) {
## make sure 0/1 are included in parameter vector ...
np <- length(params)
extend <- rep(extend,length.out=2)
lower <- rep(lower,length.out=np)
upper <- rep(upper,length.out=np)
slicep <- sort(unique(c(0,1,seq(-extend[1],1+extend[2], length=nt))))
slicepars <- t(sapply(slicep, function(x) (1 - x) * params + x * params2))
OK <- apply(slicepars,1,function(x) all(x>=lower & x<=upper))
if (any(!OK)) {
warning("some parameter sets outside of bounds were removed")
slicep <- slicep[OK]
slicepars <- slicepars[OK,]
}
v <- apply(slicepars, 1, fun)
slices <- list(data.frame(var1="trans",x=slicep,z=v))
r <- list(slices=slices,params=params,params2=params2,dim=1)
class(r) <- "slice"
r
}
slice1D <- function(params,fun,nt=101,
lower=-Inf,
upper=Inf,
verbose=TRUE,
tranges=NULL,
fun_args=NULL,
...) {
npv <- length(params)
if (is.null(pn <- names(params))) pn <- seq(npv)
if (is.null(tranges)) {
tranges <- get_all_trange(params,fun,
rep(lower,length.out=npv),
rep(upper,length.out=npv),
...)
}
slices <- vector("list",npv)
for (i in 1:npv) {
tvec <- seq(tranges[i,1],tranges[i,2],length=nt)
if (verbose) cat(pn[i],"\n")
vtmp <- vapply(tvec,
function(t) {
do.call(fun,c(list(mkpar(params,t,i)),fun_args))
}, numeric(1))
slices[[i]] <- data.frame(var1=pn[i],x=tvec,z=vtmp)
}
r <- list(slices=slices,ranges=tranges,params=params,dim=1)
class(r) <- "slice"
r
}
## OLD slice method
## should probably roll this in as an option to profile
## include attribute, warning? draw differently (leave off
## conf. limit lines)
## slice <- function(fitted, ...) UseMethod("slice")
## 1D slicing implemented as in profile
sliceOld <- function (fitted, which = 1:p, maxsteps = 100,
alpha = 0.01, zmax = sqrt(qchisq(1 - alpha/2, p)),
del = zmax/5, trace = FALSE,
tol.newmin=0.001, ...)
{
onestep <- function(step)
{
bi <- B0[i] + sgn * step * del * std.err[i]
fix <- list(bi)
names(fix) <- p.i
call$fixed <- c(fix,eval(call$fixed))
call$eval.only = TRUE
pfit <- try(eval(call), silent=TRUE) ##
if(inherits(pfit, "try-error")) return(NA)
else {
zz <- 2*(pfit@min - fitted@min)
ri <- pv0
ri[, names(pfit@coef)] <- pfit@coef
ri[, p.i] <- bi
if (zz > -tol.newmin)
zz <- max(zz, 0)
else stop("profiling has found a better solution, so original fit had not converged")
z <- sgn * sqrt(zz)
pvi <<- rbind(pvi, ri)
zi <<- c(zi, z) ## NB global set!
}
if (trace) cat(bi, z, "\n")
z
}
## Profile the likelihood around its maximum
## Based on profile.glm in MASS
summ <- summary(fitted)
std.err <- summ@coef[, "Std. Error"]
Pnames <- names(B0 <- fitted@coef)
pv0 <- t(as.matrix(B0))
p <- length(Pnames)
prof <- vector("list", length = length(which))
names(prof) <- Pnames[which]
call <- fitted@call
call$minuslogl <- fitted@minuslogl
for (i in which) {
zi <- 0
pvi <- pv0
p.i <- Pnames[i]
for (sgn in c(-1, 1)) {
if (trace)
cat("\nParameter:", p.i, c("down", "up")[(sgn + 1)/2 + 1], "\n")
step <- 0
z <- 0
## This logic was a bit frail in some cases with
## high parameter curvature. We should probably at least
## do something about cases where the mle2 call fails
## because the parameter gets stepped outside the domain.
## (We now have.)
call$start <- as.list(B0)
lastz <- 0
while ((step <- step + 1) < maxsteps && abs(z) < zmax) {
z <- onestep(step)
if(is.na(z)) break
lastz <- z
}
if(abs(lastz) < zmax) {
## now let's try a bit harder if we came up short
for(dstep in c(0.2, 0.4, 0.6, 0.8, 0.9)) {
z <- onestep(step - 1 + dstep)
if(is.na(z) || abs(z) > zmax) break
}
} else if(length(zi) < 5) { # try smaller steps
mxstep <- step - 1
step <- 0.5
while ((step <- step + 1) < mxstep) onestep(step)
}
}
si <- order(pvi[, i])
prof[[p.i]] <- data.frame(z = zi[si])
prof[[p.i]]$par.vals <- pvi[si,, drop=FALSE]
}
list(profile = prof, summary = summ)
}
## * is it possible to set up the 2D vectors so they include
## the baseline value? maybe not easily ...
slice2D <- function(params,
fun,
nt=31,
lower=-Inf,
upper=Inf,
cutoff=10,
verbose=TRUE,
tranges=NULL,
...) {
npv <- length(params)
if (is.null(pn <- names(params))) pn <- seq(npv)
if (is.null(tranges)) {
tranges <- get_all_trange(params,fun,
rep(lower,length.out=npv),
rep(upper,length.out=npv),
cutoff=cutoff,
...)
}
slices <- list()
for (i in 1:(npv-1)) {
slices[[i]] <- vector("list",npv)
for (j in (i+1):npv) {
if (verbose) cat("param",i,j,"\n")
t1vec <- seq(tranges[i,1],tranges[i,2],length=nt)
t2vec <- seq(tranges[j,1],tranges[j,2],length=nt)
mtmp <- matrix(nrow=nt,ncol=nt)
for (t1 in seq_along(t1vec)) {
for (t2 in seq_along(t2vec)) {
mtmp[t1,t2] <- fun(mkpar(params,c(t1vec[t1],t2vec[t2]),
c(i,j)))
}
}
slices[[i]][[j]] <- data.frame(var1=pn[i],var2=pn[j],
expand.grid(x=t1vec,y=t2vec),
z=c(mtmp))
}
}
r <- list(slices=slices,ranges=tranges,params=params,dim=2)
class(r) <- "slice"
r
}
## flatten slice:
## do.call(rbind,lapply(slices,do.call,what=rbind))
slices_apply <- function(s,FUN,...) {
for (i in seq_along(s)) {
for (j in seq_along(s[[i]])) {
if (!is.null(s[[i]][[j]])) {
s[[i]][[j]] <- FUN(s[[i]][[j]],...)
}
}
}
s
}
xyplot.slice <- function(x,data,type="l",scale.min=TRUE,...) {
allslice <- do.call(rbind,x$slices)
## correct ordering
allslice$var1 <- factor(allslice$var1,
levels=unique(as.character(allslice$var1)))
if (scale.min) allslice$z <- allslice$z-min(allslice$z)
pfun <- function(x1,y1,...) {
panel.xyplot(x1,y1,...)
if (is.null(x$params2)) {
## regular 1D slice
panel.abline(v=x$params[panel.number()],col="gray")
} else {
## 'transect' slice
panel.abline(v=c(0,1),col="gray")
panel.abline(h=y1[x1 %in% c(0,1)],col="gray")
}
}
xyplot(z~x|var1,data=allslice,type=type,
scales=list(x=list(relation="free")),
panel=pfun,...)
}
splom.slice <- function(x,
data,
scale.min=TRUE,
at=NULL,
which.x=NULL,
which.y=NULL,
dstep=4,
contour=FALSE,...) {
if (x$dim==1) stop("can't do splom on 1D slice object")
smat <- t(x$ranges[,1:2])
if (scale.min) {
## FIXME: something more elegant to flatten slice list?
all.z <- unlist(sapply(x$slices,
function(x) {
sapply(x,
function(x)
if (is.null(x)) NULL else x[["z"]])
}))
min.z <- min(all.z[is.finite(all.z)])
## round up to next multiple of 'dstep'
max.z <- dstep * ((max(all.z[is.finite(all.z)])-
min.z) %/% dstep + 1)
if (missing(at)) {
at <- seq(0,max.z,by=dstep)
}
scale.z <- function(X) {
X$z <- X$z-min.z
X
}
x$slices <- slices_apply(x$slices,scale.z)
}
up0 <- function(x1, y, groups, subscripts, i, j, ...) {
## browser()
sl <- x$slices[[j]][[i]]
with(sl,panel.levelplot(x=x,y=y,z=z,contour=contour,
at=if (!is.null(at)) at else pretty(z),
subscripts=seq(nrow(sl))))
panel.points(x$params[j],x$params[i],pch=16)
mm <- matrix(sl$z,nrow=length(unique(sl$x)))
## FIXME: more robust ...
wmin <- which(mm==min(mm),arr.ind=TRUE)
xmin <- unique(sl$x)[wmin[1]]
ymin <- unique(sl$y)[wmin[2]]
panel.points(xmin,ymin,pch=1)
}
lp0 <- function(...) {
}
## FIXME: use ?draw.colorkey to add a legend ...
## FIXME: make diagonal panel text smaller ???
splom(smat,lower.panel=lp0,diag.panel=diag.panel.splom,
upper.panel=up0,...)
}
## generic profiling code???
## either need (1) optimizers with 'masks' or (2)
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.