R: Profile-likelihood based confidence intervals

Question

I am using the function plkhci from library Bhat to construct Profile-likelihood based confidence intervals and I got this warning:

Warning message: In dqstep(list(label = x$label, est = btrf(xt, x$low, x$upp), low = x$low, : oops: unable to find stepsize, use default

when i run

r <- dfp(x,f=nlogf)

Can I ignore this warning as I still can get the output?

Following is the complete coding:

library(Bhat)
beta0<--8
beta1<-0.03
gamma<-0.0105
alpha<-0.05
n<-100
u<-runif(n)
u
x<-rnorm(n)
x
c<-rexp(100,1/1515)
c
t1<-(1/gamma)*log(1-((gamma/(exp(beta0+beta1*x)))*(log(1-u))))
t1
t<-pmin(t1,c)
t
delta<-1*(t1>c)
delta
length(delta)
cp<-length(delta[delta==1])/n
cp
delta[delta==1]<-ifelse(rbinom(length(delta[delta==1]),1,0.5),1,2)
delta
deltae<-ifelse(delta==0, 1,0)
deltar<-ifelse(delta==1, 1,0)
deltai<-ifelse(delta==2, 1,0)
dat=data.frame(t,delta, deltae,deltar,deltai,x)
dat$interval[delta==2] <- as.character(cut(dat$t[delta==2], breaks=seq(0, 600, 100)))
labs <- cut(dat$t[delta==2], breaks=seq(0, 600, 100))
dat$lower[delta==2]<-as.numeric( sub("\\((.+),.*", "\\1", labs) )
dat$upper[delta==2]<-as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs) )  
data0<-dat[which(dat$delta==0),]#uncensored data 
data1<-dat[which(dat$delta==1),]#right censored data 
data2<-dat[which(dat$delta==2),]#interval censored data 
nlogf<-function(para)
{
  b0<-para[1]
  b1<-para[2]
  g<-para[3]
  e<-sum((b0+b1*data0$x)+g*data0$t+(1/g)*exp(b0+b1*data0$x)*(1-exp(g*data0$t)))
  r<-sum((1/g)*exp(b0+b1*data1$x)*(1-exp(g*data1$t)))
  i<-sum(log(exp((1/g)*exp(b0+b1*data2$x)*(1-exp(g*data2$lower)))-exp((1/g)*exp(b0+b1*data2$x)*(1-exp(g*data2$upper)))))
  l<-e+r+i
  return(-l)
}
x <- list(label=c("beta0","beta1","gamma"),est=c(-8,0.03,0.0105),low=c(-10,0,0),upp=c(10,1,1))
r <- dfp(x,f=nlogf)
x$est <- r$est
plkhci(x,nlogf,"beta0")
plkhci(x,nlogf,"beta1")
plkhci(x,nlogf,"gamma")

Answer 1

I am giving you a super long answer, but it will help you see that you can chase down your own error messages (most of the time, sometimes this means of looking at functions will not work). It is good to see what is happening inside a method when it throws an warning because sometimes it is fine and sometimes you need to fix your data.

This function is REALLY involved! You can look at it by typing dfp into the R command line (NO TRAILING PARENTHESES) and it will print out the whole function.

17 lines from the end, you will see an assignment:

del <- dqstep(x, f, sens = 0.01)

You can see that this calls the function dqstep, which is reflected in your warning.

You can see this function by typing dqstep into the command line of R again. In reading through this function, also long but not so tedious, there is this section of boolean logic:

if (r < 0 | is.na(r) | b == 0) {
        warning("oops: unable to find stepsize, use default")
        cat("problem with ", x$label[i], "\n")
        break
    }  

This is the culprit, it returns the message you are getting. The line right above it spells out how r is calculated. You are feeding this function your default x from the prior function plus a sensitivity equations (which I assume dfp generates, it is huge and ugly, so I did not untangle all of it). When the previous nested function returns either an r value lower than Zero, and r value of NA or a b value of ZERO, that message is displayed.

The second error tells you that it was likely b==0 because b is in the denominator and it returned and infinity value, so NO STEP SIZE IS RETURNED FROM THIS NESTED FUNCTION to the variable del in dfp.

The step is fed into THIS equation:

h <- logit.hessian(x, f, del, dapprox = FALSE, nfcn)  

which you can look into by typing logit.hessian into the R commandline.

When you do, you see that del is a step size in a logit scale, with a default value of del=rep(0.002, length(x$est))...which the function set for you because running the function dqstep returned no value.

So, you now get to decide if using that step size in the calculation of your confidence interval seems right or if there is a problem with your data which needs resolving to make this work better for you.

When I ran it, line by line, I got this message:

Error in if (denom <= 0) { : missing value where TRUE/FALSE needed

at this line of code: r <- dfp(x,f=nlogf(x))

Which makes me think I was correct.

That is how I chase down issues I have with messages from packages when I get a message like yours.