SEL {SEL}

Semiparametric Elicitation of a bounded parameter

Description

Fits a distribution to quantiles (stated for example by an expert) using B-splines with a Brier entropy penalty.

There exists a print and a summary method to display details of the fitted SEL object. The fitted density (or cdf) can be displayed with the plot method. Different SEL objects can be displayed in one plot with the comparePlot function. The fitted density (or cdf) can be evaluated with the predict method. The coef method extracts the coefficients of the fitted B-spline basis. The quantSEL function calculates quantiles of the fitted distribution and the rvSEL function generates random variables from a fitted SEL object.

Usage

  SEL(x, alpha, bounds = c(0, 1), d = 4, inknts = x, N, gamma, Delta,
      fitbnds = c(1e-8, 10)*diff(bounds), pistar = NULL,
      constr = c("none", "unimodal", "decreasing", "increasing"),
      mode)

Arguments

x

Numeric vector of quantiles.

alpha

Numeric vector determining the levels of the quantiles specified in x.

bounds

Vector containing the bounds for the expert's density.

N

Number of equally spaced inner knots (ignored if inknts is specified).

d

Degree of the spline (the order of the spline is d+1), values larger than 20 are to be avoided; they can lead to numerical instabilities.

inknts

Vector specifying the inner knots. If equal to NULL, the B-spline basis reduces to the Bernstein polynomial basis. Per default the inner knots are located on the values specified via x.

gamma

Parameter controlling the weight of the negative entropy in the objective function to be optimized (see Bornkamp and Ickstadt (2009)).

Delta

The code calculates the gamma parameter achieving a certain overall L2 distance between the specified quantiles and the fitted distribution function (only if gamma is not specified).

fitbnds

Numeric of length 2 for the bisection search algorithm searching for gamma giving the error Delta (only relevant if gamma is missing).

pistar

Prior density function used in Brier divergence (see Bornkamp and Ickstadt (2009). If NULL Brier entropy is used.

constr

Character vector specifying, which shape constraint should be used, should be one of "none", "unimodal", "decreasing", "increasing".

mode

Numerical value needed when constraint = "unimodal". The code selects the knot average closest to mode as the location of the mode for the finite dimensional constraints on the coefficients of the B-spline in the optimization (note that the final density will only have a mode close to the value specified via mode).

Value

An SEL object containing the following entries

constr

Character determining the shape constraint used.

inknts

Inner knots.

nord

Order of the spline.

bounds

Bounds for the elicited quantity.

xalpha

List containing the specified quantiles.

Omega

Penalty matrix used for calculating Brier entropy.

coefs

Fitted coefficients of the B-spline basis.

pistar

Function handed over to SEL via the pistar argument (if any).

dplus

The d vector, only necessary when pistar is specified (see Bornkamp and Ickstadt (2009).)

Author(s)

Bjoern Bornkamp

References

Bornkamp, B. and Ickstadt, K. (2009). A Note on B-Splines for Semiparametric Elicitation. The American Statistician, 63, 373–377

O'Hagan A., Buck C. E., Daneshkhah, A., Eiser, R., Garthwaite, P., Jenkinson, D., Oakley, J. and Rakow, T. (2006), Uncertain Judgements: Eliciting Expert Probabilities, John Wiley and Sons Inc.

Dierckx, P. (1993), Curve and Surface Fitting with Splines, Clarendon Press

See Also

plot.SEL, comparePlot, quantSEL, rvSEL

Examples

### example from O'Hagan et al. (2006)
### 1st example in Bornkamp and Ickstadt (2009)
x <- c(177.5, 183.75, 190, 205, 220)
y <- c(0.175, 0.33, 0.5, 0.75, 0.95)
I <- SEL(x, y, d = 1, Delta = 0.015, bounds = c(165, 250), inknts = x)
II <- SEL(x, y, d = 14, N = 0, Delta = 0.015, bounds = c(165, 250))
III <- SEL(x, y, d = 4, Delta = 0.015, bounds = c(165, 250), inknts = x)
IV <- SEL(x, y, d = 4, Delta = 0.015, bounds = c(165, 250), inknts = x,
      constr = "u", mode = 185)
comparePlot(I, II, III, IV, type = "cdf")
comparePlot(I, II, III, IV, type = "density")



### bimodal example 
x2 <- c(0.1, 0.2, 0.5, 0.8, 0.9)
y2 <- c(0.2, 0.4, 0.45, 0.85, 0.99)
fit1 <- SEL(x2, y2, Delta=0.05, d = 4, inknts = x2)
fit2 <- SEL(x2, y2, Delta=0.05, d = 15, N = 0)
comparePlot(fit1, fit2, superpose = TRUE)



### sample from SEL object and evaluate density
xxx <- rvSEL(50000, fit1)
hist(xxx, breaks=100, freq=FALSE)
curve(predict(fit1, newdata=x), add=TRUE)



### illustrate shrinkage against uniform dist.
gma01 <- SEL(x2, y2, gamma = 0.1)
gma02 <- SEL(x2, y2, gamma = 0.2)
gma04 <- SEL(x2, y2, gamma = 0.4)
gma10 <- SEL(x2, y2, gamma = 1.0)
comparePlot(gma01, gma02, gma04, gma10, superpose = TRUE)



### including shape constraints
x <- c(177.5, 183.75, 190, 205, 220)
y <- c(0.175, 0.33, 0.5, 0.75, 0.95)
unconstr1 <- SEL(x, y, Delta = 0.05, bounds = c(165, 250))
unconstr2 <- SEL(x, y, d = 10, N = 0, Delta = 0.05, bounds = c(165,250))
unimod1 <- SEL(x, y, Delta = 0.05, bounds = c(165, 250),
            constr = "unimodal", mode = 185)
unimod2 <- SEL(x, y, d = 10, N = 0, Delta = 0.05, bounds = c(165, 250),
           constr =  "unimodal", mode = 185)
comparePlot(unconstr1, unconstr2, unimod1, unimod2)



### shrinkage against another distribution
pr <- function(x) dbeta(x, 2, 2)
pr01  <- SEL(x2, y2, gamma = 0.1, d = 3, pistar = pr)
pr03  <- SEL(x2, y2, gamma = 0.3, d = 3, pistar = pr)
pr12 <- SEL(x2, y2, gamma = 1.2, pistar = pr)
comparePlot(pr01, pr03, pr12, superpose = TRUE)



### 2nd example from Bornkamp and Ickstadt (2009)
# theta
# "true" density
pmixbeta <- function(x, a1, b1, a2, b2){
  0.3*pbeta(x/20,a1,b1)+0.7*pbeta(x/20,a2,b2)
}
dmixbeta <- function(x, a1, b1, a2, b2){
  out <- 0.3*dbeta(x/20,a1,b1)+0.7*dbeta(x/20,a2,b2)
  out/20
}
x <- c(2,10,15)
a1 <- 1;a2 <- 10;b1 <- 15;b2 <- 5
mu <- pmixbeta(x, a1, b1, a2, b2)
set.seed(1) # simulate experts statements
y <- rnorm(length(mu), mu, sqrt(mu*(1-mu)*0.05^2))
thet <- SEL(x, y, d = 4, Delta = 0.03, bounds = c(0, 20), inknts = x)
plot(thet, ylim = c(0,0.25))
curve(dmixbeta(x, a1, b1, a2, b2), add=TRUE, col="red")
abline(h=0.05, lty = 2)
legend("topright", c("true density", "elicited density", "uniform density"), 
       col=c(2,1,1), lty=c(1,1,2))

# sigma
# "true" density
dtriang <- function(x,m,a,b){
  inds <- x < m;res <- numeric(length(x))
  res[inds] <- 2*(x[inds]-a)/((b-a)*(m-a))
  res[!inds] <- 2*(b-x[!inds])/((b-a)*(b-m))
  res
}
ptriang <- function(x,m,a,b){
  inds <- x < m;res <- numeric(length(x))
  res[inds] <- (x[inds]-a)^2/((b-a)*(m-a))
  res[!inds] <- 1-(b-x[!inds])^2/((b-a)*(b-m))
  res
}
x <- c(1,2,4)
mu <- ptriang(x, 1, 0, 5)
set.seed(1) # simulate experts statements
y <- rnorm(length(mu), mu, sqrt(mu*(1-mu)*0.05^2))
sig <- SEL(x, y, d = 4, Delta = 0.03, bounds = c(0, 5),
       inknts = x, mode = 1, constr="unimodal")
plot(sig, ylim=c(0,0.4))
curve(dtriang(x, 1, 0, 5), add=TRUE, col="red")
abline(h=0.2, lty = 2)
legend("topright", c("true density", "elicited density", "uniform density"), 
       col=c(2,1,1), lty=c(1,1,2))


## Not run: 
### generate some random elicitations
numb <- max(rnbinom(1, mu = 4, size = 1), 1)
x0 <- runif(1, -1000, 1000)
x1 <- x0+runif(1, 0, 1000)
xs <- sort(runif(numb, x0, x1))
y <- runif(numb+1)
ys <- (cumsum(y)/sum(y))[1:numb]
fit1 <- SEL(xs, ys, bounds = c(x0, x1))
fit2 <- SEL(xs, ys, d = 1, inknts = xs, bounds = c(x0, x1))
fit3 <- SEL(xs, ys, d = 15, N = 0, bounds = c(x0, x1))
comparePlot(fit1, fit2, fit3, type="cdf")

## End(Not run)