Hi there,
I had a problem when I hoped to get confidence intervals for the
parameters I got using mle() of stats4 package. This problem would not
appear if ``fixed'' option was not used. The following mini-example will
demo the problem:
x <- c(100, 56, 32, 18, 10, 1)
r <- c(18, 17, 10, 6, 4, 3)
n <- c(18, 22, 17, 21, 23, 20)
loglik.1 <- function(alpha, beta, c) {
x <- log10(x)
P <- c + (1-c) * pnorm(alpha + beta * x)
control <- which(x == -Inf)
if (length(control) != 0) P[control] <- c
P <- pmax(pmin(P,1),0)
-(sum(r * log(P)) + sum((n - r)* log(1-P)))
}
loglik.2 <- function(alpha, beta) {
x <- log10(x)
P <- pnorm(alpha + beta * x)
P <- pmax(pmin(P,1),0)
-(sum(r * log(P)) + sum((n - r)* log(1-P)))
}
library(stats4)
fit.1 <- mle(loglik.1, start = list(alpha = 0, beta = 0, c = 0), method
= "BFGS", fixed = list(c=0))
fit.2 <- mle(loglik.2, start = list(alpha = 0, beta = 0), method =
"BFGS", fixed = list())
> confint(fit.1)
Profiling...
Error in approx(sp$y, sp$x, xout = cutoff) :
need at least two non-NA values to interpolate
In addition: Warning message:
In approx(sp$y, sp$x, xout = cutoff) : collapsing to unique 'x' values
> confint(fit.2)
Profiling...
2.5 % 97.5 %
alpha -2.5187909 -1.144600
beta 0.9052395 1.876322
The version I test the above code is 2.11.1 and 2.13.1.
I hope to know what's the matter? and how to avoid the error, and get
the correct confidence intervals for the parameters? Any suggestions
will be really appreciated.
P.S.: I noticed that there was a file named mle.R.rej in the source
directory of stats4. A broken patch? Thanks!
Regards,
Jinsong
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
