################################
## p-value for nutrition study example
1 - pt(3.01, df = 15.99)
x=seq(-3.5,3.5,by=0.01)
plot(x,dt(x,9),type="l",xaxt="n",xlab="",ylab="",yaxt="n")
polygon(c(qt(0.975,9),seq(qt(0.975,9),4, by=0.01),4, qt(0.975,9)),c(0,dt(seq(qt(0.975,9),4, by=0.01),9),0, 0),col="pink")
polygon(c(-qt(0.975,9),seq(-qt(0.975,9),-4, by=-0.01),-4,-qt(0.975,9)),c(0,dt(seq(-qt(0.975,9),-4, by=-0.01),9),0,0),col="pink")
text(0,0.1,"Acceptance",cex=2)
text(-3,0.05,"Rejection",cex=2)
text(3,0.05,"Rejection",cex=2)
axis(1, at=c(0), labels=expression(bar(x)-bar(y)), cex.axis=2)


################################
## Read the two-sample in R
xA=c(7.53, 7.48, 8.08, 8.09, 10.15, 8.4, 10.88, 6.13, 7.9)
xB=c(9.21, 11.51, 12.79, 11.85, 9.97, 8.79, 9.69, 9.68, 9.19)
## A two sample Welch t-test
t.test(xB, xA)


################################
## Read the two samples
x1=c(.7,-1.6,-.2,-1.2,-1,3.4,3.7,.8,0,2)
x2=c(1.9,.8,1.1,.1,-.1,4.4,5.5,1.6,4.6,3.4)
## Take the difference to get a paired t-test
dif=x2-x1
## Calculate the test and results
t.test(dif)
## Another way to calculate the paired setup
t.test(x2, x1, paired=TRUE)


################################
## Q-Q plot for each sample
par(mfrow=c(1,2))
qqnorm(xA, main="Hospital A")
qqline(xA)
qqnorm(xB, main="Hospital B")
qqline(xB)

################################
## Multiple (simulated) Q-Q plots for each sample
require(MESS)
fit1 <- lm(xA ~ 1)
qqnorm.wally <- function(x, y, ...) { qqnorm(y, ...); qqline(y, ...)}
wallyplot(fit1, FUN=qqnorm.wally, main="")
## Multiple (simulated) Q-Q plots for each sample
fit1 <- lm(xB ~ 1)
qqnorm.wally <- function(x, y, ...) { qqnorm(y, ...); qqline(y, ...)}
wallyplot(fit1, FUN=qqnorm.wally, main="")
power.t.test(n = 40, delta = 4, sd = 12.21, 
             type = "one.sample")
power.t.test(power = .80, delta = 4, sd = 12.21, 
             type = "one.sample")


################################
## Power calculation for two-sample
power.t.test(n = 10, delta = 2, sd = 1, sig.level = 0.05)


################################
## Sample size calculation for two-sample
power.t.test(power = 0.90, delta = 2, sd = 1, sig.level = 0.05)


################################
## Detectable effect size calculation for two-sample
power.t.test(power = 0.90, n = 10, sd = 1, sig.level = 0.05)