Nothing
################################################################################
### profile-method for class "twinSIR" to calculate the profile log-likelihood
### (normalized) as well as profile likelihood based confidence intervals
###
### Copyright (C) 2009 Michael Hoehle, 2014 Sebastian Meyer
###
### This file is part of the R package "surveillance",
### free software under the terms of the GNU General Public License, version 2,
### a copy of which is available at https://www.R-project.org/Licenses/.
################################################################################
######################################################################
# Function to compute likelihood based confidence interval, basically
# the two solutions to
# f(\theta) = l(\theta)-l(\hat{theta)) + 1/2 dchisq(1-alpha,df=1)=0
# are found.
#
#
# Parameters:
# logliktilde - normalized likelihood function(theta, ...)
# theta.hat - the MLE
# lower - search interval [lower,theta.hat] for f=0
# upper - search interval [theta.hat,upper] for f=0
# alpha - confidence level (see Equation 2.6 in Pawitan (2003)
# ... - additional arguments passed to function logliktilde
######################################################################
likelihood.ci <- function (logliktilde, theta.hat, lower, upper,
alpha = 0.05, ...)
{
# Highest Likelihood interval -- target function
f <- function(theta, ...) {
logliktilde(theta, ...) + 1/2*qchisq(1-alpha, df=1)
}
# Compute upper and lower boundary numerically
hl.lower <- uniroot(f, interval = c(lower, theta.hat), ...)$root
hl.upper <- uniroot(f, interval = c(theta.hat, upper), ...)$root
return(c(hl.lower,hl.upper))
}
######################################################################
# Function to compute estimated and profile likelihood based
# confidence intervals. Heavy computations might be necessary!
#
#Params:
# fitted - output from a fit with twinSIR
# profile - list with 4D vector as entries - format:
# c(index, lower, upper, grid size)
# where index is the index in the coef vector
# lower and upper are the parameter limits (can be NA)
# grid size is the grid size of the equally spaced grid
# between lower and upper (can be 0)
# alpha - (1-alpha)% profile likelihood CIs are computed.
# If alpha <= 0 then no CIs are computed
# control - control object to use for optim in the profile loglik computations
#
# Returns:
# list with profile loglikelihood evaluations on the grid
# and highest likelihood and wald confidence intervals
######################################################################
profile.twinSIR <- function (fitted, profile, alpha = 0.05,
control = list(fnscale = -1, factr = 1e1, maxit = 100), ...)
{
## Check that input is ok
profile <- as.list(profile)
if (length(profile) == 0L) {
stop("nothing to do")
}
lapply(profile, function(one) {
if (length(one) != 4L) {
stop("each profile entry has to be of form ",
"'c(index, lower, upper, grid size)'")
}})
if (is.null(fitted[["model"]])) {
stop("'fitted' must contain the model component")
}
px <- ncol(fitted$model$X)
pz <- ncol(fitted$model$Z)
## Control of the optim procedure
if (is.null(control[["fnscale",exact=TRUE]])) { control$fnscale <- -1 }
if (is.null(control[["factr",exact=TRUE]])) { control$factr <- 1e1 }
if (is.null(control[["maxit",exact=TRUE]])) { control$maxit <- 100 }
## Estimated normalized likelihood function
ltildeestim <- function(thetai,i) {
theta <- theta.ml
theta[i] <- thetai
with(fitted$model,
.loglik(theta, X=X, Z=Z, survs=survs, weights=weights)) - loglik.theta.ml
}
## Profile normalized likelihood function
ltildeprofile <- function(thetai,i)
{
emptyTheta <- rep(0, length(theta.ml))
# Likelihood l(theta_{-i}) = l(theta_i, theta_i)
ltildethetaminusi <- function(thetaminusi) {
theta <- emptyTheta
theta[-i] <- thetaminusi
theta[i] <- thetai
with(fitted$model,
.loglik(theta, X=X, Z=Z, survs=survs, weights=weights)) - loglik.theta.ml
}
# Score function of all params except thetaminusi
stildethetaminusi <- function(thetaminusi) {
theta <- emptyTheta
theta[-i] <- thetaminusi
theta[i] <- thetai
with(fitted$model,
.score(theta, X=X, Z=Z, survs=survs, weights=weights))[-i]
}
# Call optim using L-BFGS-B. For harder constrains we need constr.Optim
lower <- if (fitted$method == "L-BFGS-B") {
c(rep(0,px),rep(-Inf,pz))[-i]
} else {
-Inf
}
upper <- if (fitted$method == "L-BFGS-B") {
c(rep(Inf,px),rep(Inf,pz))[-i]
} else {
Inf
}
resOthers <- tryCatch(with(fitted$model,
optim(theta.ml[-i], fn = ltildethetaminusi, gr = stildethetaminusi,
method = fitted$method, control = control,
lower = lower, upper = upper)),
warning = function(w) print(w), error = function(e) list(value=NA))
resOthers$value
}
## Initialize
theta.ml <- coef(fitted)
loglik.theta.ml <- c(logLik(fitted))
se <- sqrt(diag(vcov(fitted)))
resProfile <- list()
## Perform profile computations for all requested parameters
cat("Evaluating the profile log-likelihood on a grid ...\n")
for (i in 1:length(profile))
{
cat("i= ",i,"/",length(profile),"\n")
#Index of the parameter in the theta vector
idx <- profile[[i]][1]
#If no borders are given use those from wald intervals (unconstrained)
if (is.na(profile[[i]][2])) profile[[i]][2] <- theta.ml[idx] - 3*se[idx]
if (is.na(profile[[i]][3])) profile[[i]][3] <- theta.ml[idx] + 3*se[idx]
#Evaluate profile loglik on a grid (if requested)
if (profile[[i]][4] > 0) {
thetai.grid <- seq(profile[[i]][2],profile[[i]][3],length.out=profile[[i]][4])
resProfile[[i]] <- matrix(NA, nrow = length(thetai.grid), ncol = 4L,
dimnames = list(NULL, c("grid","profile","estimated","wald")))
for (j in 1:length(thetai.grid)) {
cat("\tj= ",j,"/",length(thetai.grid),"\n")
resProfile[[i]][j,] <- c(thetai.grid[j],
ltildeprofile(thetai.grid[j],idx),
ltildeestim(thetai.grid[j],idx),
#9 June 2009: Bug discovered by L. Held. as part of paper revision. C.f. Pawitan p.63
- 1/2*(1/se[idx]^2)*(thetai.grid[j] - theta.ml[idx])^2)
}
}
}
#9 June 2009. This did not work.
# names(resProfile) <- names(theta.ml)[sapply(profile, function(x) x[4L]) > 0]
names(resProfile) <- names(theta.ml)[sapply(profile, function(x) x[1L])]
## Profile likelihood intervals
ciProfile <- matrix(NA, nrow = length(profile), ncol = 6L,
dimnames = list(NULL, c("idx","hl.low","hl.up","wald.low","wald.up","mle")))
ciProfile[,"idx"] <- sapply(profile, "[", 1L)
ciProfile[,"mle"] <- theta.ml[ciProfile[,"idx"]]
rownames(ciProfile) <- names(theta.ml)[ciProfile[,"idx"]]
if (alpha > 0) {
cat("Computing profile likelihood-based confidence intervals ...\n")
lower <- if (fitted$method == "L-BFGS-B") {
c(rep(0,px),rep(-Inf,pz))
} else {
-Inf
}
for (i in seq_along(profile))
{
cat(i,"/", length(profile),"\n")
#Index of the parameter in the theta vector
idx <- profile[[i]][1]
#Compute highest likelihood intervals
ci.hl <- tryCatch(
likelihood.ci(ltildeprofile, theta.hat = theta.ml[idx],
lower = max(lower[idx], theta.ml[idx]-5*se[idx]),
upper = theta.ml[idx]+5*se[idx], alpha = alpha, i = idx),
warning = function(w) print(w),
error = function(e) rep(NA,2))
#Wald intervals based on expected fisher information
ci.wald <- theta.ml[idx] + c(-1,1) * qnorm(1-alpha/2) * se[idx]
ciProfile[i,2:5] <- c(ci.hl, ci.wald)
}
}
res <- list(lp=resProfile, ci.hl=ciProfile, profileObj=profile)
class(res) <- "profile.twinSIR"
return(res)
}
######################################################################
## Plot the result of the profiler
## Parameters:
## x - the result of calling profile() on a "twinSIR" object
## which - names of selected parameters, NULL meaning all available
## conf.level - level for the horizontal line for -qchisq(,df=1)/2
## legend - logical indicating whether to add a legend to the plot,
## or numeric vector of indexes of plots where to add the legend
######################################################################
plot.profile.twinSIR <- function(x, which = NULL, conf.level = 0.95,
xlab = which, ylab = "normalized log-likelihood",
legend = TRUE, par.settings = list(), ...)
{
## extract relevant components of 'x'
lp <- x$lp[!vapply(X=x$lp, FUN=is.null, FUN.VALUE=FALSE, USE.NAMES=FALSE)]
mle <- x$ci.hl[,"mle"]
## check arguments
which <- if (is.null(which)) {
names(lp)
} else {
match.arg(which, names(lp), several.ok = TRUE)
}
xlab <- rep_len(xlab, length(which))
if (is.logical(legend))
legend <- which(legend)
if (is.list(par.settings)) {
par.defaults <- list(mfrow = sort(n2mfrow(length(which))),
mar = c(5,5,1,1)+.1, las = 1)
par.settings <- modifyList(par.defaults, par.settings)
opar <- do.call("par", par.settings)
on.exit(par(opar))
}
## loop over parameters
for (i in seq_along(which)) {
coefname <- which[i]
matplot(lp[[coefname]][,1L], lp[[coefname]][,-1L], type = "l",
col = 1:3, lty = 1:3, xlab = xlab[i], ylab = ylab)
if (i %in% legend) {
legend(x = "bottomright", legend = c("profile","estimated","Wald"),
col = 1:3, lty = 1:3)
}
## some lines which help interpretation
segments(x0=mle[coefname], y0=par("usr")[3L], y1=0,
lty=2, col="darkgray")
abline(h=-1/2*qchisq(conf.level, df=1), lty=2, col="darkgray")
}
}
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.