[R] confint() in stats4 package

Jinsong Zhao jszhao at yeah.net
Wed Aug 3 16:18:41 CEST 2011


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



More information about the R-help mailing list