R version 4.2.2 Patched (2022-11-10 r83330) -- "Innocent and Trusting"
Copyright (C) 2022 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.

> pdf("g13.pdf");options(width=64)
> #setwd("C:\\Users\\kolassa\\Class555")
> setwd("~/Taught1/960-555/Data")
> #*************************************************************/
> # Data from http://lib.stat.cmu.edu/datasets/Arsenic         */
> # reformatted into ASCII on the course home page.  Data re-  */
> # flect arsenic levels in toenail clippings; covariates in-  */
> # clude age, sex (1=M), categorical measures of quantities   */
> # used for drinking and cooking, arsenic in the water, and   */
> # arsenic in the nails.  To make arsenic.dat from Arsenic, do*/
> #antiword Arsenic|awk '((NR>39)&&(NR<61)){print}'>arsenic.dat*/
> #*************************************************************/
> arsenic<-read.table("arsenic.dat")
> names(arsenic)<-c("age","sex","drink","cook","water",
+    "nails")
> arsenic$status<-(arsenic$water>=.0001)+0
> arsenic$water[arsenic$water<.0001]<-.0001
> arsenic$lnail<-log(arsenic$nails)
> arsenic$lwater<-log(arsenic$water)
> library(boot)#For boot and boot.ci
> #Define the function to be applied to data sets to get
> #the parameter to be bootstrapped.  
> boot.ci(boot(arsenic$nails,function(x,index) 
+    return(median(x[index])),9999))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 9999 bootstrap replicates

CALL : 
boot.ci(boot.out = boot(arsenic$nails, function(x, index) return(median(x[index])), 
    9999))

Intervals : 
Level      Normal              Basic         
95%   ( 0.0141,  0.2750 )   ( 0.0400,  0.2310 )  

Level     Percentile            BCa          
95%   ( 0.119,  0.310 )   ( 0.118,  0.277 )  
Calculations and Intervals on Original Scale
> library(MultNonParam)#For exactquantileci
> exactquantileci(arsenic$nails)
$cis
      [,1]
[1,] 0.118
[2,] 0.358

$pvals
 [1]           NA           NA           NA           NA
 [5]           NA           NA           NA           NA
 [9]           NA           NA           NA           NA
[13]           NA           NA           NA 9.536743e-07

> logscale<-function(x,index) return(log(sd(x[index])))
> sdbootsamp<-boot(arsenic$nails,logscale,9999)
> sdoutput<-boot.ci(sdbootsamp)
> plot(density(sdbootsamp$t))
> boot.ci(boot(arsenic$nails,logscale,99999))$bca
     conf                                        
[1,] 0.95 11112.82 99951.97 -1.647932 -0.09432326
> library(NonparametricHeuristic)#for fun.testboot
> t1<-fun.testboot(function(data,indices) return(median(data[indices])),
+    mcsamp=1000, sampsize=20, dists=c("rexp","runif","rcauchy"), 
+    true=c(log(2),.5,0),alpha=.1)
> t1
$cilen
           Normal     Basic Percentile       BCa
rexp    0.7774421 0.7436943  0.7436943 0.7068057
runif   0.3453235 0.3379974  0.3379974 0.3335132
rcauchy 1.4303642 1.3321451  1.3321451 1.3442736

$ok
        Normal Basic Percentile   BCa
rexp     0.839 0.750      0.889 0.891
runif    0.829 0.734      0.886 0.882
rcauchy  0.908 0.824      0.870 0.873

> twinbrain<-as.data.frame(scan("IQ_Brain_Size",
+   what=list(CCMIDSA=0,FIQ=0,HC=0,ORDER=0,PAIR=0,SEX=0,
+   TOTSA=0, TOTVOL=0,WEIGHT=0),skip=27,nmax=20))
> fir<-twinbrain[twinbrain$ORDER==1,]
> fir$v1<-fir$TOTVOL
> sec<-twinbrain[twinbrain$ORDER==2,]
> sec$v2<-sec$TOTVOL
> brainpairs<-merge(fir,sec,by="PAIR")[,c("v1","v2")]
> brainpairs$diff<-brainpairs$v2-brainpairs$v1
> cat("\nParametric Bootstrap for median difference\n")

Parametric Bootstrap for median difference
> qmed<-function(x,indices) return(median(x[indices]))
> ran.diff<-function(x,ests)
+   return(ests[1]+ests[2]*rnorm(length(x)))
> bootout<-boot(brainpairs$diff,qmed,R=999,sim="parametric",
+   ran.gen=ran.diff,mle=c(mean(brainpairs$diff),
+   sd(brainpairs$diff)))
> boot.ci(bootout,type=c("basic","norm"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = bootout, type = c("basic", "norm"))

Intervals : 
Level      Normal              Basic         
95%   (-29.354,  52.330 )   (-31.585,  52.688 )  
Calculations and Intervals on Original Scale
> rho<-function(d,idx) return(cor(d[idx,1],d[idx,2]))
> bootout<-boot(brainpairs[,c("v1","v2")],rho,9999)
> hist(bootout$t,freq=FALSE)
> legend(-.4,6,lty=1:2,
+    legend=c("Estimate","Confidence Bounds"))
> abline(v=bootout$t0)
> abline(v=sort(bootout$t)[(bootout$R+1)*c(.025,.975)],lty=2)
> boot.ci(bootout)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 9999 bootstrap replicates

CALL : 
boot.ci(boot.out = bootout)

Intervals : 
Level      Normal              Basic         
95%   ( 0.7457,  1.1082 )   ( 0.8515,  1.1506 )  

Level     Percentile            BCa          
95%   ( 0.6783,  0.9773 )   ( 0.5513,  0.9722 )  
Calculations and Intervals on Original Scale
> #Second component is an estimate of its scale, squared.
> regparam<-function(d,idx) return(summary(
+    lm(d[idx,2]~d[idx,1]))$coefficients[2,1:2]^(1:2))
> boot.ci(boot(brainpairs[,c("v1","v2")],regparam,9999))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 9999 bootstrap replicates

CALL : 
boot.ci(boot.out = boot(brainpairs[, c("v1", "v2")], regparam, 
    9999))

Intervals : 
Level      Normal              Basic             Studentized     
95%   ( 0.2837,  1.1408 )   ( 0.1705,  0.9215 )   ( 0.4806,  1.0082 )  

Level     Percentile            BCa          
95%   ( 0.6320,  1.3830 )   ( 0.6001,  1.2114 )  
Calculations and Intervals on Original Scale
> medianiqr2<-function(x,index) return(c(median(x[index]),
+ (quantile(x[index],.75)-quantile(x[index],.25))^2))
> medianvar<-function(x,index) return(c(median(x[index]),
+ var(x[index])))
> boot.ci(boot(arsenic$nails,medianiqr2,99999),type="stud")
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 99999 bootstrap replicates

CALL : 
boot.ci(boot.out = boot(arsenic$nails, medianiqr2, 99999), type = "stud")

Intervals : 
Level    Studentized     
95%   ( 0.0330,  0.2606 )  
Calculations and Intervals on Original Scale
> boot.ci(boot(arsenic$nails,medianvar,99999),type="stud")
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 99999 bootstrap replicates

CALL : 
boot.ci(boot.out = boot(arsenic$nails, medianvar, 99999), type = "stud")

Intervals : 
Level    Studentized     
95%   (-0.0753,  0.3075 )  
Calculations and Intervals on Original Scale
> residuals<-resid(lm(brainpairs$v2~brainpairs$v1))
> fittedv<-brainpairs$v2-residuals
> regparam<-function(data,indices,fittedv,residuals){
+  y<-fittedv+residuals[indices]
+  return(summary(lm(y~data[,1]))$coefficients[2,1:2]^(1:2))
+ }
> library(boot)#For boot and boot.ci
> bootout<-boot(brainpairs[,c("v1","v2")],regparam,9999,
+   fittedv=fittedv,residuals=residuals)
> boot.ci(bootout,type="stud")
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 9999 bootstrap replicates

CALL : 
boot.ci(boot.out = bootout, type = "stud")

Intervals : 
Level    Studentized     
95%   ( 0.5062,  1.0498 )  
Calculations and Intervals on Original Scale
> cat("\n Fixed X bootstrap adjusting for variance  \n")

 Fixed X bootstrap adjusting for variance  
> modelfit<-lm(brainpairs$v2~brainpairs$v1)
> hat<-lm.influence(modelfit)$hat
> adjresid<-resid(modelfit)/sqrt(1-hat)
> newregparam<-function(data,indices,fittedv,adjresid,hat){
+  y<-fittedv+sqrt(1-hat)*adjresid[indices]
+  return(summary(lm(y~data[,1]))$coefficients[2,1:2]^(1:2))
+ }
> bootout<-boot(brainpairs[,c("v1","v2")],newregparam,9999,
+   fittedv=fittedv,adjresid=adjresid,hat=hat)
> boot.ci(bootout)$student[4:5]
[1] 0.5682148 0.9899201
> #**********************************************************/
> # Chicken weightdata from Example K of Cox & Snell (1981).*/
> # Variables represent protein source, protein level, fish */
> # solubles level.  Two responses representing weight at   */
> # 16wk in g follow.                                       */
> #**********************************************************/
> temp1<-temp<-as.data.frame(scan("chicken.dat",
+     what= list(source="", lev=0,fish=0,weight=0,w2=0)))
> temp1$weight<-temp1$w2;temp$h<-0;temp1$h<-1
> temp$w2<-NULL;temp1$w2<-NULL
> chicken<-rbind(temp,temp1)
> meandiff<-function(data,indices){
+   fitv<-fitted(lm(data[,1]~as.factor(data[,2])))
+   residv<-data[,1]-fitv
+   y<-fitv+residv[indices]
+   return(diff(range(unlist(lapply(split(y,data[,2]),
+      "mean")))))
+ }
> cat("Bootstrap CI for maximal difference in means\n")
Bootstrap CI for maximal difference in means
> boot.ci(boot(cbind(chicken$weight,chicken$lev),meandiff,9999))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 9999 bootstrap replicates

CALL : 
boot.ci(boot.out = boot(cbind(chicken$weight, chicken$lev), meandiff, 
    9999))

Intervals : 
Level      Normal              Basic         
95%   (-14.1, 709.3 )   (-32.6, 681.5 )  

Level     Percentile            BCa          
95%   (104.5, 818.6 )   ( 73.3, 752.2 )  
Calculations and Intervals on Original Scale
> library(bootstrap)#For jackknife
> jackknife(arsenic$nails,median)
$jack.se
[1] 0.1224907

$jack.bias
[1] 0.4033333

$jack.values
 [1] 0.2220 0.2220 0.2220 0.2220 0.1665 0.1665 0.2220 0.2220
 [9] 0.1665 0.2220 0.2220 0.1665 0.1665 0.1665 0.1665 0.1665
[17] 0.1665 0.2220 0.1665 0.2220 0.2135

$call
jackknife(x = arsenic$nails, theta = median)

> jackknife(arsenic$nails,mean)
$jack.se
[1] 0.1062373

$jack.bias
[1] 0

$jack.values
 [1] 0.37880 0.37885 0.37980 0.37885 0.37090 0.36685 0.38075
 [8] 0.37685 0.36925 0.37950 0.38110 0.34315 0.35890 0.27215
[15] 0.34220 0.37130 0.36310 0.37770 0.37100 0.37800 0.37600

$call
jackknife(x = arsenic$nails, theta = mean)

> jackknife(arsenic$nails,mean,trim=0.25)
$jack.se
[1] 0.04604028

$jack.bias
[1] -0.02450216

$jack.values
 [1] 0.2216 0.2217 0.2217 0.2217 0.2058 0.1977 0.2217 0.2177
 [9] 0.2025 0.2217 0.2217 0.1977 0.1977 0.1977 0.1977 0.2066
[17] 0.1977 0.2194 0.2060 0.2200 0.2160

$call
jackknife(x = arsenic$nails, theta = mean, trim = 0.25)

> jackknife(brainpairs$diff,median)
$jack.se
[1] 54

$jack.bias
[1] 0

$jack.values
 [1] 28 28 28 -8 -8 -8 28 28 -8 -8

$call
jackknife(x = brainpairs$diff, theta = median)

> t2<-testjack(list(mean,median,mean),dists=list(rexp,rnorm),
+   others=list(NULL,NULL,list(trim=0.25)),nsamp=100000)
> t2
, , Expectation of Statistic

               mean       median          mean
rexp   1.0004483028 0.7385140765  0.8100448677
rnorm -0.0008802304 0.0003883707 -0.0004074243

, , Expectation of Bias Estimate

              mean       median          mean
rexp  1.137979e-19  0.088473850 -6.027974e-02
rnorm 2.125557e-19 -0.003364226  6.698722e-05

, , Parameter

            mean       median        mean
rexp  0.99040938 0.6930099309 0.746744242
rnorm 0.01149362 0.0007123214 0.002838278

>