where f [a,df1,df2] is the (a-th percentile for the F distribution with the specified degees of freedom. A function to generate this confidence interval could be written as follows.

binom.conf <- function(Y,n,alpha){

L <- Y / (Y+ (n - Y + 1) * qf(1 - alpha/2,2*(n - Y + 1),2*Y))

U <- 1 - ((n-Y) / (n - Y  + (Y + 1) * qf( 1 - alpha/2,2*(Y+1),2*(n - Y))))

return(c(L,U))}

Page 164, #2: The sign test is simply the binomial test where D = # of times X > Y, n = # non-ties, and p=0.5. If you are given the data set instead of the raw data, the command to get the # of times X > Y is: sum(X > Y). The command to find out the number of non-ties is: length(X) - sum(X==Y).

Notes on Homework Four:

1) Splus has the built in function t.test to perform the t-test (alternative is either "greater", "less", or "two.sided".)

x<-c(-30,-25,-20,-15,0,1,1,1,2,2,2,3)

t.test(x,alternative="less",mu=0)

The following function, once entered, will perform the quantile test:

quantile.test<-function(dset,pstar,xstar,alternative="two.sided"){

T1<-sum(dset<=xstar)

T2<-sum(dset< xstar)

n<-length(dset)

if (alternative=="less"){   

  pval<-1-pbinom(T2-1,n,pstar)}

else if (alternative=="greater"){

  pval<-pbinom(T1,n,pstar)}

else if(alternative=="two.sided"){

  pval<-2*min(1-pbinom(T2-1,n,pstar),pbinom(T1,n,pstar))}

return(T1,T2,pval)}

Similarly, we can make a function to construct a confidence interval for the various quantiles. One important function in doing this is sort which simply puts the data in order. A second is qbinom(alpha,n,p) which gives the "y" value corresponding to the smallest probability for that binomial that is larger than alpha. Note that this is slightly different than what you would get by taking the _closest_ one. Taking the closest may give a percent less than the 1-alpha level you set, this always gives a greater one.

quantile.interval<-function(dset,pstar,conf){   

alpha <-  1 - conf

orderstat <- sort(dset)

n <- length(dset)

r <- (qbinom((alpha/2),n,pstar)-1) + 1

s <- (qbinom((1-alpha/2),n,pstar)-1) +1

low<-orderstat[r]

hi<-orderstat[s]

realconf<-pbinom(s-1,n,pstar)-pbinom(r-1,n,pstar)

return(r,s,low,hi,realconf)}

3) To compare the powers we need to repeatedly generate a data set and see if either or both of the tests reject the null hypothesis or not. At the end we simply take the percentage of times each tests rejected and thats our estimate of the power for each. How accurate are these estimates? Well they are just percentages, and we have the formula's for finding confidence intervals for percentages. The following code would do the simulation with 100 different data sets, each data set would have 15 observations, and the power is calculated using an alpha level of 0.10. To see what the function rnorm does, simply use help(rnorm).

powert<-0

powermed<-0

for (i in 1:100){

x<-rnorm(15,0.25,1)

powert<-powert+(t.test(x,alternative="greater",mu=0)$p.value<=0.10)

powermed<-powermed+((quantile.test(x,0.5,0,alternative="greater")$pval)<=0.10)

}

For the homework problem, do the simulation 400 times, with a sample size of 20, and an alpha of 0.05.

Notes on Homework Five:

1) The function t.test will not only do the one sample t-test, put also the paired t-test. We saw on homework three how to conduct the sign test, and there is a built in function for the signed rank test, wilcox.test. You can use help(t.test) or help(wilcox.test) to get the details on these two functions. The code to run all three of these tests for the data in example 1 on page 355.

x<-c(86,71,77,68,91,72,77,91,70,71,88,87)

y<-c(88,77,76,64,96,72,65,90,65,80,81,72)



t.test(y,x,alternative="less",mu=0,paired=T)



d<-sum(y>x)

n<-length(x)-sum(x==y)

binom.test(d, n, p=0.5, alternative="less")



wilcox.test(y,x,alternative="less",mu=0,paired=T,exact=T,correct=F)

2) To construct the Hodges-Lehmann confidence interval, we simply tell S-Plus to do exactly what we would do by hand. The only complication is that the table that the S-Plus funtion that duplicates Table A12 only goes in one direction (give it the value and it gives the p). So a function is needed to "get the table value" correctly. The function below simply starts with 0 and goes up to the maximum value of the statistic, and sees which one corresponds most closely to p. You can try it with a few values and it should return the values in Table A12. Note that this function even works if p>0.5.

qsignrank<-function(p,n){

if (p<0.5){    

val<-0

for (i in 0:(n*(n+1)/4)){

if (psignrank(i,n)<=p){

val<-i

}}

p1<-psignrank(val,n)

p2<-psignrank(val+1,n)

val<-ceiling(val+max(0,(p-p1)/(p2-p1)))}

else if (p==0.5){val<-(n*(n+1)/4)}

else {val<-(n*(n+1)/2)-qsignrank((1-p),n)}

return(val)

}

HL.interval<-function(x,y,conflevel){

D <- sort(y-x)

alpha <- 1-conflevel

theavgs <- matrix(0,nrow=length(D),ncol=length(D))

wL<- qsignrank(alpha/2,length(D))    

wU<- length(D)*(length(D)+1)/2 - wL +1

for (i in 1:length(D)){

  for (j in 1:length(D)){

     theavgs[i,j] <- (D[i]+D[j])/2 }}

L<-sort(theavgs[lower.tri(theavgs)==F])[wL]

U<-sort(theavgs[lower.tri(theavgs)==F])[wU]

return(D,theavgs,L,U)}

Note that it is impossible to construct a 90%, 95%, or 99% CI for this data using the Hodges-Lehmann estimator. We can however construct an 80% one. (You do not need to answer the questions associated with the data set in the book, but part c does point this out.)

3) The set-up here will be similar to the one on Homework 4. One trick we can use, is that if we were to have two vectors x and y of length 20, then runif(20,x,y) would generate 20 uniforms, the first going from x1 to y1, the second from x2 to y2, etc... The following code would run the simulation with n=25, alpha=0.10, and 300 simulation runs. Note that this may take five minutes or so to run.

 

alphat<-0

asign<-0

asrnk<-0

for (i in 1:300){

  E<-runif(25,0,1)

  X<-runif(25,10-E,10+E)

  Y<-runif(25,10-E,10+E)

  alphat<-alphat+(t.test(Y,X,alternative="greater",

     mu=0,paired=T)$p.value<=0.10)

  d<-sum(Y>X)

  n<-length(X)-sum(X==Y)

  asign<-asign+(binom.test(d,n,p=0.5,alternative="greater")$p.value<=0.10)

  asrnk<-asrnk+(wilcox.test(Y,X,alternative="greater",

     mu=0,paired=T,exact=T,correct=F)$p.value<=0.10)

}

4) wilcox.test can also be used for the Mann-Whitney test, simply used paired=F instead of paired=T. This is also how to get the two-sample t-test using t.test.

Notes on Homework Six:

1) The function that performs the F test for two variances is var.test. The following code would analyze the data in example 1 on page 304.

pres<-c(10.8,11.1,10.4,10.1,11.3)

new<-c(10.8,10.5,11.0,10.9,10.8,10.7,10.8)



var.test(pres,new,alternative="greater")

conover.test<-function(x,y,alternative="two.sided"){

u<-abs(x-mean(x))

v<-abs(y-mean(y))

n<-length(x)

m<-length(y)

N<-n+m

R<-rank(c(u,v))

Tstat<-sum(R[1:n]^2)

Rbar2<-(1/N)*(sum(R^2))

R4<-sum(R^4)

T1<-(Tstat-n*Rbar2)/sqrt(((n*m*R4)/(N*(N-1)))-((n*m*(Rbar2^2))/(N-1)))

if (alternative=="less"){

 pval<-pnorm(T1)}

else if (alternative=="greater"){

 pval<-1-pnorm(T1)}

else{

 pval<-2*min(1-pnorm(T1),pnorm(T1))}

return(T1,Tstat,pval)}

2) The function cor.test both calculates the statistic and performs the test for both Spearman's Rho and Kendall's Tau. Using help(cor.test) explains how to use it. Note that it may give slightly different values from the book due to the way it deals with ties.

3) The command lm(y~x) performs a linear regression, for predicting y from x, and the output can be made to include the estimate of the slope and it's standard error. Simply use z<-lm(y~x) and summary(z) to see these values. The following function would make the confidence interval for the slope of the line predicting y from x, using the standard linear regression model would be:

normalslope<-function(x,y,conf.level=.95){

coeffs<-summary(lm(y~x))$coefficients

df<-length(y)-2

bhat<-coeffs[2,1]

alpha<-1-conf.level

clow<-coeffs[2,1]-qt(1-alpha/2,df)*coeffs[2,2]

chi<-coeffs[2,1]+qt(1-alpha/2,df)*coeffs[2,2]

return(bhat,conf.level,clow,chi)}

The following function performs the sequence of steps listed on page 335 for constructing the nonparametric confidence interval using the large sample approximation (equation 11 on page 317).

npslope<-function(x,y,conf.level=.95){

n<-length(x)

maxN<-n*(n-1)/2

alpha<-1-conf.level

S<-rep(0,maxN)

N<-0

for (j in 2:n){

 for (i in 1:(j-1)){

  if (x[i]!=x[j]){

    N<-N+1

    S[N]<-(y[j]-y[i])/(x[j]-x[i])}

}}

S<-sort(S[1:N])

bhat<-median(S)

wval<-qnorm(1-alpha/2)*sqrt(n*(n-1)*(2*n+5)/18)

r<-floor(0.5*(N-wval))

s<-ceiling(0.5*(N+wval)+1)

clow<-S[r]

chi<-S[s]

return(wval,r,s,S,bhat,conf.level,clow,chi)}

4 and 5) The trickiest part in using the functions for ANOVA and its nonparametric counterparts are entering the data correctly. To enter the data from example 1 on page 291 and run the ANOVA and Kruskal-Wallis test we would use the following code.

meth<-factor(c(rep("1",9),rep("2",10),rep("3",7),rep("4",8)))

obs<-c(83,91,94,89,89,96,91,92,90,91,90,81,83,84,83,88,91,89,84,

101,100,91,93,96,95,94,78,82,81,77,79,81,80,81)



summary(aov(obs~meth))



kruskal.test(obs,meth)

homeowner<-factor(c(rep(1:12,4)))

grass<-factor(c(rep(1,12),rep(2,12),rep(3,12),rep(4,12)))

obs<-c(4,4,3,3,4,2,1,2,3.5,4,4,3.5,3,2,1.5,1,2,2,3,4,1,1,2,1,

2,3,1.5,2,1,2,2,1,2,3,3,2,1,1,4,4,3,4,4,3,3.5,2,1,3.5)



friedman.test(obs,grass,homeowner)

summary(aov(obs~grass+homeowner))

Notes on Homework Seven:

Page 249 #1: The function to perform the Kolmogorov-Smirnov test is simply ks.gof. You also need to tell it which distribution you are testing. help(ks.gof) will give you the list of the ones that will work. The function for performing the Chi-squared goodness of fit test is chisq.gof. Again you need to tell it what distribution you want to use. You also need to tell it which groups you are going to put the data in, unless you want to use the automatic setting, and also how many parameters you are having estimated. help(chisq.gof) provides many of the details.

The following code first generates 50 random variables that are uniformly distributed between -2 and 2, and then does both the ks.gof and chisq.gof tests that the data is from a distribution that is uniformly distributed between -2 and 2. (So it should accept most of the time.) It uses the ten equally spaced groups for the chi-squared test, and doesn't use any estimated parameters.

It then tests if the data is normal or not... without knowing what the parameters are. (So in this case we want it to reject, but remember they do not have a lot of power.) According to the help section, it is capable of estimating the values for the normal case. In the chi-square case for this test it uses the five groups that the computer chooses itself, we need to give it the estimated values, and need to tell it how many we estimated.

x<-runif(50,min=-2,max=2)



ks.gof(x,distribution="uniform",min=-2,max=2)

chisq.gof(x,cut.points=c(-2,-1.6,-1.2,-0.8,-0.4,

0,0.4,0.8,1.2,1.6,2.0),distribution="uniform",

min=-2,max=2)



ks.gof(x,distribution="normal")

chisq.gof(x,n.classes=5,distribution="normal",

mean=mean(x),sd=sqrt(var(x)),n.param.est=2)

One extra note on this problem. If you try to enter all the data by typing it in all at once with no returns, it will lose a bunch of the numbers. This is because there is a limit to how long a line can be and still be copied. So, all you need to do is to put a few return in after a comma about half way through it if you are typing it straight into S-Plus (not recommended), or cut and paste it over in parts if you are entering it into word first.

Page 441 #1: We can use code similar to what we used to graph the EDF at the very top of this page to graph the confidence bands. The d you need to enter is the value from the table that goes with the confidence level you want. The function works in the same way as the edf function near the top of this page. It just needs to be done three times: once for the edf, once for the lower bound, and once for the upper bound.

edfconf<-function(y,d=0){

x<-sort(y)

plot(c(min(x)-1,max(x)+1),c(0,1),type="n",xlab="x",ylab="p")

lines(c(x[1]-1,x[1]),c(0,0),lty=1)

lines(c(x[1],x[1]),c(0,1/length(x)),lty=2)

for (i in 1:(length(x)-1)){

lines(c(x[i],x[i+1]),c(i/length(x),i/length(x)),lty=1)

lines(c(x[i+1],x[i+1]),c(i/length(x),(i+1)/length(x)),lty=2)}

lines(c(x[length(x)],x[length(x)]+1),c(1,1),lty=1)



lines(c(x[1]-1,x[1]),c(0,0),lty=1,lwd=2)

lines(c(x[1],x[1]),c(0,max(0,1/length(x)-d)),lty=1,lwd=2)

for (i in 1:(length(x)-1)){

lines(c(x[i],x[i+1]),c(max(0,i/length(x)-d),max(0,i/length(x)-d)),lty=1,lwd=2)

lines(c(x[i+1],x[i+1]),c(max(0,i/length(x)-d),max(0,(i+1)/length(x)-d)),

lty=1,lwd=2)}

lines(c(x[length(x)],x[length(x)]+1),c(1-d,1-d),lty=1,lwd=2)



lines(c(x[1]-1,x[1]),c(0+d,0+d),lty=1,lwd=2)

lines(c(x[1],x[1]),c(0+d,min(1,1/length(x)+d)),lty=1,lwd=2)

for (i in 1:(length(x)-1)){

lines(c(x[i],x[i+1]),c(min(1,i/length(x)+d),min(1,i/length(x)+d)),lty=1,lwd=2)

lines(c(x[i+1],x[i+1]),c(min(1,i/length(x)+d),min(1,(i+1)/length(x)+d)),

lty=1,lwd=2)}

lines(c(x[length(x)],x[length(x)]+1),c(1,1),lty=1,lwd=2)

}

Page 216, #4 The function that conducts the chi-squared tests in Section 4.2 is chisq.test. To run this test you need to enter the data in the form of a matrix. The following code is for example one on page 202. (nrow=2 means there are 2 rows, and byrow=T means that we are forming the data matrix one row at a time)

 

row1<-c(6,14,17,9)

row2<-c(30,32,17,3)

data<-matrix(c(row1,row2),nrow=2,byrow=T)

chisq.test(data,correct=F)

Page 176, #1: McNemar's test follows the same format as the chi-squared tests. The function to use is mcnemar.test. (again, using the correct=F option if you are not using the continuity correction)

Page 195, #1: The chi-squared test for differences in probabilities in section 4.1 is the same is the same as the one used in section 4.2. The only difference is that the one here is a chi-square, and the discussion on 181 uses a standard normal. Remember though that a chi-square with 1 degree of freedom is a standard normal squared. If you use chisq.test to do example 1 on page 183, you will get that the statistic is 1.6116 (= -1.2695 squared) and the p-value is 0.2043.

Notes on Homework Eight:

Page 238, #1: It may be just as easy to do this one by hand... chisq.test could be used to calculate T, and then using (3) on page 230 to calculate cramer's statistic. Unfortunately chisq.test does not work if one of the rows or columns is all zero. If we define 0/0 to be zero as described on page 329, the following code will calculate both T and Cramer's statistic.

cramer.stat<-function(x){

nr<-length(x[,1])

nc<-length(x[1,])

ex<-matrix(0,ncol=nc,nr=nrow)

Tstat<-0

for (i in 1:nr){

   for (j in 1:nc){

	ex[i,j]<-(sum(x[i,])*sum(x[,j])/sum(x))

       if((!(sum(x[i,])==0))&&(!(sum(x[,j])==0))){

       Tstat<-Tstat+(x[i,j]-ex[i,j])^2/ex[i,j]}

}}

N<-sum(x)

q<-min(length(x[,1]),length(x[1,]))

cramer<-sqrt(Tstat/(N*(q-1)))

return(Tstat,cramer)

} 

yule.stat<-function(x){

a<-x[1,1]

b<-x[1,2]

c<-x[2,1]

d<-x[2,2]

alpha<-a*d/(b*c)

phi<-(a*d-b*c)/sqrt((a+b)*(c+d)*(a+c)*(b+d))

yuleq<-(a*d-b*c)/(a*d+b*c)

if(alpha==Inf){yulecollig<-1}

else {yulecollig<-(sqrt(alpha)-1)/(sqrt(alpha)+1)}

return(alpha,phi,yuleq,yulecollig)

}

H-L interval for MWW rank sum test:

This is similar to the one for the signed rank test, and just follows the instructions on page 281-283.

HL2.interval<-function(x,y,conflevel){

alpha <- 1-conflevel

thematrix <- matrix(0,nrow=length(x),ncol=length(y))

x<-sort(x)

y<-sort(y)

big<-max(length(x),length(y))

if(big<50){

wL<-qwilcox(alpha/2,length(x),length(y))-length(x)*(length(x)+1)/2}

else {

n<-length(x)

m<-length(y)

N<-length(x)+length(y)

zp<-qnorm(alpha/2)

wp<-round(n*(N+1)/2+zp*sqrt(n*m*(N+1)/12))

wL<-wp-length(x)*(length(x)+1)/2}

wU<- length(x)*length(y)+1-wL

mat1<-matrix(rep(x,length(y)),byrow=T,nrow=length(y))

mat2<-t(matrix(rep(y,length(x)),byrow=T,nrow=length(x)))

thematrix<-mat1-mat2

L<-sort(thematrix)[wL]

U<-sort(thematrix)[wU]

return(wL,wU,thematrix,L,U)}