twinstim_profile: Profile Likelihood Computation and Confidence Intervals for...

In surveillance: Temporal and Spatio-Temporal Modeling and Monitoring of Epidemic Phenomena

Profile Likelihood Computation and Confidence Intervals for twinstim objects

Description

Function to compute estimated and profile likelihood based confidence intervals for twinstim objects. Computations might be cumbersome!

WARNING: the implementation is not well tested, simply uses optim (ignoring optimizer settings from the original fit), and does not return the complete set of coefficients at each grid point.

Usage

## S3 method for class 'twinstim'
profile(fitted, profile, alpha = 0.05,
        control = list(fnscale = -1, maxit = 100, trace = 1),
        do.ltildeprofile=FALSE, ...)

Arguments

fitted	an object of class "twinstim".
profile	a list with elements being numeric vectors of length 4. These vectors must have the form c(index, lower, upper, gridsize). index: index of the parameter to be profiled in the vector coef(fitted). lower, upper: lower/upper limit of the grid on which the profile log-likelihood is evaluated. Can also be NA in which case lower/upper equals the lower/upper bound of the respective 0.3 % Wald confidence interval (+-3*se). gridsize: grid size of the equally spaced grid between lower and upper. Can also be 0 in which case the profile log-likelihood for this parameter is not evaluated on a grid.
alpha	(1-\alpha)\% profile likelihood based confidence intervals are computed. If alpha <= 0, then no confidence intervals are computed. This is currently not implemented.
control	control object to use in optim for the profile log-likelihood computations. It might be necessary to control maxit or reltol in order to obtain results in finite time.
do.ltildeprofile	If TRUE calculate profile likelihood as well. This might take a while, since an optimisation for all other parameters has to be performed. Useful for likelihood based confidence intervals. Default: FALSE.
...	unused (argument of the generic).

Value

list with profile log-likelihood evaluations on the grid, and – not implemented yet – highest likelihood and Wald confidence intervals. The argument profile is also returned.

Author(s)

Michael Höhle

Examples

# profiling takes a while
## Not run: 
#Load the twinstim model fitted to the IMD data
data("imdepi", "imdepifit")
# for profiling we need the model environment
imdepifit <- update(imdepifit, model=TRUE)

#Generate profiling object for a list of parameters for the new model
names <- c("h.(Intercept)","e.typeC")
coefList <- lapply(names, function(name) {
  c(pmatch(name,names(coef(imdepifit))),NA,NA,11)
})

#Profile object (necessary to specify a more loose convergence
#criterion). Speed things up by using do.ltildeprofile=FALSE (the default)
prof <- profile(imdepifit, coefList,
  control=list(reltol=0.1, REPORT=1), do.ltildeprofile=TRUE)

#Plot result for one variable
par(mfrow=c(1,2))
for (name in names) {
  with(as.data.frame(prof$lp[[name]]),
       matplot(grid,cbind(profile,estimated,wald),
               type="l",xlab=name,ylab="loglik"))
  legend(x="bottomleft",c("profile","estimated","wald"),lty=1:3,col=1:3)
}

## End(Not run)

surveillance documentation built on Nov. 4, 2024, 3 a.m.