# Chapter 9 Script
# November 13, 2007
# Alan T. Arnholt

# Example 9.3
ALPHA <- round(1 - pnorm(2, 1, 2), 2)
ALPHA

BETA <- round(pnorm(2,4,2), 2)
BETA

# Partial code for Figure 9.1
curve(dnorm(x,1,2),-5,10,ylab="")
lines(c(1,1),c(0,dnorm(1,1,2)))
lines(c(2,2),c(0,dnorm(2,1,2)))
lines(c(4,4),c(0,dnorm(4,4,2)))
lines(c(-5,10),c(0,0),lwd=3)
curve(dnorm(x,4,2),add=TRUE)


#  Consider the graph below to help with the Example 9.4
#  Graphical Representation of Power for H0: theta=2 versus H1: theta=10
curve(dexp(x,rate=2),0,3,xlim=c(0,2),ylim=c(0,10),col="blue",ylab="")
cv <- qexp(.05,2)
lines(c(cv,cv),c(0,dexp(cv,2)),lwd=2)
lines(c(cv,cv),c(0,dexp(cv,10)),col="red")
curve(dexp(x,rate=10),0,2,add=TRUE,col="red")
lines(c(0,2),c(0,0),lwd=3)
xs <- seq(cv,0,-.0001)
polygon(x=c(0,cv,cv,xs,0),y=c(0,0,dexp(cv,10),
dexp(xs,10),0),density=-1,col="red")

# Example 9.5
# part a.
PtypeIerror <- pnorm(36,40,6/sqrt(9)) + 1 - pnorm(44,40,6/sqrt(9))
PtypeIerror
# part b.
PtypeIerror <- pnorm(38,40,6/sqrt(36)) + 1 - pnorm(42,40,6/sqrt(36))
PtypeIerror
# part c. Power(Mu) - Figure 9.2
mu <- seq(30,50,.01)
power9 <- 1-pnorm(44,mu,6/sqrt(9)) + pnorm(36,mu,6/sqrt(9))
power36 <- 1-pnorm(42,mu,6/sqrt(36)) + pnorm(38,mu,6/sqrt(36))
plot(mu, power9, type="l", ylab=expression(Power(mu)), xlab=expression(mu),
ylim=c(0,1))
lines(mu,power36,type="l")
arrows(32, 0.6 , 34.2, .78, lwd=2, length=0.05)
arrows(32, 0.35 , 37, .78, lwd=2, length=0.05)
arrows(40, 0.4 , 40, 0.06, lwd=2, length=0.05)
text(32,0.58,expression(n==9))
text(32.3,0.33,expression(n==36))
text(40,0.45,expression(alpha==0.045))
lines(c(30,50),c(0.0455,0.0455),lty=2)
lines(c(30,50),c(0.0,0.0))

# Example 9.6
# part a.
ALPHA <- 1 - pnorm(2.036,1,1)
round(ALPHA, 3)
BETA <-  pnorm(2.036,2,1)
round(BETA, 3)
POWER <- 1- BETA
round(POWER, 3)
# part b.
ALPHA <- pnorm(1.300,1,1) - pnorm(1.100,1,1) + (1-pnorm(2.461,1,1))
round(ALPHA, 3)
BETA <- pnorm(1.100,2,1) + pnorm(2.461,2,1) - pnorm(1.300,2,1)
round(BETA, 3)

# Example 9.7
# part a.
library(PASWR)
xbar <- 56/30
cv <- qnorm(.95)
zobs <- (xbar - 1.8)/(2/sqrt(30))
pvalue <- 1-pnorm(zobs)
c(xbar,cv,zobs,pvalue)
# Or
zsum.test(mean.x=xbar,sigma.x=2,n.x=30,mu=1.8,alternative="greater")

# Example 9.8
# part a.
xbar <- 100/25
S <- sqrt((600-25*4^2)/24)
cvs <- qt(c(.025,.975),24)
cvs
tobs <- (xbar-2.5)/(S/sqrt(25))
pvalue <- 2*(1-pt(tobs,24))
c(xbar,cvs,tobs,pvalue)
tsum.test(mean.x=xbar,s.x=S,n.x=25,mu=2.5)
# part b.
delta <- (4-2.5)
lambda <- delta/(2.5/sqrt(25))
lambda
pt(qt(0.025,24),24,lambda)+(1-pt(qt(0.975,24),24,lambda))
# Or
1-pf(qt(.975,24)^2,1,24,lambda^2)
# Or
power.t.test(n=25,delta=1.5,sd=2.5,type="one.sample")
# part c.  Simulating the power of part b.
set.seed(13)
nvar <- rnorm(25 * 20000, 4, 2.5)
nvarmat <- matrix(nvar, 20000, 25) # 20000 by 25 Matrix
xbar <- apply(nvarmat, 1, mean)
S <- apply(nvarmat, 1, sd) # Change sd to stdev for S-PLUS
tstar <- (xbar - 2.5)/(S/5)
hist(tstar, xlim = c(-4, 8), nclass = "Scott", col = 13,
xlab = "", probability = T, ylim = c(0, 0.4),
main = "Central and Simulated Non-Central t-Distributions")
crit.tu <- qt(0.975, 24)
crit.tl <- qt(0.025, 24)
x <- seq(-4, 8, 0.05)
y <- dt(x, 24)
lines(x, y, lwd = 2)
lines(c(-4, 8), c(0, 0), lwd = 3)
lines(c(crit.tl, crit.tl), c(0, dt(crit.tl, 24)), lwd = 2)
lines(c(crit.tu, crit.tu), c(0, dt(crit.tu, 24)), lwd = 2)
SPU <- length(tstar[tstar > crit.tu])/length(tstar)
SPL <- length(tstar[tstar < crit.tl])/length(tstar)
Power <- SPU + SPL
Power

# Example 9.9
yield <- scan()
2.5 3.0 3.1 4.0 1.2 5.0 4.1 3.9 3.2 3.3 2.8 4.1 2.7 2.9 3.7

xbar <- mean(yield)
S <- sd(yield)       # S <- stdev(yield) for S-PLUS
cv <- qt(.95,14)
tobs <- (xbar - 2)/(S/sqrt(15))
pvalue <- 1-pt(tobs,14)
c(xbar, S, cv, tobs, pvalue)
# Or
t.test(yield, mu=2, alternative="greater")
rm("yield")

# Example 9.10
# part a.
zobs <- (20)/sqrt(100^2/64 + 108^2/144)
cv <- qnorm(.90)
pvalue <- 1-pnorm(zobs)
c(zobs,cv,pvalue)
#
# part b.
sig <- sqrt(100^2/64 + 108^2/144)
sig
CV <- qnorm(.90,0,sig)
CV
BETA <- pnorm(CV,40,sig)
BETA
POWER <- 1 - BETA
POWER

# Example 9.11
# part a.
attach(Stschool)
round(qt(0.975,24), 2) # Critical Value
mX <- mean(X,na.rm=TRUE)
mY <- mean(Y,na.rm=TRUE)
sX <- sd(X,na.rm=TRUE) # stdev(X,na.rm=T) for S-PLUS
sY <- sd(Y,na.rm=TRUE) # stdev(Y,na.rm=T) for S-PLUS
sp <- sqrt((10*sX^2 + 14*sY^2)/24) # Pooled stdev
tobs <- (mX - mY)/(sp*sqrt(1/11 + 1/15)) # t obs
round(c(mX, mY, sX, sY, sp, tobs), 2)
2*(1 - pt(tobs,24)) # P-value
#
t.test(X,Y, var.equal=TRUE)
# part b.
delta <- 10
sig <- sqrt(9^2/11 + 9^2/15)
lambda <- delta/sig
lambda
Lcv <- qt(.025,24)
Ucv <- qt(.975,24)
Lpower <- pt(Lcv,24,lambda)
Upower <- 1-pt(Ucv,24,lambda)
POWER <- Lpower + Upower
POWER
# Or
1-pf(qt(.975,24)^2,1,24,lambda^2)
detach(Stschool)

# Example 9.12
attach(Water)
X <- X[!is.na(X)]
Y <- Y[!is.na(Y)]
nX <- length(X); nY <- length(Y); mX <- mean(X)
mY <- mean(Y); sX2 <- var(X); sY2 <- var(Y)
nu <- (sX2/nX + sY2/nY)^2/((sX2/nX)^2/(nX-1) + (sY2/nY)^2/(nY-1))
round(qt(0.05, nu), 2)                   # Critical Value
tobs <- (mX - mY)/sqrt(sX2/nX + sY2/nY)  # t observed
round(c(mX, mY, sX2, sY2, nu, tobs), 2)
pt(tobs, nu)                             # P-value
#
# Or
t.test(X, Y, var.equal=FALSE, alternative="less")

# Example 9.13
#
qt(.95,9)         # Critical Value
library(lattice)  # Not needed for S-PLUS
attach(barley)
yieldMor32 <- yield[year == "1932" & site == "Morris"]
yieldCro32 <- yield[year == "1932" & site == "Crookston"]
d <- yieldMor32 - yieldCro32
EDA(d)            # Figure 9.10
t.test(d, alternative = "greater")
# OR
t.test(yieldMor32, yieldCro32, paired=TRUE, alternative="greater")
detach(barley)

# Example 9.14
#
qchisq(0.95,19) # Critical Value
attach(Washer)
s2 <- var(diameters)
s2
ChiObs <- 19*s2/0.004 # Standardized Test Statistic’s Value
ChiObs
1-pchisq(ChiObs,19) # P-value

# Example 9.15
#
attach(Bac)
EDA(X)
EDA(Y)
shapiro.test(Y)
# Hypothesis Test
qf(.05,9,9) # Critical Value
var(X)
var(Y)
Fobs <- var(X)/var(Y)
Fobs                  # Standardized Test Statistic’s Value
pf(Fobs,9,9)          # P-value
# Or
var.test(X, Y, alternative="less")

# Example 9.16
#
# part a.
# Likelihood method
probs <- dbinom(0:500,500,.2)
pvalue <- sum(probs[probs <= dbinom(90,500,.2)])
pvalue
binom.test(x=90,n=500,p=0.2)    # R
#
# Tail balancing method
pless <- pbinom(90,500,.2)
pgrea <- 1 - pbinom(89,500,.2)
c(pless,pgrea)
probs <- pbinom(0:500,500,.2)
pl <- pbinom(90,500,.2)
pl
pr <- (1-min(probs[probs > 1- pl]))
pr
pval <- pl + pr
pvals <- c(pl, pr, pval)
pvals
# part b. Confidence Interval using equation 9.8
yobs <- 90
n <- 500
alpha <- 0.05
FL <- qf(alpha/2,2*yobs,2*(n-yobs+1))
FU <- qf((1-alpha/2),2*(yobs+1),2*(n-yobs))
CI <- c( (1+(n-yobs+1)/(yobs*FL))^-1, (1+(n-yobs)/((yobs+1)*FU))^-1 )
CI

# Example 9.17
Y <- 90
n <- 500
p <- Y/n
PI <- 0.2
zobs <- (p - PI)/sqrt((PI * (1 - PI))/n)
pval <- 2 * pnorm(zobs)
zobsC <- (p - PI + 1/(2 * n))/sqrt((PI * (1 - PI))/n)
pvalC <- 2 * pnorm(zobsC)
round(c(zobs, pval, zobsC, pvalC), 4)
#
# part b. CI using equation 8.44 without continuity correction
alpha <- .05
z <- qnorm(1-alpha/2)
ll <- (p+z^2/(2*n) - z*sqrt(p*(1-p)/n +z^2/(4*n^2))) / (1+z^2/n)
ul <- (p+z^2/(2*n) + z*sqrt(p*(1-p)/n +z^2/(4*n^2))) / (1+z^2/n)
c(ll,ul)
# Or
prop.test(x=90,n=500, p=.2, correct=FALSE)
# CI using equation 8.44 with continuity correction
pl <- p - 1/(2*n)
pl
pu <- p + 1/(2*n)
pu
ll <- (pl+z^2/(2*n) - z*sqrt(pl*(1-pl)/n +z^2/(4*n^2))) / (1+z^2/n)
ul <- (pu+z^2/(2*n) + z*sqrt(pu*(1-pu)/n +z^2/(4*n^2))) / (1+z^2/n)
c(ll,ul)
prop.test(x=90,n=500, p=.2, correct=TRUE)

#
# Example 9.18
JD <- matrix(c(1,5,8,2),nrow=2,
dimnames=list(Youth=c("Delinquent","Non-delinquent"),
Glasses=c("Yes","No")))
JD # Output given for R, S-PLUS is slightly different
p <- dhyper(0:6,9,7,6)
# associated probabilities for 7 possible 2X2 tables
round(p,5)
pobs <- dhyper(1,9,7,6)
pval <- sum(p[p<=pobs])
pval
fisher.test(JD) # Output is for R, S-PLUS is slightly different.

#
# Example 9.19
HA <- matrix(c(104,189,10933,10845),nrow=2,
dimnames=list(Treatment=c("Aspirin","Placebo"),
Outcome=c("Heart attack","No heart attack")))
HA # Output is for R, S-PLUS is slightly different.
pval <- phyper(104,104+10933,189+10845,104+189) # x,m,n,k
pval
fisher.test(HA,alternative="less") # Only in R

#
# Example 9.20
# Without continuity correction
x <- 104
nx <- 11037
y <- 189
ny <- 11034
px <- x/nx
py <- y/ny
p <- (x+y)/(nx+ny)
zobs <- (px-py)/sqrt(p*(1-p)*(1/nx+1/ny))
zobs
prop.test(x=c(x,y),n=c(nx,ny),correct=FALSE)$stat^.5 * (-1)
prop.test(x=c(x,y),n=c(nx,ny),correct=FALSE)
# With continuity correction
zobs <- ((px-py)+(1/2*(1/nx+1/ny)))/sqrt(p*(1-p)*(1/nx+1/ny))
zobs
prop.test(x=c(x,y),n=c(nx,ny),correct=TRUE)$stat^.5 *(-1)
prop.test(x=c(x,y),n=c(nx,ny),correct=TRUE)