BiomedEngX/CircosHeatmap-aardio
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 3 Wiki Activity
CircosHeatmap-aardio/dist/lib/r-library/boot/tests/Examples/boot-Ex.Rout.save

2009 lines
74 KiB
Plaintext
Raw Permalink Normal View History

首次上传
2025-01-12 00:52:51 +08:00
R Under development (unstable) (2019-12-19 r77606) -- "Unsuffered Consequences"
Copyright (C) 2019 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.
Natural language support but running in an English locale
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.
> pkgname <- "boot"
> source(file.path(R.home("share"), "R", "examples-header.R"))
> options(warn = 1)
> library('boot')
>
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
> cleanEx()
> nameEx("Imp.Estimates")
> ### * Imp.Estimates
>
> flush(stderr()); flush(stdout())
>
> ### Name: Imp.Estimates
> ### Title: Importance Sampling Estimates
> ### Aliases: Imp.Estimates imp.moments imp.prob imp.quantile imp.reg
> ### Keywords: htest nonparametric
>
> ### ** Examples
>
> # Example 9.8 of Davison and Hinkley (1997) requires tilting the
> # resampling distribution of the studentized statistic to be centred
> # at the observed value of the test statistic, 1.84. In this example
> # we show how certain estimates can be found using resamples taken from
> # the tilted distribution.
> grav1 <- gravity[as.numeric(gravity[,2]) >= 7, ]
> grav.fun <- function(dat, w, orig) {
+ strata <- tapply(dat[, 2], as.numeric(dat[, 2]))
+ d <- dat[, 1]
+ ns <- tabulate(strata)
+ w <- w/tapply(w, strata, sum)[strata]
+ mns <- as.vector(tapply(d * w, strata, sum)) # drop names
+ mn2 <- tapply(d * d * w, strata, sum)
+ s2hat <- sum((mn2 - mns^2)/ns)
+ c(mns[2] - mns[1], s2hat, (mns[2] - mns[1] - orig)/sqrt(s2hat))
+ }
> grav.z0 <- grav.fun(grav1, rep(1, 26), 0)
> grav.L <- empinf(data = grav1, statistic = grav.fun, stype = "w",
+ strata = grav1[,2], index = 3, orig = grav.z0[1])
> grav.tilt <- exp.tilt(grav.L, grav.z0[3], strata = grav1[, 2])
> grav.tilt.boot <- boot(grav1, grav.fun, R = 199, stype = "w",
+ strata = grav1[, 2], weights = grav.tilt$p,
+ orig = grav.z0[1])
> # Since the weights are needed for all calculations, we shall calculate
> # them once only.
> grav.w <- imp.weights(grav.tilt.boot)
> grav.mom <- imp.moments(grav.tilt.boot, w = grav.w, index = 3)
> grav.p <- imp.prob(grav.tilt.boot, w = grav.w, index = 3, t0 = grav.z0[3])
> unlist(grav.p)
t0 raw rat reg
1.8401182 1.0000000 0.9778246 0.9767292
> grav.q <- imp.quantile(grav.tilt.boot, w = grav.w, index = 3,
+ alpha = c(0.9, 0.95, 0.975, 0.99))
> as.data.frame(grav.q)
alpha raw rat reg
1 0.900 3.048484 1.170707 1.210631
2 0.950 3.237935 1.565928 1.576607
3 0.975 3.629448 1.895876 1.923657
4 0.990 3.629448 2.258157 2.295230
>
>
>
> cleanEx()
> nameEx("abc.ci")
> ### * abc.ci
>
> flush(stderr()); flush(stdout())
>
> ### Name: abc.ci
> ### Title: Nonparametric ABC Confidence Intervals
> ### Aliases: abc.ci
> ### Keywords: nonparametric htest
>
> ### ** Examples
>
> ## Don't show:
> op <- options(digits = 5)
> ## End(Don't show)
> # 90% and 95% confidence intervals for the correlation
> # coefficient between the columns of the bigcity data
>
> abc.ci(bigcity, corr, conf=c(0.90,0.95))
conf
[1,] 0.90 0.95815 0.99173
[2,] 0.95 0.94937 0.99307
>
> # A 95% confidence interval for the difference between the means of
> # the last two samples in gravity
> mean.diff <- function(y, w)
+ { gp1 <- 1:table(as.numeric(y$series))[1]
+ sum(y[gp1, 1] * w[gp1]) - sum(y[-gp1, 1] * w[-gp1])
+ }
> grav1 <- gravity[as.numeric(gravity[, 2]) >= 7, ]
> ## IGNORE_RDIFF_BEGIN
> abc.ci(grav1, mean.diff, strata = grav1$series)
[1] 0.95000 -6.70758 -0.39394
> ## IGNORE_RDIFF_END
> ## Don't show:
> options(op)
> ## End(Don't show)
>
>
>
> cleanEx()
> nameEx("boot")
> ### * boot
>
> flush(stderr()); flush(stdout())
>
> ### Name: boot
> ### Title: Bootstrap Resampling
> ### Aliases: boot boot.return c.boot
> ### Keywords: nonparametric htest
>
> ### ** Examples
>
> ## Don't show:
> op <- options(digits = 5)
> ## End(Don't show)
> # Usual bootstrap of the ratio of means using the city data
> ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
> boot(city, ratio, R = 999, stype = "w")
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = city, statistic = ratio, R = 999, stype = "w")
Bootstrap Statistics :
original bias std. error
t1* 1.5203 0.049024 0.21583
>
>
> # Stratified resampling for the difference of means. In this
> # example we will look at the difference of means between the final
> # two series in the gravity data.
> diff.means <- function(d, f)
+ { n <- nrow(d)
+ gp1 <- 1:table(as.numeric(d$series))[1]
+ m1 <- sum(d[gp1,1] * f[gp1])/sum(f[gp1])
+ m2 <- sum(d[-gp1,1] * f[-gp1])/sum(f[-gp1])
+ ss1 <- sum(d[gp1,1]^2 * f[gp1]) - (m1 * m1 * sum(f[gp1]))
+ ss2 <- sum(d[-gp1,1]^2 * f[-gp1]) - (m2 * m2 * sum(f[-gp1]))
+ c(m1 - m2, (ss1 + ss2)/(sum(f) - 2))
+ }
> grav1 <- gravity[as.numeric(gravity[,2]) >= 7,]
> boot(grav1, diff.means, R = 999, stype = "f", strata = grav1[,2])
STRATIFIED BOOTSTRAP
Call:
boot(data = grav1, statistic = diff.means, R = 999, stype = "f",
strata = grav1[, 2])
Bootstrap Statistics :
original bias std. error
t1* -2.8462 0.007469 1.5528
t2* 16.8462 -1.475719 6.7147
>
> # In this example we show the use of boot in a prediction from
> # regression based on the nuclear data. This example is taken
> # from Example 6.8 of Davison and Hinkley (1997). Notice also
> # that two extra arguments to 'statistic' are passed through boot.
> nuke <- nuclear[, c(1, 2, 5, 7, 8, 10, 11)]
> nuke.lm <- glm(log(cost) ~ date+log(cap)+ne+ct+log(cum.n)+pt, data = nuke)
> nuke.diag <- glm.diag(nuke.lm)
> nuke.res <- nuke.diag$res * nuke.diag$sd
> nuke.res <- nuke.res - mean(nuke.res)
>
> # We set up a new data frame with the data, the standardized
> # residuals and the fitted values for use in the bootstrap.
> nuke.data <- data.frame(nuke, resid = nuke.res, fit = fitted(nuke.lm))
>
> # Now we want a prediction of plant number 32 but at date 73.00
> new.data <- data.frame(cost = 1, date = 73.00, cap = 886, ne = 0,
+ ct = 0, cum.n = 11, pt = 1)
> new.fit <- predict(nuke.lm, new.data)
>
> nuke.fun <- function(dat, inds, i.pred, fit.pred, x.pred)
+ {
+ lm.b <- glm(fit+resid[inds] ~ date+log(cap)+ne+ct+log(cum.n)+pt,
+ data = dat)
+ pred.b <- predict(lm.b, x.pred)
+ c(coef(lm.b), pred.b - (fit.pred + dat$resid[i.pred]))
+ }
>
> nuke.boot <- boot(nuke.data, nuke.fun, R = 999, m = 1,
+ fit.pred = new.fit, x.pred = new.data)
> # The bootstrap prediction squared error would then be found by
> mean(nuke.boot$t[, 8]^2)
[1] 0.084223
> # Basic bootstrap prediction limits would be
> new.fit - sort(nuke.boot$t[, 8])[c(975, 25)]
[1] 6.1488 7.2619
>
>
> # Finally a parametric bootstrap. For this example we shall look
> # at the air-conditioning data. In this example our aim is to test
> # the hypothesis that the true value of the index is 1 (i.e. that
> # the data come from an exponential distribution) against the
> # alternative that the data come from a gamma distribution with
> # index not equal to 1.
> air.fun <- function(data) {
+ ybar <- mean(data$hours)
+ para <- c(log(ybar), mean(log(data$hours)))
+ ll <- function(k) {
+ if (k <= 0) 1e200 else lgamma(k)-k*(log(k)-1-para[1]+para[2])
+ }
+ khat <- nlm(ll, ybar^2/var(data$hours))$estimate
+ c(ybar, khat)
+ }
>
> air.rg <- function(data, mle) {
+ # Function to generate random exponential variates.
+ # mle will contain the mean of the original data
+ out <- data
+ out$hours <- rexp(nrow(out), 1/mle)
+ out
+ }
>
> air.boot <- boot(aircondit, air.fun, R = 999, sim = "parametric",
+ ran.gen = air.rg, mle = mean(aircondit$hours))
>
> # The bootstrap p-value can then be approximated by
> sum(abs(air.boot$t[,2]-1) > abs(air.boot$t0[2]-1))/(1+air.boot$R)
[1] 0.446
> ## Don't show:
> options(op)
> ## End(Don't show)
>
>
>
> cleanEx()
> nameEx("boot.array")
> ### * boot.array
>
> flush(stderr()); flush(stdout())
>
> ### Name: boot.array
> ### Title: Bootstrap Resampling Arrays
> ### Aliases: boot.array
> ### Keywords: nonparametric
>
> ### ** Examples
>
> # A frequency array for a nonparametric bootstrap
> city.boot <- boot(city, corr, R = 40, stype = "w")
> boot.array(city.boot)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 0 0 0 2 1 1 3 1
[2,] 2 2 0 2 0 0 1 0 1 2
[3,] 0 1 1 1 2 1 2 1 0 1
[4,] 1 0 1 0 0 1 2 1 1 3
[5,] 0 3 0 1 0 0 2 0 2 2
[6,] 0 0 1 1 0 1 4 1 2 0
[7,] 1 3 1 2 0 0 0 1 2 0
[8,] 0 2 2 2 0 1 0 1 0 2
[9,] 3 0 1 0 2 2 1 1 0 0
[10,] 1 1 0 1 3 0 0 1 2 1
[11,] 0 1 2 1 3 1 0 0 0 2
[12,] 1 1 3 1 0 0 0 2 1 1
[13,] 0 1 0 1 1 3 1 1 1 1
[14,] 2 0 0 1 1 1 0 1 3 1
[15,] 1 0 0 0 0 0 5 0 4 0
[16,] 0 0 1 1 2 1 0 1 4 0
[17,] 1 0 1 0 5 0 0 1 0 2
[18,] 1 0 1 0 3 2 1 1 0 1
[19,] 1 2 2 2 0 1 0 1 1 0
[20,] 1 1 1 1 2 0 1 0 1 2
[21,] 2 0 2 1 1 1 0 1 0 2
[22,] 2 0 1 1 1 0 1 2 1 1
[23,] 1 2 1 1 0 1 1 1 0 2
[24,] 3 0 1 0 1 0 1 1 0 3
[25,] 1 2 1 2 1 1 0 1 1 0
[26,] 2 0 2 0 1 3 0 1 0 1
[27,] 2 1 0 2 1 1 0 1 2 0
[28,] 1 1 2 0 1 0 1 0 3 1
[29,] 1 2 1 1 1 3 1 0 0 0
[30,] 2 0 0 2 1 0 0 2 1 2
[31,] 0 1 1 1 2 0 3 1 0 1
[32,] 0 0 1 1 1 4 0 0 0 3
[33,] 1 0 3 2 0 0 0 2 1 1
[34,] 0 1 0 2 4 1 1 1 0 0
[35,] 0 1 1 0 2 1 1 0 3 1
[36,] 3 0 1 0 2 0 1 1 2 0
[37,] 2 1 1 0 1 0 2 2 1 0
[38,] 1 0 2 0 1 1 3 1 1 0
[39,] 1 1 0 2 0 1 1 1 2 1
[40,] 0 1 2 3 1 1 0 2 0 0
>
> perm.cor <- function(d,i) cor(d$x,d$u[i])
> city.perm <- boot(city, perm.cor, R = 40, sim = "permutation")
> boot.array(city.perm, indices = TRUE)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 7 8 2 6 9 5 4 3 10 1
[2,] 1 9 8 6 5 3 4 7 2 10
[3,] 10 1 6 9 4 5 2 8 3 7
[4,] 4 5 7 6 3 10 9 8 2 1
[5,] 10 3 6 7 5 2 8 9 1 4
[6,] 4 10 8 2 5 6 7 3 9 1
[7,] 5 8 10 3 2 4 7 1 6 9
[8,] 3 8 4 9 2 7 1 5 10 6
[9,] 6 5 1 3 7 2 4 10 9 8
[10,] 4 9 6 10 8 1 7 3 5 2
[11,] 5 10 1 8 6 2 3 4 9 7
[12,] 7 6 1 3 10 2 8 4 5 9
[13,] 10 8 3 6 9 7 1 5 2 4
[14,] 6 7 10 1 5 4 8 2 9 3
[15,] 9 1 5 3 2 8 6 10 7 4
[16,] 9 5 7 6 4 10 1 3 8 2
[17,] 1 8 3 10 4 5 9 7 6 2
[18,] 1 10 5 8 2 4 7 9 6 3
[19,] 8 2 7 4 9 3 10 6 5 1
[20,] 10 4 8 6 5 9 1 3 2 7
[21,] 2 1 10 4 5 7 8 3 6 9
[22,] 3 1 10 5 2 4 8 6 9 7
[23,] 8 4 1 10 9 6 2 7 5 3
[24,] 2 3 7 8 1 6 4 5 10 9
[25,] 4 8 9 5 10 6 2 7 3 1
[26,] 1 9 2 10 6 5 4 8 7 3
[27,] 9 10 7 1 3 8 6 2 4 5
[28,] 6 7 2 3 4 1 8 9 10 5
[29,] 5 2 1 9 4 8 7 6 10 3
[30,] 7 4 6 2 5 1 9 10 3 8
[31,] 8 5 1 6 7 9 4 3 2 10
[32,] 3 10 4 5 7 6 1 9 8 2
[33,] 7 4 9 10 3 5 6 1 2 8
[34,] 10 6 2 8 7 3 4 9 1 5
[35,] 7 5 6 1 10 9 2 4 3 8
[36,] 9 1 8 6 10 5 2 7 3 4
[37,] 5 6 8 2 4 10 3 9 7 1
[38,] 5 10 4 3 1 2 8 9 7 6
[39,] 6 9 7 3 5 1 10 2 8 4
[40,] 9 3 8 2 6 7 5 1 10 4
>
>
>
> cleanEx()
> nameEx("boot.ci")
> ### * boot.ci
>
> flush(stderr()); flush(stdout())
>
> ### Name: boot.ci
> ### Title: Nonparametric Bootstrap Confidence Intervals
> ### Aliases: boot.ci
> ### Keywords: nonparametric htest
>
> ### ** Examples
>
> # confidence intervals for the city data
> ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
> city.boot <- boot(city, ratio, R = 999, stype = "w", sim = "ordinary")
> boot.ci(city.boot, conf = c(0.90, 0.95),
+ type = c("norm", "basic", "perc", "bca"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates
CALL :
boot.ci(boot.out = city.boot, conf = c(0.9, 0.95), type = c("norm",
"basic", "perc", "bca"))
Intervals :
Level Normal Basic
90% ( 1.116, 1.826 ) ( 1.054, 1.742 )
95% ( 1.048, 1.894 ) ( 0.956, 1.775 )
Level Percentile BCa
90% ( 1.299, 1.987 ) ( 1.292, 1.957 )
95% ( 1.266, 2.085 ) ( 1.254, 2.058 )
Calculations and Intervals on Original Scale
>
> # studentized confidence interval for the two sample
> # difference of means problem using the final two series
> # of the gravity data.
> diff.means <- function(d, f)
+ { n <- nrow(d)
+ gp1 <- 1:table(as.numeric(d$series))[1]
+ m1 <- sum(d[gp1,1] * f[gp1])/sum(f[gp1])
+ m2 <- sum(d[-gp1,1] * f[-gp1])/sum(f[-gp1])
+ ss1 <- sum(d[gp1,1]^2 * f[gp1]) - (m1 * m1 * sum(f[gp1]))
+ ss2 <- sum(d[-gp1,1]^2 * f[-gp1]) - (m2 * m2 * sum(f[-gp1]))
+ c(m1 - m2, (ss1 + ss2)/(sum(f) - 2))
+ }
> grav1 <- gravity[as.numeric(gravity[,2]) >= 7, ]
> grav1.boot <- boot(grav1, diff.means, R = 999, stype = "f",
+ strata = grav1[ ,2])
> boot.ci(grav1.boot, type = c("stud", "norm"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates
CALL :
boot.ci(boot.out = grav1.boot, type = c("stud", "norm"))
Intervals :
Level Normal Studentized
95% (-5.897, 0.190 ) (-7.370, -0.052 )
Calculations and Intervals on Original Scale
>
> # Nonparametric confidence intervals for mean failure time
> # of the air-conditioning data as in Example 5.4 of Davison
> # and Hinkley (1997)
> mean.fun <- function(d, i)
+ { m <- mean(d$hours[i])
+ n <- length(i)
+ v <- (n-1)*var(d$hours[i])/n^2
+ c(m, v)
+ }
> air.boot <- boot(aircondit, mean.fun, R = 999)
> boot.ci(air.boot, type = c("norm", "basic", "perc", "stud"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates
CALL :
boot.ci(boot.out = air.boot, type = c("norm", "basic", "perc",
"stud"))
Intervals :
Level Normal Basic
95% ( 32.2, 184.5 ) ( 21.2, 169.4 )
Level Studentized Percentile
95% ( 45.6, 295.8 ) ( 46.8, 195.0 )
Calculations and Intervals on Original Scale
>
> # Now using the log transformation
> # There are two ways of doing this and they both give the
> # same intervals.
>
> # Method 1
> boot.ci(air.boot, type = c("norm", "basic", "perc", "stud"),
+ h = log, hdot = function(x) 1/x)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates
CALL :
boot.ci(boot.out = air.boot, type = c("norm", "basic", "perc",
"stud"), h = log, hdot = function(x) 1/x)
Intervals :
Level Normal Basic
95% ( 4.015, 5.491 ) ( 4.093, 5.521 )
Level Studentized Percentile
95% ( 3.922, 5.808 ) ( 3.845, 5.273 )
Calculations and Intervals on Transformed Scale
>
> # Method 2
> vt0 <- air.boot$t0[2]/air.boot$t0[1]^2
> vt <- air.boot$t[, 2]/air.boot$t[ ,1]^2
> boot.ci(air.boot, type = c("norm", "basic", "perc", "stud"),
+ t0 = log(air.boot$t0[1]), t = log(air.boot$t[,1]),
+ var.t0 = vt0, var.t = vt)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates
CALL :
boot.ci(boot.out = air.boot, type = c("norm", "basic", "perc",
"stud"), var.t0 = vt0, var.t = vt, t0 = log(air.boot$t0[1]),
t = log(air.boot$t[, 1]))
Intervals :
Level Normal Basic
95% ( 4.070, 5.435 ) ( 4.093, 5.521 )
Level Studentized Percentile
95% ( 3.922, 5.808 ) ( 3.845, 5.273 )
Calculations and Intervals on Original Scale
>
>
>
> cleanEx()
> nameEx("censboot")
> ### * censboot
>
> flush(stderr()); flush(stdout())
>
> ### Name: censboot
> ### Title: Bootstrap for Censored Data
> ### Aliases: censboot cens.return
> ### Keywords: survival
>
> ### ** Examples
>
> library(survival)
Attaching package: ‘survival’
The following object is masked from ‘package:boot’:
aml
> # Example 3.9 of Davison and Hinkley (1997) does a bootstrap on some
> # remission times for patients with a type of leukaemia. The patients
> # were divided into those who received maintenance chemotherapy and
> # those who did not. Here we are interested in the median remission
> # time for the two groups.
> data(aml, package = "boot") # not the version in survival.
> aml.fun <- function(data) {
+ surv <- survfit(Surv(time, cens) ~ group, data = data)
+ out <- NULL
+ st <- 1
+ for (s in 1:length(surv$strata)) {
+ inds <- st:(st + surv$strata[s]-1)
+ md <- min(surv$time[inds[1-surv$surv[inds] >= 0.5]])
+ st <- st + surv$strata[s]
+ out <- c(out, md)
+ }
+ out
+ }
> aml.case <- censboot(aml, aml.fun, R = 499, strata = aml$group)
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
>
> # Now we will look at the same statistic using the conditional
> # bootstrap and the weird bootstrap. For the conditional bootstrap
> # the survival distribution is stratified but the censoring
> # distribution is not.
>
> aml.s1 <- survfit(Surv(time, cens) ~ group, data = aml)
> aml.s2 <- survfit(Surv(time-0.001*cens, 1-cens) ~ 1, data = aml)
> aml.cond <- censboot(aml, aml.fun, R = 499, strata = aml$group,
+ F.surv = aml.s1, G.surv = aml.s2, sim = "cond")
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
>
>
> # For the weird bootstrap we must redefine our function slightly since
> # the data will not contain the group number.
> aml.fun1 <- function(data, str) {
+ surv <- survfit(Surv(data[, 1], data[, 2]) ~ str)
+ out <- NULL
+ st <- 1
+ for (s in 1:length(surv$strata)) {
+ inds <- st:(st + surv$strata[s] - 1)
+ md <- min(surv$time[inds[1-surv$surv[inds] >= 0.5]])
+ st <- st + surv$strata[s]
+ out <- c(out, md)
+ }
+ out
+ }
> aml.wei <- censboot(cbind(aml$time, aml$cens), aml.fun1, R = 499,
+ strata = aml$group, F.surv = aml.s1, sim = "weird")
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
Warning in min(surv$time[inds[1 - surv$surv[inds] >= 0.5]]) :
no non-missing arguments to min; returning Inf
>
> # Now for an example where a cox regression model has been fitted
> # the data we will look at the melanoma data of Example 7.6 from
> # Davison and Hinkley (1997). The fitted model assumes that there
> # is a different survival distribution for the ulcerated and
> # non-ulcerated groups but that the thickness of the tumour has a
> # common effect. We will also assume that the censoring distribution
> # is different in different age groups. The statistic of interest
> # is the linear predictor. This is returned as the values at a
> # number of equally spaced points in the range of interest.
> data(melanoma, package = "boot")
> library(splines)# for ns
> mel.cox <- coxph(Surv(time, status == 1) ~ ns(thickness, df=4) + strata(ulcer),
+ data = melanoma)
> mel.surv <- survfit(mel.cox)
> agec <- cut(melanoma$age, c(0, 39, 49, 59, 69, 100))
> mel.cens <- survfit(Surv(time - 0.001*(status == 1), status != 1) ~
+ strata(agec), data = melanoma)
> mel.fun <- function(d) {
+ t1 <- ns(d$thickness, df=4)
+ cox <- coxph(Surv(d$time, d$status == 1) ~ t1+strata(d$ulcer))
+ ind <- !duplicated(d$thickness)
+ u <- d$thickness[!ind]
+ eta <- cox$linear.predictors[!ind]
+ sp <- smooth.spline(u, eta, df=20)
+ th <- seq(from = 0.25, to = 10, by = 0.25)
+ predict(sp, th)$y
+ }
> mel.str <- cbind(melanoma$ulcer, agec)
>
> # this is slow!
> mel.mod <- censboot(melanoma, mel.fun, R = 499, F.surv = mel.surv,
+ G.surv = mel.cens, cox = mel.cox, strata = mel.str, sim = "model")
> # To plot the original predictor and a 95% pointwise envelope for it
> mel.env <- envelope(mel.mod)$point
> th <- seq(0.25, 10, by = 0.25)
> plot(th, mel.env[1, ], ylim = c(-2, 2),
+ xlab = "thickness (mm)", ylab = "linear predictor", type = "n")
> lines(th, mel.mod$t0, lty = 1)
> matlines(th, t(mel.env), lty = 2)
>
>
>
> cleanEx()
detaching ‘package:splines’, ‘package:survival’
> nameEx("control")
> ### * control
>
> flush(stderr()); flush(stdout())
>
> ### Name: control
> ### Title: Control Variate Calculations
> ### Aliases: control
> ### Keywords: nonparametric
>
> ### ** Examples
>
> # Use of control variates for the variance of the air-conditioning data
> mean.fun <- function(d, i)
+ { m <- mean(d$hours[i])
+ n <- nrow(d)
+ v <- (n-1)*var(d$hours[i])/n^2
+ c(m, v)
+ }
> air.boot <- boot(aircondit, mean.fun, R = 999)
> control(air.boot, index = 2, bias.adj = TRUE)
[1] -11.25369
> air.cont <- control(air.boot, index = 2)
> # Now let us try the variance on the log scale.
> air.cont1 <- control(air.boot, t0 = log(air.boot$t0[2]),
+ t = log(air.boot$t[, 2]))
>
>
>
> cleanEx()
> nameEx("cv.glm")
> ### * cv.glm
>
> flush(stderr()); flush(stdout())
>
> ### Name: cv.glm
> ### Title: Cross-validation for Generalized Linear Models
> ### Aliases: cv.glm
> ### Keywords: regression
>
> ### ** Examples
>
> # leave-one-out and 6-fold cross-validation prediction error for
> # the mammals data set.
> data(mammals, package="MASS")
> mammals.glm <- glm(log(brain) ~ log(body), data = mammals)
> (cv.err <- cv.glm(mammals, mammals.glm)$delta)
[1] 0.4918650 0.4916571
> (cv.err.6 <- cv.glm(mammals, mammals.glm, K = 6)$delta)
[1] 0.4911574 0.4887828
>
> # As this is a linear model we could calculate the leave-one-out
> # cross-validation estimate without any extra model-fitting.
> muhat <- fitted(mammals.glm)
> mammals.diag <- glm.diag(mammals.glm)
> (cv.err <- mean((mammals.glm$y - muhat)^2/(1 - mammals.diag$h)^2))
[1] 0.491865
>
>
> # leave-one-out and 11-fold cross-validation prediction error for
> # the nodal data set. Since the response is a binary variable an
> # appropriate cost function is
> cost <- function(r, pi = 0) mean(abs(r-pi) > 0.5)
>
> nodal.glm <- glm(r ~ stage+xray+acid, binomial, data = nodal)
> (cv.err <- cv.glm(nodal, nodal.glm, cost, K = nrow(nodal))$delta)
[1] 0.1886792 0.1886792
> (cv.11.err <- cv.glm(nodal, nodal.glm, cost, K = 11)$delta)
[1] 0.2075472 0.2057672
>
>
>
> cleanEx()
> nameEx("empinf")
> ### * empinf
>
> flush(stderr()); flush(stdout())
>
> ### Name: empinf
> ### Title: Empirical Influence Values
> ### Aliases: empinf
> ### Keywords: nonparametric math
>
> ### ** Examples
>
> # The empirical influence values for the ratio of means in
> # the city data.
> ratio <- function(d, w) sum(d$x *w)/sum(d$u*w)
> empinf(data = city, statistic = ratio)
[1] -1.04367815 -0.58417763 -0.37092459 -0.18958996 0.03164142 0.10544878
[7] 0.09236345 0.20365074 1.02178280 0.73381132
> city.boot <- boot(city, ratio, 499, stype="w")
> empinf(boot.out = city.boot, type = "reg")
1 1 1 1 1 1
-1.188052737 -0.702204565 -0.447209091 -0.378938644 0.007498044 0.162261064
1 1 1 1
0.116272627 0.242543405 1.166357773 1.021472124
>
> # A statistic that may be of interest in the difference of means
> # problem is the t-statistic for testing equality of means. In
> # the bootstrap we get replicates of the difference of means and
> # the variance of that statistic and then want to use this output
> # to get the empirical influence values of the t-statistic.
> grav1 <- gravity[as.numeric(gravity[,2]) >= 7,]
> grav.fun <- function(dat, w) {
+ strata <- tapply(dat[, 2], as.numeric(dat[, 2]))
+ d <- dat[, 1]
+ ns <- tabulate(strata)
+ w <- w/tapply(w, strata, sum)[strata]
+ mns <- as.vector(tapply(d * w, strata, sum)) # drop names
+ mn2 <- tapply(d * d * w, strata, sum)
+ s2hat <- sum((mn2 - mns^2)/ns)
+ c(mns[2] - mns[1], s2hat)
+ }
>
> grav.boot <- boot(grav1, grav.fun, R = 499, stype = "w",
+ strata = grav1[, 2])
>
> # Since the statistic of interest is a function of the bootstrap
> # statistics, we must calculate the bootstrap replicates and pass
> # them to empinf using the t argument.
> grav.z <- (grav.boot$t[,1]-grav.boot$t0[1])/sqrt(grav.boot$t[,2])
> empinf(boot.out = grav.boot, t = grav.z)
1 1 1 1 1 1 1
-2.7187022 -0.9515832 -2.4507504 -1.3409041 0.6167772 -1.2768844 -1.3438236
1 1 1 1 1 1 2
-0.2313133 -1.0504975 -3.1949223 1.2809344 3.7306567 8.9310126 2.7195701
2 2 2 2 2 2 2
4.1111581 3.1632367 1.1507436 -2.2930450 -2.7667452 -2.3645554 -0.2002981
2 2 2 2 2
2.0526618 0.3615873 -1.9421623 -2.1031282 -1.8890234
>
>
>
> cleanEx()
> nameEx("envelope")
> ### * envelope
>
> flush(stderr()); flush(stdout())
>
> ### Name: envelope
> ### Title: Confidence Envelopes for Curves
> ### Aliases: envelope
> ### Keywords: dplot htest
>
> ### ** Examples
>
> # Testing whether the final series of measurements of the gravity data
> # may come from a normal distribution. This is done in Examples 4.7
> # and 4.8 of Davison and Hinkley (1997).
> grav1 <- gravity$g[gravity$series == 8]
> grav.z <- (grav1 - mean(grav1))/sqrt(var(grav1))
> grav.gen <- function(dat, mle) rnorm(length(dat))
> grav.qqboot <- boot(grav.z, sort, R = 999, sim = "parametric",
+ ran.gen = grav.gen)
> grav.qq <- qqnorm(grav.z, plot.it = FALSE)
> grav.qq <- lapply(grav.qq, sort)
> plot(grav.qq, ylim = c(-3.5, 3.5), ylab = "Studentized Order Statistics",
+ xlab = "Normal Quantiles")
> grav.env <- envelope(grav.qqboot, level = 0.9)
> lines(grav.qq$x, grav.env$point[1, ], lty = 4)
> lines(grav.qq$x, grav.env$point[2, ], lty = 4)
> lines(grav.qq$x, grav.env$overall[1, ], lty = 1)
> lines(grav.qq$x, grav.env$overall[2, ], lty = 1)
>
>
>
> cleanEx()
> nameEx("exp.tilt")
> ### * exp.tilt
>
> flush(stderr()); flush(stdout())
>
> ### Name: exp.tilt
> ### Title: Exponential Tilting
> ### Aliases: exp.tilt
> ### Keywords: nonparametric smooth
>
> ### ** Examples
>
> # Example 9.8 of Davison and Hinkley (1997) requires tilting the resampling
> # distribution of the studentized statistic to be centred at the observed
> # value of the test statistic 1.84. This can be achieved as follows.
> grav1 <- gravity[as.numeric(gravity[,2]) >=7 , ]
> grav.fun <- function(dat, w, orig) {
+ strata <- tapply(dat[, 2], as.numeric(dat[, 2]))
+ d <- dat[, 1]
+ ns <- tabulate(strata)
+ w <- w/tapply(w, strata, sum)[strata]
+ mns <- as.vector(tapply(d * w, strata, sum)) # drop names
+ mn2 <- tapply(d * d * w, strata, sum)
+ s2hat <- sum((mn2 - mns^2)/ns)
+ c(mns[2]-mns[1], s2hat, (mns[2]-mns[1]-orig)/sqrt(s2hat))
+ }
> grav.z0 <- grav.fun(grav1, rep(1, 26), 0)
> grav.L <- empinf(data = grav1, statistic = grav.fun, stype = "w",
+ strata = grav1[,2], index = 3, orig = grav.z0[1])
> grav.tilt <- exp.tilt(grav.L, grav.z0[3], strata = grav1[,2])
> boot(grav1, grav.fun, R = 499, stype = "w", weights = grav.tilt$p,
+ strata = grav1[,2], orig = grav.z0[1])
STRATIFIED WEIGHTED BOOTSTRAP
Call:
boot(data = grav1, statistic = grav.fun, R = 499, stype = "w",
strata = grav1[, 2], weights = grav.tilt$p, orig = grav.z0[1])
Bootstrap Statistics :
original bias std. error mean(t*)
t1* 2.846154 -0.3661063 1.705171 5.702944
t2* 2.392353 -0.3538294 1.002889 3.444050
t3* 0.000000 -0.5160619 1.314298 1.473456
>
>
>
> cleanEx()
> nameEx("glm.diag.plots")
> ### * glm.diag.plots
>
> flush(stderr()); flush(stdout())
>
> ### Name: glm.diag.plots
> ### Title: Diagnostics plots for generalized linear models
> ### Aliases: glm.diag.plots
> ### Keywords: regression dplot hplot
>
> ### ** Examples
>
> # In this example we look at the leukaemia data which was looked at in
> # Example 7.1 of Davison and Hinkley (1997)
> data(leuk, package = "MASS")
> leuk.mod <- glm(time ~ ag-1+log10(wbc), family = Gamma(log), data = leuk)
> leuk.diag <- glm.diag(leuk.mod)
> glm.diag.plots(leuk.mod, leuk.diag)
>
>
>
> cleanEx()
> nameEx("jack.after.boot")
> ### * jack.after.boot
>
> flush(stderr()); flush(stdout())
>
> ### Name: jack.after.boot
> ### Title: Jackknife-after-Bootstrap Plots
> ### Aliases: jack.after.boot
> ### Keywords: hplot nonparametric
>
> ### ** Examples
>
> # To draw the jackknife-after-bootstrap plot for the head size data as in
> # Example 3.24 of Davison and Hinkley (1997)
> frets.fun <- function(data, i) {
+ pcorr <- function(x) {
+ # Function to find the correlations and partial correlations between
+ # the four measurements.
+ v <- cor(x)
+ v.d <- diag(var(x))
+ iv <- solve(v)
+ iv.d <- sqrt(diag(iv))
+ iv <- - diag(1/iv.d) %*% iv %*% diag(1/iv.d)
+ q <- NULL
+ n <- nrow(v)
+ for (i in 1:(n-1))
+ q <- rbind( q, c(v[i, 1:i], iv[i,(i+1):n]) )
+ q <- rbind( q, v[n, ] )
+ diag(q) <- round(diag(q))
+ q
+ }
+ d <- data[i, ]
+ v <- pcorr(d)
+ c(v[1,], v[2,], v[3,], v[4,])
+ }
> frets.boot <- boot(log(as.matrix(frets)), frets.fun, R = 999)
> # we will concentrate on the partial correlation between head breadth
> # for the first son and head length for the second. This is the 7th
> # element in the output of frets.fun so we set index = 7
> jack.after.boot(frets.boot, useJ = FALSE, stinf = FALSE, index = 7)
>
>
>
> cleanEx()
> nameEx("k3.linear")
> ### * k3.linear
>
> flush(stderr()); flush(stdout())
>
> ### Name: k3.linear
> ### Title: Linear Skewness Estimate
> ### Aliases: k3.linear
> ### Keywords: nonparametric
>
> ### ** Examples
>
> # To estimate the skewness of the ratio of means for the city data.
> ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
> k3.linear(empinf(data = city, statistic = ratio))
[1] 7.831452e-05
>
>
>
> cleanEx()
> nameEx("linear.approx")
> ### * linear.approx
>
> flush(stderr()); flush(stdout())
>
> ### Name: linear.approx
> ### Title: Linear Approximation of Bootstrap Replicates
> ### Aliases: linear.approx
> ### Keywords: nonparametric
>
> ### ** Examples
>
> # Using the city data let us look at the linear approximation to the
> # ratio statistic and its logarithm. We compare these with the
> # corresponding plots for the bigcity data
>
> ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
> city.boot <- boot(city, ratio, R = 499, stype = "w")
> bigcity.boot <- boot(bigcity, ratio, R = 499, stype = "w")
> op <- par(pty = "s", mfrow = c(2, 2))
>
> # The first plot is for the city data ratio statistic.
> city.lin1 <- linear.approx(city.boot)
> lim <- range(c(city.boot$t,city.lin1))
> plot(city.boot$t, city.lin1, xlim = lim, ylim = lim,
+ main = "Ratio; n=10", xlab = "t*", ylab = "tL*")
> abline(0, 1)
>
> # Now for the log of the ratio statistic for the city data.
> city.lin2 <- linear.approx(city.boot,t0 = log(city.boot$t0),
+ t = log(city.boot$t))
> lim <- range(c(log(city.boot$t),city.lin2))
> plot(log(city.boot$t), city.lin2, xlim = lim, ylim = lim,
+ main = "Log(Ratio); n=10", xlab = "t*", ylab = "tL*")
> abline(0, 1)
>
> # The ratio statistic for the bigcity data.
> bigcity.lin1 <- linear.approx(bigcity.boot)
> lim <- range(c(bigcity.boot$t,bigcity.lin1))
> plot(bigcity.lin1, bigcity.boot$t, xlim = lim, ylim = lim,
+ main = "Ratio; n=49", xlab = "t*", ylab = "tL*")
> abline(0, 1)
>
> # Finally the log of the ratio statistic for the bigcity data.
> bigcity.lin2 <- linear.approx(bigcity.boot,t0 = log(bigcity.boot$t0),
+ t = log(bigcity.boot$t))
> lim <- range(c(log(bigcity.boot$t),bigcity.lin2))
> plot(bigcity.lin2, log(bigcity.boot$t), xlim = lim, ylim = lim,
+ main = "Log(Ratio); n=49", xlab = "t*", ylab = "tL*")
> abline(0, 1)
>
> par(op)
>
>
>
> graphics::par(get("par.postscript", pos = 'CheckExEnv'))
> cleanEx()
> nameEx("lines.saddle.distn")
> ### * lines.saddle.distn
>
> flush(stderr()); flush(stdout())
>
> ### Name: lines.saddle.distn
> ### Title: Add a Saddlepoint Approximation to a Plot
> ### Aliases: lines.saddle.distn
> ### Keywords: aplot smooth nonparametric
>
> ### ** Examples
>
> # In this example we show how a plot such as that in Figure 9.9 of
> # Davison and Hinkley (1997) may be produced. Note the large number of
> # bootstrap replicates required in this example.
> expdata <- rexp(12)
> vfun <- function(d, i) {
+ n <- length(d)
+ (n-1)/n*var(d[i])
+ }
> exp.boot <- boot(expdata,vfun, R = 9999)
> exp.L <- (expdata - mean(expdata))^2 - exp.boot$t0
> exp.tL <- linear.approx(exp.boot, L = exp.L)
> hist(exp.tL, nclass = 50, probability = TRUE)
> exp.t0 <- c(0, sqrt(var(exp.boot$t)))
> exp.sp <- saddle.distn(A = exp.L/12,wdist = "m", t0 = exp.t0)
>
> # The saddlepoint approximation in this case is to the density of
> # t-t0 and so t0 must be added for the plot.
> lines(exp.sp, h = function(u, t0) u+t0, J = function(u, t0) 1,
+ t0 = exp.boot$t0)
>
>
>
> cleanEx()
> nameEx("norm.ci")
> ### * norm.ci
>
> flush(stderr()); flush(stdout())
>
> ### Name: norm.ci
> ### Title: Normal Approximation Confidence Intervals
> ### Aliases: norm.ci
> ### Keywords: htest
>
> ### ** Examples
>
> # In Example 5.1 of Davison and Hinkley (1997), normal approximation
> # confidence intervals are found for the air-conditioning data.
> air.mean <- mean(aircondit$hours)
> air.n <- nrow(aircondit)
> air.v <- air.mean^2/air.n
> norm.ci(t0 = air.mean, var.t0 = air.v)
conf
[1,] 0.95 46.93055 169.2361
> exp(norm.ci(t0 = log(air.mean), var.t0 = 1/air.n)[2:3])
[1] 61.38157 190.31782
>
> # Now a more complicated example - the ratio estimate for the city data.
> ratio <- function(d, w)
+ sum(d$x * w)/sum(d$u *w)
> city.v <- var.linear(empinf(data = city, statistic = ratio))
> norm.ci(t0 = ratio(city,rep(0.1,10)), var.t0 = city.v)
conf
[1,] 0.95 1.167046 1.873579
>
>
>
> cleanEx()
> nameEx("plot.boot")
> ### * plot.boot
>
> flush(stderr()); flush(stdout())
>
> ### Name: plot.boot
> ### Title: Plots of the Output of a Bootstrap Simulation
> ### Aliases: plot.boot
> ### Keywords: hplot nonparametric
>
> ### ** Examples
>
> # We fit an exponential model to the air-conditioning data and use
> # that for a parametric bootstrap. Then we look at plots of the
> # resampled means.
> air.rg <- function(data, mle) rexp(length(data), 1/mle)
>
> air.boot <- boot(aircondit$hours, mean, R = 999, sim = "parametric",
+ ran.gen = air.rg, mle = mean(aircondit$hours))
> plot(air.boot)
>
> # In the difference of means example for the last two series of the
> # gravity data
> grav1 <- gravity[as.numeric(gravity[, 2]) >= 7, ]
> grav.fun <- function(dat, w) {
+ strata <- tapply(dat[, 2], as.numeric(dat[, 2]))
+ d <- dat[, 1]
+ ns <- tabulate(strata)
+ w <- w/tapply(w, strata, sum)[strata]
+ mns <- as.vector(tapply(d * w, strata, sum)) # drop names
+ mn2 <- tapply(d * d * w, strata, sum)
+ s2hat <- sum((mn2 - mns^2)/ns)
+ c(mns[2] - mns[1], s2hat)
+ }
>
> grav.boot <- boot(grav1, grav.fun, R = 499, stype = "w", strata = grav1[, 2])
> plot(grav.boot)
> # now suppose we want to look at the studentized differences.
> grav.z <- (grav.boot$t[, 1]-grav.boot$t0[1])/sqrt(grav.boot$t[, 2])
> plot(grav.boot, t = grav.z, t0 = 0)
>
> # In this example we look at the one of the partial correlations for the
> # head dimensions in the dataset frets.
> frets.fun <- function(data, i) {
+ pcorr <- function(x) {
+ # Function to find the correlations and partial correlations between
+ # the four measurements.
+ v <- cor(x)
+ v.d <- diag(var(x))
+ iv <- solve(v)
+ iv.d <- sqrt(diag(iv))
+ iv <- - diag(1/iv.d) %*% iv %*% diag(1/iv.d)
+ q <- NULL
+ n <- nrow(v)
+ for (i in 1:(n-1))
+ q <- rbind( q, c(v[i, 1:i], iv[i,(i+1):n]) )
+ q <- rbind( q, v[n, ] )
+ diag(q) <- round(diag(q))
+ q
+ }
+ d <- data[i, ]
+ v <- pcorr(d)
+ c(v[1,], v[2,], v[3,], v[4,])
+ }
> frets.boot <- boot(log(as.matrix(frets)), frets.fun, R = 999)
> plot(frets.boot, index = 7, jack = TRUE, stinf = FALSE, useJ = FALSE)
>
>
>
> cleanEx()
> nameEx("saddle")
> ### * saddle
>
> flush(stderr()); flush(stdout())
>
> ### Name: saddle
> ### Title: Saddlepoint Approximations for Bootstrap Statistics
> ### Aliases: saddle
> ### Keywords: smooth nonparametric
>
> ### ** Examples
>
> # To evaluate the bootstrap distribution of the mean failure time of
> # air-conditioning equipment at 80 hours
> saddle(A = aircondit$hours/12, u = 80)
$spa
pdf cdf
0.01005866 0.24446677
$zeta.hat
[1] -0.02580078
>
> # Alternatively this can be done using a conditional poisson
> saddle(A = cbind(aircondit$hours/12,1), u = c(80, 12),
+ wdist = "p", type = "cond")
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.090909
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.909091
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.090909
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.909091
$spa
pdf cdf
0.01005943 0.24438736
$zeta.hat
A1 A2
-0.02580805 0.89261577
$zeta2.hat
[1] 0.6931472
>
> # To use the Lugananni-Rice approximation to this
> saddle(A = cbind(aircondit$hours/12,1), u = c(80, 12),
+ wdist = "p", type = "cond",
+ LR = TRUE)
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.090909
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.909091
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.090909
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.909091
$spa
pdf cdf
0.01005943 0.24447362
$zeta.hat
A1 A2
-0.02580805 0.89261577
$zeta2.hat
[1] 0.6931472
>
> # Example 9.16 of Davison and Hinkley (1997) calculates saddlepoint
> # approximations to the distribution of the ratio statistic for the
> # city data. Since the statistic is not in itself a linear combination
> # of random Variables, its distribution cannot be found directly.
> # Instead the statistic is expressed as the solution to a linear
> # estimating equation and hence its distribution can be found. We
> # get the saddlepoint approximation to the pdf and cdf evaluated at
> # t = 1.25 as follows.
> jacobian <- function(dat,t,zeta)
+ {
+ p <- exp(zeta*(dat$x-t*dat$u))
+ abs(sum(dat$u*p)/sum(p))
+ }
> city.sp1 <- saddle(A = city$x-1.25*city$u, u = 0)
> city.sp1$spa[1] <- jacobian(city, 1.25, city.sp1$zeta.hat) * city.sp1$spa[1]
> city.sp1
$spa
pdf cdf
0.05565040 0.02436306
$zeta.hat
[1] -0.02435547
>
>
>
> cleanEx()
> nameEx("saddle.distn")
> ### * saddle.distn
>
> flush(stderr()); flush(stdout())
>
> ### Name: saddle.distn
> ### Title: Saddlepoint Distribution Approximations for Bootstrap Statistics
> ### Aliases: saddle.distn
> ### Keywords: nonparametric smooth dplot
>
> ### ** Examples
>
> # The bootstrap distribution of the mean of the air-conditioning
> # failure data: fails to find value on R (and probably on S too)
> air.t0 <- c(mean(aircondit$hours), sqrt(var(aircondit$hours)/12))
> ## Not run: saddle.distn(A = aircondit$hours/12, t0 = air.t0)
>
> # alternatively using the conditional poisson
> saddle.distn(A = cbind(aircondit$hours/12, 1), u = 12, wdist = "p",
+ type = "cond", t0 = air.t0)
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.344718
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.655282
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.344718
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.655282
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.832240
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.167760
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.832240
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.167760
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.444538
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.555462
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.444538
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.555462
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.494449
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.505551
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.494449
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.505551
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.594269
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.405731
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.594269
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.405731
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.544359
Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.455641
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.544359
Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.455641
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.519404
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.480596
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.519404
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.480596
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.561394
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.438606
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.561394
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.438606
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.290549
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.709451
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.290549
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.709451
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.019703
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.980297
Warning in dpois(y, mu, log = TRUE) : non-integer x = 11.019703
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.980297
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.748857
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.251143
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.748857
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.251143
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.478011
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.521989
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.478011
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.521989
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.207165
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.792835
Warning in dpois(y, mu, log = TRUE) : non-integer x = 10.207165
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.792835
Warning in dpois(y, mu, log = TRUE) : non-integer x = 9.936320
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.063680
Warning in dpois(y, mu, log = TRUE) : non-integer x = 9.936320
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.063680
Warning in dpois(y, mu, log = TRUE) : non-integer x = 9.665474
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.334526
Warning in dpois(y, mu, log = TRUE) : non-integer x = 9.665474
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.334526
Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.785225
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.214775
Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.785225
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.214775
Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.175822
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.824178
Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.175822
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.824178
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.566419
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.433581
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.566419
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.433581
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.957016
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.042984
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.957016
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.042984
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.347613
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.652387
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.347613
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.652387
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.738210
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.261790
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.738210
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.261790
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.128807
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.871193
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.128807
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.871193
Saddlepoint Distribution Approximations
Call :
saddle.distn(A = cbind(aircondit$hours/12, 1), u = 12, wdist = "p",
type = "cond", t0 = air.t0)
Quantiles of the Distribution
0.1% 27.4
0.5% 35.4
1.0% 39.7
2.5% 46.7
5.0% 53.5
10.0% 62.5
20.0% 75.3
50.0% 104.5
80.0% 139.0
90.0% 158.8
95.0% 175.9
97.5% 191.2
99.0% 209.6
99.5% 222.4
99.9% 249.5
Smoothing spline used 20 points in the range 9.8 to 304.7.
>
> # Distribution of the ratio of a sample of size 10 from the bigcity
> # data, taken from Example 9.16 of Davison and Hinkley (1997).
> ratio <- function(d, w) sum(d$x *w)/sum(d$u * w)
> city.v <- var.linear(empinf(data = city, statistic = ratio))
> bigcity.t0 <- c(mean(bigcity$x)/mean(bigcity$u), sqrt(city.v))
> Afn <- function(t, data) cbind(data$x - t*data$u, 1)
> ufn <- function(t, data) c(0,10)
> saddle.distn(A = Afn, u = ufn, wdist = "b", type = "cond",
+ t0 = bigcity.t0, data = bigcity)
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Warning in eval(family$initialize) :
non-integer counts in a binomial glm!
Saddlepoint Distribution Approximations
Call :
saddle.distn(A = Afn, u = ufn, wdist = "b", type = "cond", t0 = bigcity.t0,
data = bigcity)
Quantiles of the Distribution
0.1% 1.070
0.5% 1.092
1.0% 1.104
2.5% 1.122
5.0% 1.139
10.0% 1.158
20.0% 1.184
50.0% 1.237
80.0% 1.304
90.0% 1.348
95.0% 1.392
97.5% 1.436
99.0% 1.494
99.5% 1.537
99.9% 1.636
Smoothing spline used 20 points in the range 1.014 to 1.96.
>
> # From Example 9.16 of Davison and Hinkley (1997) again, we find the
> # conditional distribution of the ratio given the sum of city$u.
> Afn <- function(t, data) cbind(data$x-t*data$u, data$u, 1)
> ufn <- function(t, data) c(0, sum(data$u), 10)
> city.t0 <- c(mean(city$x)/mean(city$u), sqrt(city.v))
> saddle.distn(A = Afn, u = ufn, wdist = "p", type = "cond", t0 = city.t0,
+ data = city)
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.866400
Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.511350
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.622251
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.866400
Warning in dpois(y, mu, log = TRUE) : non-integer x = 8.511350
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.622251
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.210844
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.107208
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.681949
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.210844
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.107208
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.681949
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.038622
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.809279
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.152100
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.038622
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.809279
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.152100
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.452511
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.160314
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.387175
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.452511
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.160314
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.387175
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.159455
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.835832
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.004713
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.159455
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.835832
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.004713
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.155629
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.115412
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.728960
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.155629
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.115412
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.728960
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.205022
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.133273
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.661705
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.205022
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.133273
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.661705
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.722845
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.033105
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.244049
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.722845
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.033105
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.244049
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.787431
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.388294
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.824275
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.787431
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.388294
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.824275
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.415407
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.940756
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.643837
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.415407
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.940756
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.643837
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.043383
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.493218
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.463399
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.043383
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.493218
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.463399
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.671359
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.045680
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.282961
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.671359
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.045680
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.282961
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.299336
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.598142
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.102523
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.299336
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.598142
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.102523
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.433935
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.409277
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.156788
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.433935
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.409277
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.156788
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.360158
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.271008
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.368834
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.360158
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.271008
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.368834
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.319252
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.639889
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.040859
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.319252
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.639889
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.040859
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.278346
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.008770
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.712884
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.278346
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.008770
Warning in dpois(y, mu, log = TRUE) : non-integer x = 4.712884
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.237440
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.377650
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.384909
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.237440
Warning in dpois(y, mu, log = TRUE) : non-integer x = 1.377650
Warning in dpois(y, mu, log = TRUE) : non-integer x = 5.384909
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.196534
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.746531
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.056935
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.196534
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.746531
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.056935
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.155629
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.115412
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.728960
Warning in dpois(y, mu, log = TRUE) : non-integer x = 3.155629
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.115412
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.728960
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.734728
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.299661
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.965612
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.734728
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.299661
Warning in dpois(y, mu, log = TRUE) : non-integer x = 6.965612
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.228786
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.666383
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.104831
Warning in dpois(y, mu, log = TRUE) : non-integer x = 2.228786
Warning in dpois(y, mu, log = TRUE) : non-integer x = 0.666383
Warning in dpois(y, mu, log = TRUE) : non-integer x = 7.104831
Saddlepoint Distribution Approximations
Call :
saddle.distn(A = Afn, u = ufn, wdist = "p", type = "cond", t0 = city.t0,
data = city)
Quantiles of the Distribution
0.1% 1.216
0.5% 1.236
1.0% 1.248
2.5% 1.272
5.0% 1.301
10.0% 1.340
20.0% 1.393
50.0% 1.502
80.0% 1.618
90.0% 1.680
95.0% 1.732
97.5% 1.777
99.0% 1.830
99.5% 1.866
99.9% 1.938
Smoothing spline used 20 points in the range 1.182 to 2.061.
>
>
>
> cleanEx()
> nameEx("simplex")
> ### * simplex
>
> flush(stderr()); flush(stdout())
>
> ### Name: simplex
> ### Title: Simplex Method for Linear Programming Problems
> ### Aliases: simplex
> ### Keywords: optimize
>
> ### ** Examples
>
> # This example is taken from Exercise 7.5 of Gill, Murray and Wright (1991).
> enj <- c(200, 6000, 3000, -200)
> fat <- c(800, 6000, 1000, 400)
> vitx <- c(50, 3, 150, 100)
> vity <- c(10, 10, 75, 100)
> vitz <- c(150, 35, 75, 5)
> simplex(a = enj, A1 = fat, b1 = 13800, A2 = rbind(vitx, vity, vitz),
+ b2 = c(600, 300, 550), maxi = TRUE)
Linear Programming Results
Call : simplex(a = enj, A1 = fat, b1 = 13800, A2 = rbind(vitx, vity,
vitz), b2 = c(600, 300, 550), maxi = TRUE)
Maximization Problem with Objective Function Coefficients
x1 x2 x3 x4
200 6000 3000 -200
Optimal solution has the following values
x1 x2 x3 x4
0.0 0.0 13.8 0.0
The optimal value of the objective function is 41400.
>
>
>
> cleanEx()
> nameEx("smooth.f")
> ### * smooth.f
>
> flush(stderr()); flush(stdout())
>
> ### Name: smooth.f
> ### Title: Smooth Distributions on Data Points
> ### Aliases: smooth.f
> ### Keywords: smooth nonparametric
>
> ### ** Examples
>
> # Example 9.8 of Davison and Hinkley (1997) requires tilting the resampling
> # distribution of the studentized statistic to be centred at the observed
> # value of the test statistic 1.84. In the book exponential tilting was used
> # but it is also possible to use smooth.f.
> grav1 <- gravity[as.numeric(gravity[, 2]) >= 7, ]
> grav.fun <- function(dat, w, orig) {
+ strata <- tapply(dat[, 2], as.numeric(dat[, 2]))
+ d <- dat[, 1]
+ ns <- tabulate(strata)
+ w <- w/tapply(w, strata, sum)[strata]
+ mns <- as.vector(tapply(d * w, strata, sum)) # drop names
+ mn2 <- tapply(d * d * w, strata, sum)
+ s2hat <- sum((mn2 - mns^2)/ns)
+ c(mns[2] - mns[1], s2hat, (mns[2]-mns[1]-orig)/sqrt(s2hat))
+ }
> grav.z0 <- grav.fun(grav1, rep(1, 26), 0)
> grav.boot <- boot(grav1, grav.fun, R = 499, stype = "w",
+ strata = grav1[, 2], orig = grav.z0[1])
> grav.sm <- smooth.f(grav.z0[3], grav.boot, index = 3)
>
> # Now we can run another bootstrap using these weights
> grav.boot2 <- boot(grav1, grav.fun, R = 499, stype = "w",
+ strata = grav1[, 2], orig = grav.z0[1],
+ weights = grav.sm)
>
> # Estimated p-values can be found from these as follows
> mean(grav.boot$t[, 3] >= grav.z0[3])
[1] 0.0240481
> imp.prob(grav.boot2, t0 = -grav.z0[3], t = -grav.boot2$t[, 3])
$t0
[1] -1.840118
$raw
[1] 0.02254379
$rat
[1] 0.02409135
$reg
[1] 0.02224027
>
>
> # Note that for the importance sampling probability we must
> # multiply everything by -1 to ensure that we find the correct
> # probability. Raw resampling is not reliable for probabilities
> # greater than 0.5. Thus
> 1 - imp.prob(grav.boot2, index = 3, t0 = grav.z0[3])$raw
[1] 0.08678105
> # can give very strange results (negative probabilities).
>
>
>
> cleanEx()
> nameEx("tilt.boot")
> ### * tilt.boot
>
> flush(stderr()); flush(stdout())
>
> ### Name: tilt.boot
> ### Title: Non-parametric Tilted Bootstrap
> ### Aliases: tilt.boot
> ### Keywords: nonparametric
>
> ### ** Examples
>
> # Note that these examples can take a while to run.
>
> # Example 9.9 of Davison and Hinkley (1997).
> grav1 <- gravity[as.numeric(gravity[,2]) >= 7, ]
> grav.fun <- function(dat, w, orig) {
+ strata <- tapply(dat[, 2], as.numeric(dat[, 2]))
+ d <- dat[, 1]
+ ns <- tabulate(strata)
+ w <- w/tapply(w, strata, sum)[strata]
+ mns <- as.vector(tapply(d * w, strata, sum)) # drop names
+ mn2 <- tapply(d * d * w, strata, sum)
+ s2hat <- sum((mn2 - mns^2)/ns)
+ c(mns[2]-mns[1],s2hat,(mns[2]-mns[1]-orig)/sqrt(s2hat))
+ }
> grav.z0 <- grav.fun(grav1, rep(1, 26), 0)
> tilt.boot(grav1, grav.fun, R = c(249, 375, 375), stype = "w",
+ strata = grav1[,2], tilt = TRUE, index = 3, orig = grav.z0[1])
TILTED BOOTSTRAP
Exponential tilting used
First 249 replicates untilted,
Next 375 replicates tilted to -2.648,
Next 375 replicates tilted to 1.879.
Call:
tilt.boot(data = grav1, statistic = grav.fun, R = c(249, 375,
375), stype = "w", strata = grav1[, 2], tilt = TRUE, index = 3,
orig = grav.z0[1])
Bootstrap Statistics :
original bias std. error
t1* 2.846154 -0.2728113 2.513218
t2* 2.392353 -0.2450416 1.187570
t3* 0.000000 -0.7305799 2.152173
>
>
> # Example 9.10 of Davison and Hinkley (1997) requires a balanced
> # importance resampling bootstrap to be run. In this example we
> # show how this might be run.
> acme.fun <- function(data, i, bhat) {
+ d <- data[i,]
+ n <- nrow(d)
+ d.lm <- glm(d$acme~d$market)
+ beta.b <- coef(d.lm)[2]
+ d.diag <- boot::glm.diag(d.lm)
+ SSx <- (n-1)*var(d$market)
+ tmp <- (d$market-mean(d$market))*d.diag$res*d.diag$sd
+ sr <- sqrt(sum(tmp^2))/SSx
+ c(beta.b, sr, (beta.b-bhat)/sr)
+ }
> acme.b <- acme.fun(acme, 1:nrow(acme), 0)
> acme.boot1 <- tilt.boot(acme, acme.fun, R = c(499, 250, 250),
+ stype = "i", sim = "balanced", alpha = c(0.05, 0.95),
+ tilt = TRUE, index = 3, bhat = acme.b[1])
>
>
>
> cleanEx()
> nameEx("tsboot")
> ### * tsboot
>
> flush(stderr()); flush(stdout())
>
> ### Name: tsboot
> ### Title: Bootstrapping of Time Series
> ### Aliases: tsboot ts.return
> ### Keywords: nonparametric ts
>
> ### ** Examples
>
> lynx.fun <- function(tsb) {
+ ar.fit <- ar(tsb, order.max = 25)
+ c(ar.fit$order, mean(tsb), tsb)
+ }
>
> # the stationary bootstrap with mean block length 20
> lynx.1 <- tsboot(log(lynx), lynx.fun, R = 99, l = 20, sim = "geom")
>
> # the fixed block bootstrap with length 20
> lynx.2 <- tsboot(log(lynx), lynx.fun, R = 99, l = 20, sim = "fixed")
>
> # Now for model based resampling we need the original model
> # Note that for all of the bootstraps which use the residuals as their
> # data, we set orig.t to FALSE since the function applied to the residual
> # time series will be meaningless.
> lynx.ar <- ar(log(lynx))
> lynx.model <- list(order = c(lynx.ar$order, 0, 0), ar = lynx.ar$ar)
> lynx.res <- lynx.ar$resid[!is.na(lynx.ar$resid)]
> lynx.res <- lynx.res - mean(lynx.res)
>
> lynx.sim <- function(res,n.sim, ran.args) {
+ # random generation of replicate series using arima.sim
+ rg1 <- function(n, res) sample(res, n, replace = TRUE)
+ ts.orig <- ran.args$ts
+ ts.mod <- ran.args$model
+ mean(ts.orig)+ts(arima.sim(model = ts.mod, n = n.sim,
+ rand.gen = rg1, res = as.vector(res)))
+ }
>
> lynx.3 <- tsboot(lynx.res, lynx.fun, R = 99, sim = "model", n.sim = 114,
+ orig.t = FALSE, ran.gen = lynx.sim,
+ ran.args = list(ts = log(lynx), model = lynx.model))
>
> # For "post-blackening" we need to define another function
> lynx.black <- function(res, n.sim, ran.args) {
+ ts.orig <- ran.args$ts
+ ts.mod <- ran.args$model
+ mean(ts.orig) + ts(arima.sim(model = ts.mod,n = n.sim,innov = res))
+ }
>
> # Now we can run apply the two types of block resampling again but this
> # time applying post-blackening.
> lynx.1b <- tsboot(lynx.res, lynx.fun, R = 99, l = 20, sim = "fixed",
+ n.sim = 114, orig.t = FALSE, ran.gen = lynx.black,
+ ran.args = list(ts = log(lynx), model = lynx.model))
>
> lynx.2b <- tsboot(lynx.res, lynx.fun, R = 99, l = 20, sim = "geom",
+ n.sim = 114, orig.t = FALSE, ran.gen = lynx.black,
+ ran.args = list(ts = log(lynx), model = lynx.model))
>
> # To compare the observed order of the bootstrap replicates we
> # proceed as follows.
> table(lynx.1$t[, 1])
2 3 4 5 6 7 8 10 11 13 14
16 21 41 4 3 5 3 1 3 1 1
> table(lynx.1b$t[, 1])
2 3 4 5 6 7 9 10 11 12 13 15 16
2 5 24 3 2 8 2 3 38 8 2 1 1
> table(lynx.2$t[, 1])
2 3 4 5 6 7 9 11 12 13 15 17 24
18 20 41 3 4 1 1 4 1 3 1 1 1
> table(lynx.2b$t[, 1])
2 3 4 5 6 7 8 10 11 12 14 15
1 3 31 7 1 6 4 3 35 4 3 1
> table(lynx.3$t[, 1])
2 3 4 5 6 7 8 9 10 11 12 13 14 16
4 2 11 1 3 6 3 2 2 50 7 6 1 1
> # Notice that the post-blackened and model-based bootstraps preserve
> # the true order of the model (11) in many more cases than the others.
>
>
>
> cleanEx()
> nameEx("var.linear")
> ### * var.linear
>
> flush(stderr()); flush(stdout())
>
> ### Name: var.linear
> ### Title: Linear Variance Estimate
> ### Aliases: var.linear
> ### Keywords: nonparametric
>
> ### ** Examples
>
> # To estimate the variance of the ratio of means for the city data.
> ratio <- function(d,w) sum(d$x * w)/sum(d$u * w)
> var.linear(empinf(data = city, statistic = ratio))
[1] 0.03248701
>
>
>
> ### * <FOOTER>
> ###
> cleanEx()
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed: 44.104 0.665 45.086 0 0
> grDevices::dev.off()
null device
1
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')
Reference in New Issue Copy Permalink