# This routine draws a two-sided hypothesis test -------------------------------

options(scipen = 6)                    # Stop sci notation if <6 dec places

# Provided Values
mean1 = 3.35                           # First sample mean
mean2 = 3.09                           # Second sample mean
sd1 = 1.24                             # First sample std deviation
sd2 = 1.21                             # Second sample std deviation
n1 = 142                               # First sample size
n2 = 126                               # Second sample size

# Calculated Values

seBlend = sqrt(sd1^2/n1 + sd2^2/n2)    # Blended std. error
dMeans = max(mean1, mean2) -           # Calculate difference in means 
  min(mean1, mean2)
z = (dMeans)/seBlend                   # Calculate z score
rangeZ = ceiling(abs(z)) + 2           # Set range for x values
rangeZ = ifelse(rangeZ < 4, 4, rangeZ) # Set minimum range of +/- 4

# Hypothesis test results
p = 2 * (1 - pnorm(z))                 # p value for 2 sided z-test
ptext = round(p, 2)                    # Version of p for annotation
OneMinusP = round(1 - p, 2)            # Version of 1 - p for annotation

# Graphing values
X1 = seq(-rangeZ, rangeZ, by = .001)   # Set up sequence of x values
Y1 = dnorm(X1)                         # Set up t curve values for each x
X2 = X1                                # Second x to use for filling curve
Y2 = Y1                                # Second y to use for filling curve
cd = data.frame(X1, Y1, X2, Y2)        # Join in data frame

png(filename = "illustrations/fig-15-14-normal-curve-BW.png",
  units = "in",                        # Set measurements in inches
  res = 1200,                          # Set resolution at 1200dpi
  width = 6,                           # Width at 6 inches
  height = 4)                          # Height at 4 inches

par(mai = c(1, 1, .25, .25))           # Set margin - no title


# Plot components
plot(cd$X1, cd$Y1,                     # Plot normal distribution
  type = "l", lwd = 2,                 # Line plot with linewidth =2
  ylim = c(0, .5),                     # Set range for y axis
  xlab = "z-score",                    # Add x axis label
  ylab = "Probability",                # Add y axis label
  yaxs = "i")                          # Eliminate y axis overage


abline(h = 0)                          # Horizontal boundary line at y = 0

lines(x = c(z, z), y = c(0, .43))      # Add short vertical line at z
lines(x = c(-z, -z), y = c(0, .43))    # Add short vertical line at -z    

# Set up borders of shaded regions
xz1 = which(
  cd$X1 == max(cd$X1[cd$X1 < z]))      # Find x just to left of z
cd$Y2[xz1] = 0                         # Set y for that x to 0      
cd$X2[xz1] = cd$X2[xz1 - 1]            # Shift max x over for vertical line
xz2 = which(
  cd$X1 == min(cd$X1[cd$X1 > -z]))     # Find x just to right of -z
cd$Y2[xz2] = 0                         # Set y for that x to 0      
cd$X2[xz2] = cd$X2[xz2 + 1]            # Shift max x over for vertical line

xz3 = which(
  cd$X1 == min(cd$X1[cd$X1 > z]))      # Find x just to right of z
cd$Y2[xz3] = 0                         # Set y for that x to 0      
cd$X2[xz3] = cd$X2[xz3 + 1]            # Shift max x over for vertical line

xz4 = which(
  cd$X1 == max(cd$X1[cd$X1 < -z]))     # Find x just to left of -z
cd$Y2[xz4] = 0                         # Set y for that x to 0      
cd$X2[xz4] = cd$X2[xz4 - 1]            # Shift max x over for vertical line
    
# Do shading
polygon(cd$X2[cd$X1 < z & cd$X1 > -z], # Fill for center region
  cd$Y2[cd$X1 < z & cd$X1 > -z],
  col = "light gray")                  # Set color
  
polygon(cd$X2[cd$X1 < -z],             # Fill for left tail 
  cd$Y2[cd$X1 < -z], 
  col = "darkgray")                    # Set color
  
polygon(cd$X2[cd$X1 > z],              # Fill for right tail 
  cd$Y2[cd$X1 > z],
  col = "darkgray")                    # Set color

abline(v = 0)                          # Vertical line at x = 0    

# Add labels
text(x = z, y = .42, pos = 4,          # Add label showing z-score 1
  labels = paste("z = +", round(z, 2)))
  
text(x = -z, y = .42, pos = 2,         # Add label showing z-score 2
  labels = paste("z = -", round(z, 2)))  

label1 = paste(                        # Annotation for area under curve
  "Area under curve\nbetween z values = ", 
  OneMinusP)
  
label2 = paste(                        # Annotation for area in tails
  "Area under curve\nin both tails = ",
  ptext)

text(x = -rangeZ/1.4, y = .2,          # Add label for main pr region
  labels = label1,                     # Use text created above
  cex = .75)                           # Font size 
  
arrows(x0 = -rangeZ/1.4, x1 = -z/2,    # Arrow from text to main pr region
  y0 = .17, y1 = .1,
  length = .15)                        # Set arrowhead size
  
text(x = rangeZ/1.4, y = .2,           # Add label for tails
  labels = label2,                     # Use label created above
  cex = .75)                           # Set text size

arrows(x0 = z + (rangeZ - z)/2,        # Add arrow from text to tails label
  y0 = .17,      
  y1 = .02,
  length = .15)
  
dev.off()                              # Output the png file