pdf("g13.pdf");options(width=64)
#setwd("C:\\Users\\kolassa\\Class555")
setwd("~/Taught1/960-555/Data")
<pre>
#*************************************************************/
# 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)
#<a name=block1> Block 1</a>
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))
#<a name=block2> Block 2</a>
library(MultNonParam)#For exactquantileci
exactquantileci(arsenic$nails)
#<a name=block3> Block 3</a>
logscale<-function(x,index) return(log(sd(x[index])))
sdbootsamp<-boot(arsenic$nails,logscale,9999)
sdoutput<-boot.ci(sdbootsamp)
#<a name=block4> Block 4</a>
plot(density(sdbootsamp$t))
#<a name=block5> Block 5</a>
boot.ci(boot(arsenic$nails,logscale,99999))$bca
#<a name=block6> Block 6</a>
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
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
#<a name=block7> Block 7</a>
cat("\nParametric Bootstrap for median difference\n")
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"))
#<a name=block8> Block 8</a>
rho<-function(d,idx) return(cor(d[idx,1],d[idx,2]))
#<a name=block9> Block 9</a>
bootout<-boot(brainpairs[,c("v1","v2")],rho,9999)
#<a name=block10> Block 10</a>
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)
#<a name=block11> Block 11</a>
boot.ci(bootout)
#<a name=block12> Block 12</a>
#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))
#<a name=block13> Block 13</a>
boot.ci(boot(brainpairs[,c("v1","v2")],regparam,9999))
#<a name=block14> Block 14</a>
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")
boot.ci(boot(arsenic$nails,medianvar,99999),type="stud")
#<a name=block15> Block 15</a>
residuals<-resid(lm(brainpairs$v2~brainpairs$v1))
fittedv<-brainpairs$v2-residuals
#<a name=block16> Block 16</a>
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)
#<a name=block17> Block 17</a>
boot.ci(bootout,type="stud")
#<a name=block18> Block 18</a>
cat("\n Fixed X bootstrap adjusting for variance  \n")
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]
#**********************************************************/
# 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)
#<a name=block19> Block 19</a>
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")
boot.ci(boot(cbind(chicken$weight,chicken$lev),meandiff,9999))
#<a name=block20> Block 20</a>
library(bootstrap)#For jackknife
jackknife(arsenic$nails,median)
#<a name=block21> Block 21</a>
jackknife(arsenic$nails,mean)
#<a name=block22> Block 22</a>
jackknife(arsenic$nails,mean,trim=0.25)
#<a name=block23> Block 23</a>
jackknife(brainpairs$diff,median)
#<a name=block24> Block 24</a>
t2<-testjack(list(mean,median,mean),dists=list(rexp,rnorm),
  others=list(NULL,NULL,list(trim=0.25)),nsamp=100000)
t2
</pre>