# Rosner, chapter 8 # ================================= # significance # from chapter 7 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 8.5 # H0 : delta == 0 # H1 : delta != 0 d <- c(13, 3, -1, 9, 7, 7, 6, 4, -2, 2) dBar <- mean(d) n <- length(d) s.d.squared <- sum((d - mean(d))^2 / (n-1)) s.d <- sqrt(s.d.squared) # paired t test : test statistic t <- dBar / (s.d / sqrt(n)); t # critical value method alpha <- .05 critical.value <- qt(alpha/2, df=n-1, lower.tail=FALSE); critical.value accept.H0 <- (-critical.value < t) & (t < critical.value); accept.H0 # p value p <- 2 * pt(t, df=n-1, lower.tail=FALSE); p significance(p) # ================================= # example 8.6 gynecology # H0 : delta == 0 # H1 : delta != 0 dBar <- 4 s.d <- 8 n <- 20 # paired t test : test statistic t <- dBar / (s.d / sqrt(n)); t # critical value method alpha <- .05 critical.value <- qt(alpha/2, df=n-1, lower.tail=FALSE); critical.value accept.H0 <- (-critical.value < t) & (t < critical.value); accept.H0 # p value p <- 2 * pt(t, df=n-1, lower.tail=FALSE); p significance(p) # ================================= # example 8.7 hypertension dBar <- 4.80 s.d <- 4.566 n <- 10 alpha <- .05 e <- qt(1 - alpha/2, df=n-1) * s.d / sqrt(n) ci <- dBar + c(-e, e); ci # ================================= # example 8.8 gynecology dBar <- 4 s.d <- 8 n <- 20 alpha <- .05 e <- qt(1 - alpha/2, df=n-1) * s.d / sqrt(n) ci <- dBar + c(-e, e); ci # ================================= # example 8.9 hypertension n.oc <- 8 mean.sbp.oc <- 132.86 s.oc <- 15.34 n.non.oc <- 21 mean.sbp.non.oc <- 127.44 s.non.oc <- 18.23 # SBP.oc ~ N(mu1, sigma1^2) # SBP.non.oc ~ N(mu2, sigma2^2) # H0 : mu1 == mu2 # H1 : mu1 != mu2 # assume sigma1 == sigma2 # ================================= # variance.pooled.estimate s.squared.pooled.estimate <- function(s1, s2, n1, n2) { ((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 -2) } # ================================= # t test, two independent samples, equal variances t.two.samples.equal.variances <- function(x1Bar, x2Bar, s, n1, n2){ (x1Bar - x2Bar) / (s * sqrt(1/n1 + 1/n2)) } # ================================= # example 8.10 hypertension # common variance # from example 8.9 n.oc <- 8 mean.sbp.oc <- 132.86 s.oc <- 15.34 n.non.oc <- 21 mean.sbp.non.oc <- 127.44 s.non.oc <- 18.23 s.squared <- s.squared.pooled.estimate(s1=s.oc, s2=s.non.oc, n1=n.oc, n2=n.non.oc) s.hypertension <- sqrt(s.squared); s.hypertension x1Bar <- mean.sbp.oc x2Bar <- mean.sbp.non.oc s <- s.hypertension n1 <- n.oc n2 <- n.non.oc t <- t.two.samples.equal.variances(x1Bar, x2Bar, s, n1, n2); t # critical value alpha <- .05 n1 <- n.oc n2 <- n.non.oc critical.value <- qt(1 - alpha/2, df=n1+n2-2, lower.tail=TRUE); critical.value accept.H0 <- (-critical.value < t) & (t < critical.value); accept.H0 # p value for t > 0 n1 <- n.oc n2 <- n.non.oc p <- 2 * pt(t, df=n1+n2-2, lower.tail=FALSE); p significance(p) # ================================= # CI for mean difference, two independent samples, equal variances ci.t.two.samples.equal.variances <- function(x1Bar, x2Bar, s, n1, n2, alpha){ e <- qt(1 - alpha/2, df=n1+n2-2) * s * sqrt(1/n1 + 1/n2); x1Bar - x2Bar + c(-e, e) } # ================================= # example 8.11 hypertension # from examples 8.9 and 8.10 x1Bar <- mean.sbp.oc x2Bar <- mean.sbp.non.oc s <- s.hypertension n1 <- n.oc n2 <- n.non.oc alpha <- .05 ci <- ci.t.two.samples.equal.variances(x1Bar, x2Bar, s, n1, n2, alpha); ci # ================================= # example 8.12 cardiovascular disease # H0 : var1 == var2 # H1 : va1 != var2 # X1 ~ N(mu1, sigma1^2) # X2 ~ N(mu2, sigma2^2) n1 <- 100 x1Bar <- 207.3 s1 <- 35.6 n2 <- 74 x2Bar <-193.4 s2 <- 17.3 ratio.of.sample.variances <- (s1 / s2)^2; ratio.of.sample.variances # ================================= # F distribution xs <- seq(from=0, to=2, by=.05) plot(xs, xs, type="n", xlim=c(0,2), ylim=c(0,1.5), xlab="x", ylab="f(x)", main="Density for F distributions") curve(df(x, df1=1, df2=6), col="dark red", add=TRUE) curve(df(x, df1=4, df2=6), col="rosybrown4", add=TRUE) # ================================= # example 8.13 cardiovascular disease df1 <- 5 df2 <- 9 alpha <- .01 qf(1 - alpha, df1, df2) # ================================= # example 8.14 F df1 <- 6 df2 <- 8 alpha <- .05 qf(alpha, df1, df2) 1 / qf(1 - alpha, df2, df1) # ================================= # example 8.15 cardiovascular disease, pediatrics # example 8.12 n1 <- 100 s1 <- 35.6 n2 <- 74 s2 <- 17.3 # test statistic f <- (s1 / s2)^2; f # critical values df1 <- n1 - 1 df2 <- n2 - 1 alpha <- .05 c1 <- qf(alpha/2, df1, df2); c1 c2 <- qf(1 - alpha/2, df1, df2); c2 accept.H0 <- (c1 < f) & (f < c2); accept.H0 # p value p <- 2 * pf(f, df1, df2, lower.tail=FALSE); p significance(p) # ================================= # example 8.16 hypertension s1 <- 15.34 s2 <- 18.23 n1 <- 8 n2 <- 21 # variance ratio f <- (s2 / s1)^2; f df1 <- n1 - 1 df2 <- n2 - 1 alpha <- .05 c1 <- qf(alpha/2, df2, df1); c1 c2 <- qf(1 - alpha/2, df2, df1); c2 accept.H0 <- (c1 < f) & (f < c2); accept.H0 # p value p <- 2 * pf(f, df2, df1, lower.tail=FALSE); p significance(p) # ================================= # t test, two independent samples, unequal variances t.two.samples.unequal.variances <- function(x1Bar, x2Bar, s1, s2, n1, n2){ (x1Bar - x2Bar) / sqrt(s1^2 / n1 + s2^2 / n2) } # ================================= # df.approximate d.prime <- (s1^2 / n1 + s2^2 / n2)^2 / ( (s1^2 / n1)^2 / (n1 - 1) + (s2^2 / n2)^2 / (n2 - 1) ) df.approximate <- floor(d.prime) # ================================= # example 8.17 cardiovascular disease, pediatrics # from example 8.12 and 8.15 # H0 : var1 == var2 # H1 : va1 != var2 # X1 ~ N(mu1, sigma1^2) # X2 ~ N(mu2, sigma2^2) n1 <- 100 x1Bar <- 207.3 s1 <- 35.6 n2 <- 74 x2Bar <-193.4 s2 <- 17.3 t <- t.two.samples.unequal.variances(x1Bar, x2Bar, s1, s2, n1, n2); t # df.approximate d.prime <- (s1^2 / n1 + s2^2 / n2)^2 / ( (s1^2 / n1)^2 / (n1 - 1) + (s2^2 / n2)^2 / (n2 - 1) ) df.approximate <- floor(d.prime); df.approximate # critical values alpha <- .05 c1 <- qt(alpha/2, df=df.approximate); c1 c2 <- qt(1 - alpha/2, df=df.approximate); c2 accept.H0 <- (c1 < t) & (t < c2); accept.H0 # p value p <- 2 * pt(t, df=df.approximate, lower.tail=FALSE); p significance(p) # ================================= # example 8.18 infectious disease n1 <- 7 x1Bar <- 11.57 s1 <- 8.81 n2 <- 18 x2Bar <-7.44 s2 <- 3.70 # F test for equality of two variances f <- (s1 / s2)^2; f # critical values df1 <- n1 - 1 df2 <- n2 - 1 alpha <- .05 c1 <- qf(alpha/2, df1, df2); c1 c2 <- qf(1 - alpha/2, df1, df2); c2 accept.H0 <- (c1 < f) & (f < c2); accept.H0 # p value p <- 2 * pf(f, df1, df2, lower.tail=FALSE); p significance(p) # two sample t test with unequal variances t <- t.two.samples.unequal.variances(x1Bar, x2Bar, s1, s2, n1, n2); t # df.approximate d.prime <- (s1^2 / n1 + s2^2 / n2)^2 / ( (s1^2 / n1)^2 / (n1 - 1) + (s2^2 / n2)^2 / (n2 - 1) ) df.approximate <- floor(d.prime); df.approximate # p value p <- 2 * pt(t, df=df.approximate, lower.tail=FALSE); p significance(p) # ================================= # two sided CI for mu1 - mu2 with unequal variances ci.t.two.samples.unequal.variances <- function(x1Bar, x2Bar, s1, s2 , n1, n2, df.approximate, alpha){ e <- qt(1 - alpha/2, df= df.approximate) * sqrt( s1^2/n1 + s2^2/n2); x1Bar - x2Bar +c(-e, e) } # ================================= # example 8.19 infectious disease # from example 8.18 (with slightly different numbers) n1 <- 7 x1Bar <- 11.571 s1 <- 8.810 n2 <- 18 x2Bar <-7.444 s2 <- 3.698 alpha <- .05 df.approximate <- 6 ci <- ci.t.two.samples.unequal.variances(x1Bar, x2Bar, s1, s2 , n1, n2, df.approximate, alpha); ci # ================================= # sample size (for each group) sample.size.two.sided <- function(mu1, sigma1, mu2, sigma2, alpha, beta){ (sigma1^2 + sigma2^2) * (qnorm(1-alpha/2) + qnorm(1-beta))^2 / abs(mu2 - mu1)^2 } # ================================= # sample sizes (for each group : n2 == k * n1) sample.sizes.two.sided <- function(mu1, sigma1, mu2, sigma2, alpha, beta, k){ size1 <- (sigma1^2 + (sigma2)^2/k) * (qnorm(1-alpha/2) + qnorm(1-beta))^2 / abs(mu2 - mu1)^2; size2 <- (k*sigma1^2 + sigma2^2) * (qnorm(1-alpha/2) + qnorm(1-beta))^2 / abs(mu2 - mu1)^2; c(size1, size2) } # ================================= # example 8.28 and 8.29 hypertension # from example 8.9. p.276 # X1 ~ N(mu1, sigma1^2) # X2 ~ N(mu2, sigma2^2) # H0 : mu1 == mu2 # H1 : mu1 != mu2 mu1 <- 132.86 s1 <- 15.34 mu2 <- 127.44 s2 <- 18.23 alpha <- .05 beta <- .20 sample.size <- sample.size.two.sided(mu1, sigma1=s1, mu2, sigma2=s2, alpha, beta); sample.size ceiling(sample.size) # ================================= # example 8.30 hypertension # data of examples 8.28 and 8.29 k <- 2 sample.sizes.two.sided(mu1, sigma1=s1, mu2, sigma2=s2, alpha, beta, k) # ================================= # example 8.31 hypertension # X1 ~ N(mu1, sigma1^2) # X2 ~ N(mu2, sigma2^2) # H0 : mu1 == mu 2 # H1 : mu1 != mu 2 # anticipated : abs(mu1 - mu2) == delta # ================================= # power power <- function(sigma1, sigma2, n1, n2, alpha, delta){ z <- -qnorm(1-alpha/2) + delta / sqrt( sigma1^2/n1 + sigma2^2/n2 ); pnorm(z) } # ================================= # example 8.32 hypertension sigma1 <- 15.34 sigma2 <- 18.23 n1 <- 100 n2 <- 100 alpha <- .05 delta <- 5 power(sigma1, sigma2, n1, n2, alpha, delta) # ================================= # sample size, longitudinal study sample.size.longitudinal.two.sided <- function(sigma1, sigma2, alpha, beta, delta, rho){ var.d <- sigma1^2 + sigma2^2 - 2 * rho * sigma1 * sigma2; 2 * var.d * (qnorm(1-alpha/2) + qnorm(1-beta))^2 / delta^2 } # ================================= # power, longitudinal study power.longitudinal.two.sided <- function(sigma1, sigma2, n, alpha, delta, rho){ var.d <- sigma1^2 + sigma2^2 - 2 * rho * sigma1 * sigma2; z <- -qnorm(1-alpha/2) + sqrt(n) * delta / sqrt(2 * var.d); pnorm(z) } # ================================= # example 8.33 hypertension # H0 : mu1 == mu 2 # H1 : mu1 != mu 2 # anticipated : abs(mu1 - mu2) == delta sigma1 <- 15 sigma2 <- 15 mu1 <- -8 mu2 <- -3 alpha <- .05 beta <- .20 delta <- abs(mu1 - mu2) rho <- .70 sample.size <- sample.size.longitudinal.two.sided(sigma1, sigma2, alpha, beta, delta, rho); sample.size ceiling(sample.size) # ================================= # example 8.34 hypertension n <- 75 power <- power.longitudinal.two.sided(sigma1, sigma2, n, alpha, delta, rho); power