Revision 56d752cf6d1a8e2775e1acc46482c206cfd5b269 authored by Adrian Baddeley on 14 March 2012, 07:09:51 UTC, committed by cran-robot on 14 March 2012, 07:09:51 UTC
1 parent a4d492a
ippm.R
#
# ippm.R
#
# $Revision: 2.7 $ $Date: 2011/12/03 11:37:03 $
#
# Fisher scoring algorithm for irregular parameters in ppm trend
#
ippm <- function(...,
iScore=NULL,
start=list(),
covfunargs=start,
maxiter=20, tol=1e-4, progress=TRUE, stepfactor=1,
dbug=FALSE) {
# validate
if(!is.list(start) || length(start) == 0)
stop("start should be a list of initial values for irregular parameters")
if(!is.list(iScore) || length(iScore) != length(start))
stop("iScore should be a list of the same length as start")
stopifnot(identical(names(iScore), names(start)))
if(!all(unlist(lapply(iScore, is.function))))
stop("iScore should be a list of functions")
#
smap <- match(names(start), names(covfunargs))
if(any(is.na(smap)))
stop("variables in start should be a subset of variables in covfunargs")
covfunargs[smap] <- start
#
# fit the initial model and extract information
fit0 <- ppm(..., covfunargs=covfunargs)
lpl0 <- fit0$maxlogpl
p <- length(coef(fit0))
# examine covariates and trend
covariates <- fit0$covariates
isfun <- unlist(lapply(covariates, is.function))
covfuns <- covariates[isfun]
# determine which covariates depend on which irregular parameters
pnames <- names(start)
hasarg <- function(f,a) { a %in% names(formals(f)) }
depmat <- matrix(FALSE, nrow=length(covfuns), ncol=length(pnames))
rownames(depmat) <- names(covfuns)
colnames(depmat) <- pnames
for(j in 1:length(pnames))
depmat[,j] <- unlist(lapply(covfuns, hasarg, pnames[j]))
# find covariates that depend on ANY irregular parameter
depvar <- rownames(depmat)[apply(depmat, 1, any)]
# check that these covariates appear only in offset terms
covnames.fitted <- model.covariates(fit0, fitted=TRUE, offset=FALSE)
if(any(uhoh <- depvar %in% covnames.fitted))
stop(paste(ngettext(sum(uhoh), "The covariate", "The covariates"),
commasep(sQuote(depvar[uhoh])),
"should appear only in offset terms"))
# check that every irregular parameter to be updated appears somewhere
covnames.offset <- model.covariates(fit0, fitted=FALSE, offset=TRUE)
usearg <- apply(depmat[covnames.offset, , drop=FALSE], 2, any)
if(!all(usearg)) {
nbad <- sum(!usearg)
warning(paste("Cannot maximise over the irregular",
ngettext(nbad, "parameter", "parameters"),
commasep(sQuote(names(usearg)[!usearg])),
ngettext(nbad, "because it is", "because they are"),
"not used in any term of the model"))
# restrict
start <- start[usearg]
iScore <- iScore[usearg]
pnames <- names(start)
}
# ready
iter <- 0
param <- start
pvec <- as.numeric(param)
ndigits <- max(0, -ceiling(log10(tol))) + 1
# go
for(phase in 1:2) {
if(progress) cat(paste("Phase", phase, "\n"))
maxit <- if(phase == 1) 2 else maxiter
for(iter in 0:maxit) {
# fit model with current irregular parameters
covfunargs[smap] <- param
fit <- ppm(..., covfunargs=covfunargs)
lpl <- logLik(fit, warn=FALSE)
if(progress) {
co <- coef(fit)
cat(paste(paren(iter, "["),
paste(paste(pnames, "=", round(pvec, digits=ndigits)),
collapse=", "),
"->",
paste(paste(names(co), "=", round(co, digits=ndigits)),
collapse=", "),
"; logPL=", signif(lpl, 8), "\n"))
}
# compute model matrix and inverse fisher information
stuff <- ppm.influence(fit, what="derivatives",
iScore=iScore,
iArgs=param)
score <- stuff$deriv$score
vc <- stuff$deriv$vc
fush <- stuff$deriv$fush
if(dbug) {
cat("\nscore=\n")
print(score)
cat("\nvc=\n")
print(vc)
cat("\nfush=\n")
print(fush)
}
if(phase == 1) {
# Fisher scoring on partial matrix
subscore <- score[, -(1:p), drop=FALSE]
subfush <- fush[-(1:p), -(1:p)]
stepvec <- as.numeric(solve(subfush) %*% t(subscore))
} else {
# Fisher scoring step on full parameter vector
stepvec <- as.numeric(vc %*% t(score))[ -(1:p)]
}
if(dbug) {
cat("\nstep=\n")
print(stepvec)
}
# update parameters or exit
if(iter > 0 && all(abs(stepvec) < tol))
break
pvec <- pvec + stepfactor * stepvec
param <- as.list(pvec)
names(param) <- pnames
}
}
if(iter == maxiter)
warning(paste("Maximum number of iterations", paren(maxiter),
"reached without convergence\n"))
return(fit)
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...