plkhci: Profile-likelihood based confidence intervals

Description Usage Arguments Value Note Author(s) See Also Examples

View source: R/plkhci.R

Description

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.

Usage

1
plkhci(x, nlogf, label, prob=0.95, eps=.001, nmax=10, nfcn=0)

Arguments

x

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)

nlogf

the negative log of the density function (not necessarily normalized) for which samples are to be obtained

label

parameter for which confidence bounds are computed

prob

probability associated with the confidence interval

eps

a numerical value. Convergence results when all (logit-transformed) derivatives are smaller eps

nmax

maximum number of Newton-Raphson iterations in each direction

nfcn

number of function calls

Value

2 component vector giving lower and upper p% confidence bounds

Note

At this point, only a single parameter label can be passed to plkhci. This function is part of the Bhat exploration tool

Author(s)

E. Georg Luebeck (FHCRC)

See Also

dfp, newton, logit.hessian

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
        # 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)
        

Example output

Sat May 26 16:37:04 2018 
 
starting at 
1   fmin:  25053.43     10.00 10.00  0.01 

reducing step size 
reducing step size 
reducing step size 
reducing step size 
reducing step size 
reducing step size 
reducing step size 
reducing step size 
reducing step size 
reducing step size 
reducing step size 
reducing step size 
reducing step size 
covariance matrix is not positive-definite or nearly singular 
dfp search intermediate output: 
iter:  0   fmin:  -5722.875    nfcn:  23 
start main loop: 
DFP search completed with status  0 
 
Optimization Result: 
iter:  14   fmin:  -25518.78    nfcn:  178 

  label estimate    low         high      
1 a     19.672695   18.610849   20.779657 
2 b     0.52007907  0.45697481  0.59163367
3 c     0.051330289 0.046352266 0.05628206

 neg. log. likelihood:  -25518.78 
 
 will attempt to compute both bounds (+/- direction) 

trying lower bound ------------------------ 
starting at:    0.9657513 
initial guess:  20.22061 0.4935567 0.05059988 

begin Newton-Raphson search for profile lkh conf. bounds: 
eps value for stop criterium: 0.001 
nmax                        : 10 

 CONVERGENCE:  5  iterations 

chisquare value is:  3.84169 
confidence bound of  a  is  20.77392 
log derivatives:     -0.01737574 -0.006629493 
  label estimate  log deriv   log curv
1 a     20.7739   -55.9705    5563.48 
2 b     0.467573  -0.0173757  2883.67 
3 c     0.0497386 -0.00662949 454.947 

trying upper bound ------------------------ 
starting at:    0.9557458 
initial guess:  19.13606 0.5479867 0.05206013 

begin Newton-Raphson search for profile lkh conf. bounds: 
eps value for stop criterium: 0.001 
nmax                        : 10 

 CONVERGENCE:  5  iterations 

chisquare value is:  3.841336 
confidence bound of  a  is  18.60534 
log derivatives:     -0.019745 -0.003788647 
  label estimate  log deriv   log curv
1 a     18.6053   55.582      5801.87 
2 b     0.576527  -0.019745   3187.37 
3 c     0.0526788 -0.00378865 553.572 

[1] 18.60534 20.77392
Warning messages:
1: In idir * as.vector(dbl0) * sqrt(qchi/(dbb0 - tmp0)) :
  Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

2: In idir * as.vector(dbl0) * sqrt(qchi/(dbb0 - tmp0)) :
  Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

 neg. log. likelihood:  -25518.78 
 
 will attempt to compute both bounds (+/- direction) 

trying lower bound ------------------------ 
starting at:    0.9830178 
initial guess:  19.23707 0.5547335 0.0531601 

begin Newton-Raphson search for profile lkh conf. bounds: 
eps value for stop criterium: 0.001 
nmax                        : 10 

 CONVERGENCE:  5  iterations 

chisquare value is:  3.841255 
confidence bound of  b  is  0.5900483 
log derivatives:     0.001746996 -0.04062658 
  label estimate log deriv  log curv
1 a     18.793   0.001747   5809.89 
2 b     0.590048 -30.365    3205.41 
3 c     0.054662 -0.0406266 587.969 

trying upper bound ------------------------ 
starting at:    0.9399685 
initial guess:  20.11573 0.4875352 0.0494969 

begin Newton-Raphson search for profile lkh conf. bounds: 
eps value for stop criterium: 0.001 
nmax                        : 10 

 CONVERGENCE:  5  iterations 

chisquare value is:  3.841807 
confidence bound of  b  is  0.4555413 
log derivatives:     0.05348663 -0.06622347 
  label estimate  log deriv  log curv
1 a     20.5478   0.0534866  5562.12 
2 b     0.455541  27.253     2870.3  
3 c     0.0472593 -0.0662235 417.611 

[1] 0.4555413 0.5900483
Warning messages:
1: In idir * as.vector(dbl0) * sqrt(qchi/(dbb0 - tmp0)) :
  Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

2: In idir * as.vector(dbl0) * sqrt(qchi/(dbb0 - tmp0)) :
  Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

 neg. log. likelihood:  -25518.78 
 
 will attempt to compute both bounds (+/- direction) 

trying lower bound ------------------------ 
starting at:    1.035521 
initial guess:  19.51417 0.5453683 0.05381505 

begin Newton-Raphson search for profile lkh conf. bounds: 
eps value for stop criterium: 0.001 
nmax                        : 10 

 CONVERGENCE:  5  iterations 

chisquare value is:  3.84129 
confidence bound of  c  is  0.0559518 
log derivatives:     0.1181582 0.07858222 
  label estimate  log deriv log curv
1 a     19.3912   0.118158  5724.99 
2 b     0.568441  0.0785822 3075.5  
3 c     0.0559518 -22.2201  589.781 

trying upper bound ------------------------ 
starting at:    0.8904929 
initial guess:  19.83219 0.4959326 0.04883894 

begin Newton-Raphson search for profile lkh conf. bounds: 
eps value for stop criterium: 0.001 
nmax                        : 10 

 CONVERGENCE:  5  iterations 

chisquare value is:  3.841794 
confidence bound of  c  is  0.04592346 
log derivatives:     0.2032383 0.1572934 
  label estimate  log deriv log curv
1 a     20.0279   0.203238  5635.08 
2 b     0.46983   0.157293  2983.97 
3 c     0.0459235 15.9238   406.846 

[1] 0.04592346 0.05595180
Warning messages:
1: In idir * as.vector(dbl0) * sqrt(qchi/(dbb0 - tmp0)) :
  Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

2: In idir * as.vector(dbl0) * sqrt(qchi/(dbb0 - tmp0)) :
  Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

 neg. log. likelihood:  -25518.78 
 
 will attempt to compute both bounds (+/- direction) 

trying lower bound ------------------------ 
starting at:    0.7204909 
initial guess:  19.53959 0.5412233 0.05341652 

begin Newton-Raphson search for profile lkh conf. bounds: 
eps value for stop criterium: 0.001 
nmax                        : 10 

 CONVERGENCE:  5  iterations 

chisquare value is:  2.705379 
confidence bound of  c  is  0.05525405 
log derivatives:     0.0874495 0.05734045 
  label estimate  log deriv log curv
1 a     19.433    0.0874495 5719.05 
2 b     0.560745  0.0573404 3070.03 
3 c     0.0552541 -18.3131  577.082 

trying upper bound ------------------------ 
starting at:    0.6347411 
initial guess:  19.80648 0.4997397 0.04923941 

begin Newton-Raphson search for profile lkh conf. bounds: 
eps value for stop criterium: 0.001 
nmax                        : 10 

 CONVERGENCE:  5  iterations 

chisquare value is:  2.705787 
confidence bound of  c  is  0.04685451 
log derivatives:     0.1343002 0.1012921 
  label estimate  log deriv log curv
1 a     19.9672   0.1343    5643.64 
2 b     0.477939  0.101292  2993.28 
3 c     0.0468545 13.7563   423.373 

[1] 0.04685451 0.05525405
Warning messages:
1: In idir * as.vector(dbl0) * sqrt(qchi/(dbb0 - tmp0)) :
  Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

2: In idir * as.vector(dbl0) * sqrt(qchi/(dbb0 - tmp0)) :
  Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

Bhat documentation built on May 29, 2017, 8:14 p.m.

Related to plkhci in Bhat...