pdf("g10.pdf");options(width=64)
#setwd("C:\\Users\\kolassa\\Class555")
setwd("~/Taught1/960-555/Data")
<pre>
#**************************************************************/
#* Data on blood presure change after captopril from          */
#* http://www.stat.ucla.edu/~rgould/datasets/bloodpressure.dat*/
#* Variables are                                              */
#* before systolic bp, after systolic bp, systolic bp change  */
#* before diastolic bp, after diastolic bp, diast. bp change  */
#* Example E from Cox and Snell, Applied Statistics: Princi-  */
#* ples and examples.  First line has column headings.  Fields*/
#* are separated by tabs.  SAS will give an error with this.  */
#**************************************************************/
bp<-as.data.frame(scan('bloodpressure.dat',skip=1,
  what=list(spb=0,spa=0,spd=0,dpb=0,dpa=0,dpd=0)))
#<a name=block1> Block 1</a>
# For Hotelling and multivariate rank tests resp:
library(Hotelling)#for HotellingsT2
library(ICSNP)#for rank.ctest
cat('\n One-sample Hotelling Test\n')
HotellingsT2(bp[,c("spd","dpd")])
cat('\n Multivariate Sign Test\n')
rank.ctest(bp[,c("spd","dpd")],scores="sign")
cat('\n Multivariate Signed Rank Test\n')
rank.ctest(bp[,c("spd","dpd")])
#<a name=block2> Block 2</a>
cat('\n Size of Multivariate Tests\n')
library(NonparametricHeuristic)#For fun.comparepower
fun.comparepower(c(20,40),ndim=2,nsamp=20000,
  dist=c("rnorm","rcauchy","rlaplace"))
cat('\n Power of Multivariate Tests\n')
fun.comparepower(c(20,40),ndim=2,nsamp=20000,
  dist=c("rnorm","rcauchy","rlaplace"),hypoth=.5)
#<a name=block3> Block 3</a>
library(MultNonParam)#For shiftcr
shiftcr(bp[,c("dpd","spd")])
#<a name=block4> Block 4</a>
wheat<-as.data.frame(scan("set5.data",what=list(variety="",
  y0=0,y1=0,y2=0,y3=0,y4=0,y5=0,y6=0,y7=0,y8=0,y9=0),
  na.strings="-"))
# Observations are represented by columns rather than by
# rows.  Swap this.  New column names are in first column.
dimnames(wheat)[[1]]<-wheat[,1]
wheat<-as.data.frame(t(wheat[,-1]))
dimnames(wheat)[[1]]<-c("E1","E2","N3","N4","N5","N6","W7",
 "E8","E9","N10")
wheat$region<-factor(c("North","Other")[1+(
   substring(dimnames(wheat)[[1]],1,1)!="N")],
   c("Other","North"))
plot(wheat$Huntsman,wheat$Atou,pch=(wheat$region=="North")+1,
  main="Wheat Yields")
legend(6,5,legend=c("Other","North"),pch=1:2)
#<a name=block5> Block 5</a>
library(Hotelling)#for hotelling.test
print(hotelling.test(wheat$Huntsman+wheat$Atou~wheat$region))
#<a name=block6> Block 6</a>
t.test(wheat$Huntsman~wheat$region);t.test(wheat$Atou~wheat$region)
#<a name=block7> Block 7</a>
cat('Wheat rank test, normal theory p-values')
print(hotelling.test(rank(wheat$Huntsman)+rank(wheat$Atou)~wheat$region))
#<a name=block8> Block 8</a>
#Brute-force way to get estimate of permutation p-value for
#both T2 and the max t statistic.
cat('Permutation Tests for Wheat Data, Brute Force')
obsh<-hotelling.test(wheat$Huntsman+wheat$Atou~wheat$region)$stats$statistic
obst<-max(c(t.test(wheat$Huntsman~wheat$region)$statistic,
   t.test(wheat$Atou~wheat$region)$statistic))
out<-array(NA,c(1000,2))
dimnames(out)<-list(NULL,c("Hotelling","t test"))
for(j in seq(dim(out)[[1]])){
   newr<-sample(wheat$region,length(wheat$region))
   hto<-hotelling.test(wheat$Huntsman+wheat$Atou~newr)
   out[j,1]<-hto$stats$statistic>=obsh
   out[j,2]<-max(t.test(wheat$Huntsman~newr)$statistic,
     t.test(wheat$Atou~newr)$statistic)>=obst
}
apply(out,2,mean)
#<a name=block9> Block 9</a>
print(hotelling.test(wheat$Huntsman+wheat$Atou~wheat$region,perm=T,
  progBar=FALSE))
#<a name=block10> Block 10</a>
library(ICSNP)#For rank.ctest and HotelingsT2.
rank.ctest(cbind(wheat$Huntsman,wheat$Atou)~wheat$region)
#<a name=block11> Block 11</a>
rank.ctest(cbind(wheat$Huntsman,wheat$Atou)~wheat$region,scores="sign")
#*************************************************************/
# 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=block12> Block 12</a>
cat('\n Density estimation  \n')
#Save the density object at the same time it is plotted.
plot(a<-density(arsenic$nails),lty=1,
  main="Density for Arsenic in Nails for Various Kernels")
lines(density(arsenic$nails,kernel="epanechnikov"),lty=2)
lines(density(arsenic$nails,kernel="triangular"),lty=3)
lines(density(arsenic$nails,kernel="rectangular"),lty=4)
legend(1,2, lty=rep(1,4), legend=c("Normal","Quadratic",
   "Triangular","Rectangular"))
#<a name=block13> Block 13</a>
plot(density(arsenic$nails,bw=a$bw/10),xlab="Arsenic",type="l",
    main="Density estimate with inopportune bandwidths")
lines(density(arsenic$nails,bw=a$bw*5),lty=2)
legend(1,2,legend=paste("Band width default",c("/10","*5")),
   lty=c(1,2))
</pre>