# Rosner, chapter 7 # ================================= # example 7.10 obstetrics # H0 : mu == 120 # H1 : mu < 120 n <- 100 xBar <- 115 s <- 24 mu0 <- 120 t <- (xBar - mu0)/(s/sqrt(n)); t # critical value alpha <- .05 critical.value <- qt(p=alpha, df=n-1, lower.tail=TRUE); critical.value accept.H0 <- t >= critical.value; accept.H0 # rejection region xs <- seq(from=-3, to=3, by=.05) plot(xs, dt(xs, df=n-1), type="l", col="dark red", xlab="t", ylab="density", main="t ~ t(df=99) : rejection region") lines(x=c(critical.value, critical.value), y=c(0, dt(critical.value,df=99)), col="red") poly.xs <- seq(from=-3, to= critical.value, by=.01) polygon(x=c(-3, poly.xs, critical.value), y=c(0, dt(poly.xs, df=n-1), 0), col="rosybrown2", border="dark red") # ================================= # example 7.11 obstetrics # critical value alpha <- .01 n <- 100 critical.value <- qt(p=alpha, df=n-1, lower.tail=TRUE); critical.value # ================================= # example 7.12 obstetrics # p value critical.value <- -2.08 n <- 100 pt(critical.value, df=n-1, lower.tail=TRUE) # ================================= # significance significance <- function(p){ if (.01 <= p & p < .05) { print('results.significant') } else if (.001 <= p & p < .01) { print('results.highly.significant') } else if (p < .001) { print('results.very.highly.significant') } else { print('results.not.statistically.significant') } } # ================================= # example 7.13 cardiology # H0 : mu == 25 # H1 : mu < 25 # p value n <- 8 xBar <- 16 s <- 10 mu0 <- 25 t.statistic <- (xBar - mu0)/(s/sqrt(n)); t p <- pt(t.statistic, df=n-1, lower.tail=TRUE); p # result alpha <- .05 accept.H0 <- p >= alpha; accept.H0 significance(p) # ================================= # example 7.14 obstetrics significance(p=.020) # ================================= # example 7.15 cardiology significance(p=.019) # ================================= # example 7.16 obstetrics # H0 : mu == 120 # H1 : mu < 120 n <- 10000 xBar <- 119 s <- 24 mu0 <- 120 t <- (xBar - mu0)/(s/sqrt(n)); t # critical value alpha <- .05 critical.value <- qt(p=alpha, df=n-1, lower.tail=TRUE); critical.value accept.H0 <- t >= critical.value; accept.H0 # normal approximation p <- pnorm(t); p significance(p) # however, xBar == 119 is not very different from mu == 120 # ================================= # example 7.17 obstetrics # H0 : mu == 120 # H1 : mu < 120 n <- 10 xBar <- 110 s <- 24 mu0 <- 120 t <- (xBar - mu0)/(s/sqrt(n)); t # critical value alpha <- .05 critical.value <- qt(p=alpha, df=n-1, lower.tail=TRUE); critical.value accept.H0 <- t >= critical.value; accept.H0 # p value p <- pt(t, df=n-1, lower.tail=TRUE); p significance(p) # ================================= # example 7.18 cardiovascular disease, pediatrics # H0 : mu == 175 # H1 : mu > 175 n <- 10 xBar <- 200 s <- 50 mu0 <- 175 t <- (xBar - mu0)/(s/sqrt(n)); t # critical value alpha <- .05 critical.value <- qt(p=1-alpha, df=n-1); critical.value accept.H0 <- t < critical.value; accept.H0 # p value p <- 1 - pt(t, df=n-1); p significance(p) # rejection region xs <- seq(from=-3, to=3, by=.05) plot(xs, dt(xs, df=n-1), type="l", col="dark red", xlab="t (score == 1.58)", ylab="density", main="t ~ t(df=9) : rejection region") lines(x=c(critical.value, critical.value), y=c(0, dt(critical.value,df=9)), col="red") poly.xs <- seq(from= critical.value, to=3, by=.01) polygon(x=c(critical.value, poly.xs, 3), y=c(0, dt(poly.xs, df=n-1), 0), col="rosybrown2", border="dark red") lines(x=c(t, t), y=c(0, 0.04), lwd=2, col="blue") # ================================= # example 7.20 cardiovascular disease # H0 : mu == mu0 == 190 # H1 : mu != mu0 mu0 <- 190 n <- 100 mu <- 181.52 s <- 40 # ================================= # example 7.21 cardiovascular disease # H0 : mu == mu0 == 190 # H1 : mu != mu0 mu0 <- 190 n <- 100 xBar <- 181.52 s <- 40 t <- (xBar - mu0)/(s/sqrt(n)); t # critical values alpha <- .05 critical.value.1 <- qt(p=alpha/2, df=n-1); critical.value.1 critical.value.2 <- qt(p=1-alpha/2, df=n-1); critical.value.2 accept.H0 <- (critical.value.1 < t & t < critical.value.2); accept.H0 # ================================= # example 7.22 cardiovascular disease # p value p <- 2 * pt(t, df=n-1, lower.tail=TRUE); p significance(p) # ================================= # example 7.24 cardiovascular disease # H0 : mu == mu0 # H1 : mu != mu0 # sigma known xBar <- 181.52 mu0 <- 190 sigma <- 40 n <- 200 z <- (xBar - mu0)/(sigma/sqrt(n)); z # critical values alpha <- .05 critical.value.1 <- qnorm(p=alpha/2); critical.value.1 critical.value.2 <- qnorm(p=1-alpha/2); critical.value.2 accept.H0 <- (critical.value.1 < z & z < critical.value.2); accept.H0 # p value p <- 2 * pnorm(z, lower.tail=TRUE); p # ================================= # power (mu1 < mu0) # H0 : mu1 == mu0 # H1 : mu1 < mu0 power.lower.tail <- function(mu0, mu1, sigma, n, alpha) { pnorm(qnorm(alpha, lower.tail=TRUE) + (mu0 - mu1) * sqrt(n) /sigma) } # ================================= # example 7.26 obstetrics # H0 : mu1 == mu0 # H1 : mu1 < mu0 mu0 <- 120 mu1 <- 115 sigma <- 24 n <- 100 alpha <- .05 power.lower.tail(mu0, mu1, sigma, n, alpha) # ================================= # power (mu1 > mu0) # H0 : mu1 == mu0 # H1 : mu1 > mu0 power.upper.tail <- function(mu0, mu1, sigma, n, alpha) { 1 - pnorm(qnorm(1 - alpha, lower.tail=TRUE) + (mu0 - mu1) * sqrt(n) /sigma) } # ================================= # example 7.27 cardiovascular disease, pediatrics # H0 : mu1 == mu0 # H1 : mu1 > mu0 mu0 <- 175 mu1 <- 190 sigma <- 50 n <- 10 alpha <- .05 power.upper.tail(mu0, mu1, sigma, n, alpha) # ================================= # example 7.28 cardiovascular disease, pediatrics # H0 : mu1 == mu0 # H1 : mu1 > mu0 mu0 <- 175 mu1 <- 190 sigma <- 50 n <- 10 alpha <- .01 power.upper.tail(mu0, mu1, sigma, n, alpha) # ================================= # example 7.29 obstetrics # H0 : mu1 == mu0 # H1 : mu1 < mu0 mu0 <- 120 mu1 <- 110 sigma <- 24 n <- 100 alpha <- .05 power.lower.tail(mu0, mu1, sigma, n, alpha) # ================================= # example 7.30 cardiology # H0 : mu == mu0 # H1 : mu < mu0 mu0 <- 25 mu1 <- 20 sigma <- 10 n <- 8 alpha <- .05 power.lower.tail(mu0, mu1, sigma, n, alpha) mu0 <- 25 mu1 <- 20 sigma <- 15 n <- 8 alpha <- .05 power.lower.tail(mu0, mu1, sigma, n, alpha) # ================================= # example 7.31 obstetrics # H0 : mu1 == mu0 # H1 : mu1 < mu0 mu0 <- 120 mu1 <- 115 sigma <- 24 n <- 10 alpha <- .05 power.lower.tail(mu0, mu1, sigma, n, alpha) # ================================= # power curve mu0 <- 120 sigma <- 24 n <- 100 alpha <- .05 mus <- seq(from=106, to=120, by=.1) plot(mus, power.lower.tail(mu0, mus, sigma, n, alpha), type="l", col="dark red",xlab="birthweight", ylab="power", main="Power curve for birthweight data") # ================================= # power two-tailed # H0 : mu1 == mu0 # H1 : mu1 != mu0 power.two.tailed <- function(mu0, mu1, sigma, n, alpha) { pnorm(-qnorm(1 - alpha/2, lower.tail=TRUE) + (mu0 - mu1) * sqrt(n) /sigma) + pnorm(-qnorm(1 - alpha/2, lower.tail=TRUE) + (mu1 - mu0) * sqrt(n) /sigma) } power.two.tailed.approximation <- function(mu0, mu1, sigma, n, alpha) { pnorm(-qnorm(1 - alpha/2, lower.tail=TRUE) + abs(mu0 - mu1) * sqrt(n) /sigma) } # ================================= # example 7.32 cardiology # abs(mu0 - mu1) == 5 sigma <- 10 n <- 20 alpha <- .05 # approximation pnorm(-qnorm(1 - alpha/2, lower.tail=TRUE) + 5 * sqrt(n) /sigma) # true value pnorm(-qnorm(1 - alpha/2, lower.tail=TRUE) + 5 * sqrt(n) /sigma) + pnorm(-qnorm(1 - alpha/2, lower.tail=TRUE) - 5 * sqrt(n) /sigma) # ================================= # sample size estimation (one-sided alternative) # H0 : mu == mu0 # H1 : mu == mu1 sample.size.one.sided <- function(mu0, mu1, sigma, alpha, beta) { sigma^2 * (qnorm(1-beta) + qnorm(1-alpha))^2 / (mu0 - mu1)^2 } # ================================= # example 7.33 obstretics mu0 <- 120 mu1 <- 115 sigma <- 24 alpha <- .05 beta <- .20 sample.size.one.sided(mu0, mu1, sigma, alpha, beta) # ================================= # example 7.34 obstretics mu0 <- 120 mu1 <- 110 mu1.old <- 115 (mu0 - mu1.old)^2 / (mu0 - mu1)^2 # sample size would be 1/4 as large # ================================= # example 7.35 carsiovascular disease, pediatrics mu0 <- 175 mu1 <- 190 sigma <- 50 alpha <- .05 beta <- .10 sample.size.one.sided(mu0, mu1, sigma, alpha, beta) # ================================= # example 7.36 obstretics mu0 <- 120 mu1 <- 115 sigma <- 30 # instead of 24 alpha <- .05 beta <- .20 sample.size.one.sided(mu0, mu1, sigma, alpha, beta) mu0 <- 120 mu1 <- 115 sigma <- 24 alpha <- .01 # instead of .05 beta <- .20 sample.size.one.sided(mu0, mu1, sigma, alpha, beta) mu0 <- 120 mu1 <- 115 sigma <- 24 alpha <- .05 beta <- .10 # instead of .20 sample.size.one.sided(mu0, mu1, sigma, alpha, beta) mu0 <- 120 mu1 <- 110 # instead of 115 sigma <- 24 alpha <- .05 beta <- .20 sample.size.one.sided(mu0, mu1, sigma, alpha, beta) # ================================= # sample size estimation (two-sided alternative) # H0 : mu == mu0 # H1 : mu == mu1 sample.size.two.sided <- function(mu0, mu1, sigma, alpha, beta) { sigma^2 * (qnorm(1-beta) + qnorm(1-alpha/2))^2 / (mu0 - mu1)^2 } # ================================= # example 7.37 cardiology # abs(mu0 - mu1) == 5 sigma <- 10 alpha <- .05 beta <- .20 # two-sided test sigma^2 * (qnorm(1-beta) + qnorm(1-alpha/2))^2 / 5^2 # one-sided test sigma^2 * (qnorm(1-beta) + qnorm(1-alpha))^2 / 5^2 # ================================= # example 7.38 cardiology # sample size estimation based on CI width sample.size.CI <- function(s, L, alpha){ 4 * qnorm(1 - alpha/2)^2 * s^2 / L^2 } # ================================= # example 7.39 cardiology s <- 10 L <- 5 alpha <- .05 sample.size.CI(s, L, alpha) # ================================= # two-sided CI for mu using t ci.two.sided.t <- function(xBar, s, n, alpha){ e <- qt(p=(1 - alpha/2), df=n-1) * s / sqrt(n) xBar + c(-e, e) } # ================================= # example 7.40 cardiology # H0 : mu == mu0 == 190 # H1 : mu != mu0 xBar <- 181.52 s <- 40 n <- 100 alpha <- .05 ci.two.sided.t(xBar, s, n, alpha) # ================================= # example 7.41 cardiovascular disease # H0 : mu == mu0 # H1 : mu != mu0 mu0 <- 190 xBar <- 185 s <- 40 n <- 100 alpha <- .05 ci.two.sided.t(xBar, s, n, alpha) # this CI contains the null mean, mu0 == 190 # p value p <- 2 * pt((xBar - mu0) / (s/sqrt(n)), df=n-1); p # p > .05 # ================================= # example 7.42 cardiovascular disease # from 7.22 mu0 <- 190 xBar <- 181.52 s <- 40 n <- 100 t <- (xBar - mu0)/(s/sqrt(n)); t # p value p <- 2 * pt(t, df=n-1, lower.tail=TRUE); p # from 7.40 xBar <- 181.52 s <- 40 n <- 100 alpha <- .05 ci.two.sided.t(xBar, s, n, alpha) # ================================= # two-sided CI for mu using z ci.two.sided.z <- function(xBar, sigma, n, alpha){ e <- qnorm(p=(1 - alpha/2)) * sigma / sqrt(n) xBar + c(-e, e) } # ================================= # example 7.43 cardiovascular disease # H0 : mu == mu0 == 190 # H1 : mu != mu0 mu0 <- 190 xBar <- 181.52 sigma <- 40 n <- 200 alpha <- .05 ci.two.sided.z(xBar, sigma, n, alpha) # mu ~ N(xBar, sigma^2/n) # P(mu < 190) pnorm((mu0 - xBar)/(sigma/sqrt(n))) # ================================= # example 7.45 hypertension sigma0Squared <- 35 sSquared <- 8.178 n <- 10 alpha <- .05 Xsquared <- (n-1) * sSquared / sigma0Squared; Xsquared critical.values <- c(qchisq(alpha/2, df=n-1), qchisq(1 - alpha/2, df=n-1)); critical.values accept.H0 <- (critical.values[1] < Xsquared) & (Xsquared < critical.values[2]); accept.H0 # p value p <- 2 * pchisq(Xsquared, df=n-1); p # ================================= # example 7.48 cancer pHat <- .04 p0 <- .02 q0 <- 1 - p0 n <- 10000 alpha <- .05 # test n * p0 * q0 >= 5 z <- (pHat - p0) /sqrt(p0 * q0 / n); z ci <- c(qnorm(alpha/2), qnorm(1 - alpha/2)); ci accept.H0 <- (ci[1] < z) & (z < ci[2]); accept.H0 # p value # pHat > p0 p <- 2 * (1 - pnorm(z)); p significance(p) # ================================= # example 7.49 occupational health, cancer # H0 : p == .20 # H1 : p != .20 pHat <- 5/13 p0 <- .20 q0 <- 1 - p0 n <- 13 alpha <- .05 # test n * p0 * q0 >= 5 pHat > p0 p <- 2 * (1 - sum(dbinom(0:4, size=n, prob=p0))); p significance(p) # ================================= # power for one-sample binomial test (two-sided) power.binomial.two.sided <- function(p0, q0, p1, q1, n, alpha){ pnorm(sqrt((p0 * q0) / (p1 * q1)) * (qnorm(alpha/2) + abs(p0 - p1) * sqrt(n) / sqrt(p0 * q0))) } # assuming n*p0*q1 >= 5 # ================================= # example 7.50 cancer p0 <- .02 q0 <- 1 - p0 p1 <- .05 q1 <- 1 - p1 n <- 500 alpha <- .05 # test n * p0 * q1 >= 5 power.binomial.two.sided(p0, q0, p1, q1, n, alpha) # ================================= # sample size for one-sample binomial test (two-sided) # H0 : p == p0 # H1 : p != p0 n.binomial.two.sided <- function(p0, q0, p1, q1, n, alpha, beta){ (p0 * q0) * (qnorm(1 - alpha/2) + qnorm(1 - beta) * sqrt((p1 * q1) / (p0 * q0)))^2 / (p1 - p0)^2 } # ================================= # example 7.51 cancer p0 <- .02 q0 <- 1 - p0 p1 <- .05 q1 <- 1 - p1 n <- 500 alpha <- .05 beta <- .10 n.binomial.two.sided(p0, q0, p1, q1, n, alpha, beta) # ================================= # example 7.52 and 7.53 cancer # H0 : mu == 3.3 # H1 : mu != 3.3 x <- 4 # use Table 8, p.827 # from: example 6.57 exact CI for Poisson parameter # load package epitools library(epitools) pois.exact(x=4, pt=8418*10, conf.level=.95) # => # x pt rate lower upper conf.level # 1 4 84180 4.751722e-05 1.294713e-05 0.000121663 0.95 # values for lambda: # 1.294713e-05 and 0.000121663 events per 8418*10 person-years is 95% CI for lambda # corresponding values for mu == lambda * T: 1.294713e-05 * 8418 * 10 0.000121663* 8418 * 10 # 1.089889 and 10.24159 # these are exactly the values for x == 4 alpha == .95 in table 8, p.827 # there is also a poisson.test() in the stats package accept.H0 <- (1.089889 < 4) & (4 < 10.24159); accept.H0 # ================================= # example 7.54 cancer # p value method x <- 4 mu0 <- 3.3 p <- 2 * (1 - ppois(3, lambda=mu0)); p significance(p) # ================================= # example 7.55 occupational health smr <- 100 * 4 / 3.3; smr # ================================= # example 7.56 occupational health library(epitools) x <- 21 mu0 <- 18.1 pt <- 8418 * 10 pois.exact(x=21, pt=8418*10, conf.level=.95) # x pt rate lower upper conf.level # 1 21 84180 0.0002494654 0.0001544233 0.0003813344 0.95 ci.lower <- 0.0001544233 * pt; ci.lower ci.upper <- 0.0003813344 * pt; ci.upper accept.H0 <- (ci.lower < x) & (x < ci.upper); accept.H0 # p value p <- 2 * (1 - ppois(20, lambda=mu0)); p significance(p) smr <- 100 * x / mu0; smr # ================================= # example 7.57 occupational health # large sample test # H0 : mu == mu0 == 18.1 # H1 : mu != mu0 x <- 21 XSquared <- (x - mu0)^2 / mu0; XSquared # XSquared ~ chi_1^2 under H0 alpha <-.05 critical.value <- qchisq(p=(1 - alpha), df=1); critical.value accept.H0 <- (XSquared < critical.value); accept.H0 # p value p <- 1 - pchisq(XSquared, df=1); p # ================================= # example 7.58 endocrinology # H0 : mu == mu0 == 0 # H1 : mu != mu0 mu0 <- 0 xBar <- -5.0 n <- 41 se <- 2.0 alpha <- .05 t <- (xBar - mu0) / se; t # t ~ t(df=n-1) p <- 2 * pt(t, df=n-1); p significance(p)