R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Semiparametric Elicitation of a bounded parameter

SEL	R Documentation

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

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)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(SEL)
Loading required package: splines
Loading required package: quadprog
Loading required package: lattice
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/SEL/SEL.Rd_%03d_medium.png", width=480, height=480)
> ### Name: SEL
> ### Title: Semiparametric Elicitation of a bounded parameter
> ### Aliases: SEL
> ### Keywords: misc
> 
> ### ** 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: 
> ##D ### generate some random elicitations
> ##D numb <- max(rnbinom(1, mu = 4, size = 1), 1)
> ##D x0 <- runif(1, -1000, 1000)
> ##D x1 <- x0+runif(1, 0, 1000)
> ##D xs <- sort(runif(numb, x0, x1))
> ##D y <- runif(numb+1)
> ##D ys <- (cumsum(y)/sum(y))[1:numb]
> ##D fit1 <- SEL(xs, ys, bounds = c(x0, x1))
> ##D fit2 <- SEL(xs, ys, d = 1, inknts = xs, bounds = c(x0, x1))
> ##D fit3 <- SEL(xs, ys, d = 15, N = 0, bounds = c(x0, x1))
> ##D comparePlot(fit1, fit2, fit3, type="cdf")
> ## End(Not run)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>