function to find prob
*100% confidence intervals using
profile-likelihood. Numerical solutions are obtained via a modified
Newton-Raphson algorithm. The method is described in Venzon and Moolgavkar,
Journal of the Royal Statistical Society, Series C vol 37, no.1, 1988, pp.
87-94.
plkhci(x, nlogf, label, prob = 0.95, eps = 0.001, nmax = 10, nfcn = 0)
a list with components 'label' (of mode character), 'est' (the parameter vector with the initial guess), 'low' (vector with lower bounds), and 'upp' (vector with upper bounds)
the negative log of the density function (not necessarily normalized) for which samples are to be obtained
parameter for which confidence bounds are computed
probability associated with the confidence interval
a numerical value. Convergence results when all
(logit-transformed) derivatives are smaller eps
maximum number of Newton-Raphson iterations in each direction
number of function calls
2 component vector giving lower and upper p% confidence bounds
# NOT RUN {
# generate some Poisson counts on the fly
dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50))
data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05))))
# neg. log-likelihood of Poisson model with 'linear-quadratic' mean:
nlogf <- function (x) {
ds <- data[, 1]
y <- data[, 2]
g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds))
return(sum(g - y * log(g)))
}
# for example define
x <- list(label=c("a","b","c"),est=c(10.,10.,.01),low=c(0,0,0),upp=c(100,20,.1))
# get MLEs using dfp:
r <- dfp(x,f=nlogf)
x$est <- r$est
plkhci(x,nlogf,"a")
plkhci(x,nlogf,"b")
plkhci(x,nlogf,"c")
# e.g. 90% confidence bounds for "c"
plkhci(x,nlogf,"c",prob=0.9)
# }
Run the code above in your browser using DataLab