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("g05.pdf");options(width=64)
> #setwd("C:\\Users\\kolassa\\Class555")
> setwd("~/Taught1/960-555/Data")
> #********************************************************/
> #* Yarn strength data from Example Q of Cox & Snell     */
> #* (1981).  Variables represent strength of two types of*/
> #* yarn collected from six different bobbins.           */
> #********************************************************/
> yarn<-as.data.frame(scan("yarn.dat",what=
+     list(strength=0, bobbin=0,type="")))
> yarn$ab<-pmin(rank(yarn$strength),rank(-yarn$strength))
> library(NonparametricHeuristic)#For siegel.tukey.ranks
> yarn$st<-round(siegel.tukey.ranks(yarn$strength),2)
> yarnranks<-yarn[order(yarn$strength),
+    c("strength","type","ab","st")]
> library(DescTools)#For SiegelTukeyTest
> SiegelTukeyTest(strength~type,data=yarn)

	Siegel-Tukey-test for equal variability

data:  strength by type
ST = 608.67, p-value = 0.7179
alternative hypothesis: true ratio of scales is not equal to 1

> yarnsplit<-split(yarn$strength,yarn$type)
> ansari.test(yarnsplit[[1]],yarnsplit[[2]])

	Ansari-Bradley test

data:  yarnsplit[[1]] and yarnsplit[[2]]
AB = 310, p-value = 0.6786
alternative hypothesis: true ratio of scales is not equal to 1

> #*************************************************************/
> # 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(NonparametricHeuristic)#For invertsigntest
> invertsigntest(split(arsenic$nails,arsenic$sex))
[1] -0.158  0.278
> wilcox.test(nails~sex,data=arsenic,conf.int=TRUE)

	Wilcoxon rank sum test with continuity correction

data:  nails by sex
W = 50, p-value = 0.9135
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -0.2780188  0.1580615
sample estimates:
difference in location 
          -0.008916766 

> par(mfrow=c(1,1))
> plot(range(yarn$strength),c(0,1),type="n",
+   main="Yarn Strength",xlab="Yarn Strength",
+   ylab="Probability")
> yarnsplit<-split(yarn$strength,yarn$type)
> lines(ecdf(yarnsplit[[1]]),col=1)
> lines(ecdf(yarnsplit[[2]]),col=2)
> legend(17,.2,lty=c(1,1),col=c(1,2),legend=c("A","B"))
> ks.test(yarnsplit[[1]],yarnsplit[[2]])

	Exact two-sample Kolmogorov-Smirnov test

data:  yarnsplit[[1]] and yarnsplit[[2]]
D = 0.29167, p-value = 0.2166
alternative hypothesis: two-sided

> #Library CvM2SL2Test must be installed from the archives.
> #library(devtools)
> #install_version("CvM2SL2Test",version="2.0-1")
> library(CvM2SL2Test)#For cvmts.test and cmvts.pval
> cvmstat<-cvmts.test(yarnsplit[[1]],yarnsplit[[2]])
> cvmts.pval(cvmstat,length(yarnsplit[[1]]),
+   length(yarnsplit[[2]]))
[1] 0.09474231
> maize<-as.data.frame(scan("T58.1",what=list(exno=0,tabno=0,
+   lineno=0,loc="",block="",plot=0,trt="",ears=0, wght=0)))
> maize$wght[maize$wght==-9999]<-NA
> #Location TEAN has no tied values.  R treats ranks of
> #missing values nonintuitively.  Remove missing values.
> tean<-maize[(maize$loc=="TEAN")&(!is.na(maize$wght)),]
> cat('\n Kruskal Wallis H Test for Maize Data  \n')

 Kruskal Wallis H Test for Maize Data  
> kruskal.test(split(tean$wght,tean$trt))

	Kruskal-Wallis rank sum test

data:  split(tean$wght, tean$trt)
Kruskal-Wallis chi-squared = 9.256, df = 11, p-value =
0.5983

> #Alternative R syntax:
> #kruskal.test(tean$wght,tean$trt)
> #Note that treatment is already a factor.
> anova(lm(wght~trt,data=tean))
Analysis of Variance Table

Response: wght
          Df  Sum Sq Mean Sq F value Pr(>F)
trt       11  9.0072 0.81883  0.7248 0.7045
Residuals 23 25.9831 1.12970               
> anova(lm(rank(wght,na.last=NA)~trt,data=tean))
Analysis of Variance Table

Response: rank(wght, na.last = NA)
          Df  Sum Sq Mean Sq F value Pr(>F)
trt       11  971.87  88.352  0.7821 0.6547
Residuals 23 2598.13 112.962               
> library(NonparametricHeuristic)#For showmultigroupscoretest
> showmultigroupscoretest(c(3,3,4))
> showmultigroupscoretest(c(3,3,4),fun.givescore(1:10,sv="ns"))
> library(exactRankTests)#Gives savage, normal scores
> library(NonparametricHeuristic)#For genmultscore
> cat("Other scoring schemes: Normal Scores\n")
Other scoring schemes: Normal Scores
> genmultscore(tean$wght,tean$trt,
+    cscores(tean$wght,type="Normal"))

	General Rank Score Multi-Group Test

data:  NA
= 3.5747, df = 11, p-value = 0.9808

> cat("Other scoring schemes: Savage Scores\n")
Other scoring schemes: Savage Scores
> genmultscore(tean$wght,tean$trt,
+    cscores(tean$wght,type="Savage"))

	General Rank Score Multi-Group Test

data:  NA
= 3.2046, df = 11, p-value = 0.9877

> #date()
> #aov.P(tean$wght[!is.na(tean$wght)],
> #    tean$trt[!is.na(tean$wght)])
> #date()
> obsp<-anova(lm(wght~trt,data=tean))[[4]][1]
> out<-rep(NA,10000)
> #Monte Carlo approximation to the permutation distribution
> for(i in seq(length(out))){
+   out[i]<-anova(lm(sample(wght)~trt,data=tean))[[4]][1]
+ }
> mean(out>=obsp)
[1] 0.7151
>