R version 4.2.1 (2022-06-23) -- "Funny-Looking Kid"
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("g10.pdf");options(width=64)
> #setwd("C:\\Users\\kolassa\\Class555")
> setwd("~/Taught1/960-555/Data")
> #**************************************************************/
> #* 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)))
> # For Hotelling and multivariate rank tests resp:
> library(Hotelling)#for HotellingsT2
> library(ICSNP)#for rank.ctest
> cat('\n One-sample Hotelling Test\n')

 One-sample Hotelling Test
> HotellingsT2(bp[,c("spd","dpd")])

	Hotelling's one sample T2-test

data:  bp[, c("spd", "dpd")]
T.2 = 31.803, df1 = 2, df2 = 13, p-value = 9.839e-06
alternative hypothesis: true location is not equal to c(0,0)

> cat('\n Multivariate Sign Test\n')

 Multivariate Sign Test
> rank.ctest(bp[,c("spd","dpd")],scores="sign")

	Marginal One Sample Sign Test

data:  bp[, c("spd", "dpd")]
T = 15, df = 2, p-value = 0.0005531
alternative hypothesis: true location is not equal to c(0,0)

> cat('\n Multivariate Signed Rank Test\n')

 Multivariate Signed Rank Test
> rank.ctest(bp[,c("spd","dpd")])

	Marginal One Sample Signed Rank Test

data:  bp[, c("spd", "dpd")]
T = 11.636, df = 2, p-value = 0.002973
alternative hypothesis: true location is not equal to c(0,0)

> cat('\n Size of Multivariate Tests\n')

 Size of Multivariate Tests
> library(NonparametricHeuristic)#For fun.comparepower
> fun.comparepower(c(20,40),ndim=2,nsamp=20000,
+   dist=c("rnorm","rcauchy","rlaplace"))
, , two.sided, 0, 20

                rnorm rcauchy rlaplace
T2             0.0500 0.01740  0.04615
Sign test      0.0433 0.04030  0.04450
Sign rank test 0.0405 0.03965  0.04165

, , two.sided, 0, 40

                 rnorm rcauchy rlaplace
T2             0.05135 0.01565  0.04715
Sign test      0.04720 0.04865  0.04765
Sign rank test 0.04720 0.04900  0.04755

> cat('\n Power of Multivariate Tests\n')

 Power of Multivariate Tests
> fun.comparepower(c(20,40),ndim=2,nsamp=20000,
+   dist=c("rnorm","rcauchy","rlaplace"),hypoth=.5)
, , two.sided, 0, 20

                rnorm rcauchy rlaplace
T2             0.7461 0.07935  0.48015
Sign test      0.7026 0.26015  0.54975
Sign rank test 0.5264 0.31900  0.54610

, , two.sided, 0, 40

                 rnorm rcauchy rlaplace
T2             0.97775 0.08810  0.80095
Sign test      0.96985 0.51165  0.89585
Sign rank test 0.87210 0.63990  0.89705

> library(MultNonParam)#For shiftcr
> shiftcr(bp[,c("dpd","spd")])
> 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)
> library(Hotelling)#for hotelling.test
> print(hotelling.test(wheat$Huntsman+wheat$Atou~wheat$region))
Test stat:  13.24 
Numerator df:  2 
Denominator df:  7 
P-value:  0.03279 
> t.test(wheat$Huntsman~wheat$region);t.test(wheat$Atou~wheat$region)

	Welch Two Sample t-test

data:  wheat$Huntsman by wheat$region
t = 3.4925, df = 7.7873, p-value = 0.008522
alternative hypothesis: true difference in means between group Other and group North is not equal to 0
95 percent confidence interval:
 0.4698619 2.3221381
sample estimates:
mean in group Other mean in group North 
              6.432               5.036 


	Welch Two Sample t-test

data:  wheat$Atou by wheat$region
t = 3.4809, df = 6.7947, p-value = 0.01075
alternative hypothesis: true difference in means between group Other and group North is not equal to 0
95 percent confidence interval:
 0.4652694 2.4747306
sample estimates:
mean in group Other mean in group North 
              6.762               5.292 

> cat('Wheat rank test, normal theory p-values')
Wheat rank test, normal theory p-values> print(hotelling.test(rank(wheat$Huntsman)+rank(wheat$Atou)~wheat$region))
Test stat:  14.797 
Numerator df:  2 
Denominator df:  7 
P-value:  0.0256 
> #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')
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)
Hotelling    t test 
    0.023     0.008 
> print(hotelling.test(wheat$Huntsman+wheat$Atou~wheat$region,perm=T,
+   progBar=FALSE))
Test stat:  5.7926 
Numerator df:  2 
Denominator df:  7 
Permutation P-value:  0.0159 
Number of permutations : 10000 
> library(ICSNP)#For rank.ctest and HotelingsT2.
> rank.ctest(cbind(wheat$Huntsman,wheat$Atou)~wheat$region)

	Marginal Two Sample Rank Sum Test

data:  cbind(wheat$Huntsman, wheat$Atou) by wheat$region
T = 6.4908, df = 2, p-value = 0.03895
alternative hypothesis: true location difference is not equal to c(0,0)

> rank.ctest(cbind(wheat$Huntsman,wheat$Atou)~wheat$region,scores="sign")

	Marginal Two Sample Median Test

data:  cbind(wheat$Huntsman, wheat$Atou) by wheat$region
T = 7, df = 2, p-value = 0.0302
alternative hypothesis: true location difference is not equal to c(0,0)

> #*************************************************************/
> # 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)
> cat('\n Density estimation  \n')

 Density estimation  
> #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"))
> 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))
>