# Rosner, chapter 6 # ================================= # example 6.12 hypertension # random selection sample(0:999, 20) # ================================= # example 6.15 diabetes # random selection sample(0:9, 5) # ================================= # example 6.16 obstetrics # repeat 5 times: sample(0:999, 10) # ================================= # example 6.24 gynecology data <-c(97.2, 96.8, 97.4, 97.4, 97.3, 97.0, 97.1, 97.3, 97.2, 97.3) mean(data) sd(data) length(data) se <- sd(data) / sqrt(length(data)); se ci <- mean(data) + se * c(-2,2); ci # ================================= # example 6.30 X ~ qt(p=.95, df=23) qt(p=.95, df=23) # ================================= # example 6.31 CI for mean birthweight # sample 1 xBar <- 116.90 s <- 21.7 n <- 10 alpha <- .05 # 95% CI tstar <- qt(p=1-alpha/2, df=n-1) sem <- s / sqrt(n) # standard error of the mean e <- tstar * sem ci <- xBar + c(-e, e); ci # ================================= # CI for the mean using t distribution t.ci <- function(xBar, s, n, alpha){ tstar <- qt(p=1-alpha/2, df=n-1); sem <- s / sqrt(n); e <- tstar * sem; xBar + c(-e, e)} # ================================= # example 6.32 CI for mean birthweight t.ci(xBar=116.90, s=21.70, n=10, alpha=.05) # sample 1 t.ci(xBar=132.80, s=32.62, n=10, alpha=.05) # sample 2 t.ci(xBar=117.00, s=22.44, n=10, alpha=.05) # sample 3 t.ci(xBar=106.70, s=14.13, n=10, alpha=.05) # sample 4 t.ci(xBar=111.90, s=20.46, n=10, alpha=.05) # sample 5 # ================================= # example 6.33 gynecology data <-c(97.2, 96.8, 97.4, 97.4, 97.3, 97.0, 97.1, 97.3, 97.2, 97.3) mean(data) sd(data) length(data) t.ci(xBar=97.20, s=.189, n=10, alpha=.05) # ================================= # example 6.34 CI for mean birthweight t.ci(xBar=116.90, s=21.70, n=10, alpha=.01) # sample 1 # ================================= # example 6.35 gynecology t.ci(xBar=97.20, s=.189, n=100, alpha=.05) # ================================= # example 6.36 CI for mean basal temperature t.ci(xBar=97.2, s=0.4, n=10, alpha=.05) # ================================= # example 6.37 cardiovascular diseases, pediatrics t.ci(xBar=207.3, s=30, n=100, alpha=.05) # ================================= # case study 6.6 : BMD # sem == s/sqrt(n) sem <- 0.014 n <- 41 t.ci(xBar=-0.036, s=(sem * sqrt(41)), n=41, alpha=.05) # ================================= # example 6.38 gynecology data <-c(97.2, 96.8, 97.4, 97.4, 97.3, 97.0, 97.1, 97.3, 97.2, 97.3) sd(data)^2 # ================================= # example 6.39 hypertension data <- c(-6, 3, 2, -3, 1, 0, -1, 1, 3, -2) mean(data) sd(data)^2 # chi-square distributions plot(1, type="n", xlim=c(0,10), ylim=c(0, 1.3), xlab="value", ylab="frequency", main="chi square distributions, df=1,2,5") curve(dchisq(x, df=1), type="l", col="steelblue", add=TRUE) curve(dchisq(x, df=2), type="l", col="steelblue1", add=TRUE) curve(dchisq(x, df=5), type="l", col="steelblue4", add=TRUE) # ================================= # example 6.40 hypertension alpha <- .05 qchisq((1-alpha/2), df=10) qchisq(alpha/2, df=10) # ================================= # CI for the variance using chi square distribution chiSquare.ci <- function(sSquared, n, alpha){ uLower <- qchisq(p=alpha/2, df=n-1); uUpper <- qchisq(p=1-alpha/2, df=n-1); c((n-1) * sSquared / uUpper, (n-1) * sSquared / uLower)} # ================================= # example 6.41 hypertension chiSquare.ci(sSquared=8.178, n=10, alpha=.05) # ================================= # example 6.42 and 6.43 cancer n <- 5000 p <- 28 pHat <- p / n; pHat se.estimated <- sqrt(pHat * (1-pHat) / n); se.estimated # ================================= # example 6.44 cancer # liklihood == p^30 * (1-p)^70 # ================================= # example 6.46 MLE, cancer # log(liklihood) == 30 * log(p) + 70 log(1-p) # d/dp log(liklihood) == 0 # 30/p == 70/(1-p) # 30*(1-p) == 70*p # p == 30/100 == .3 # ================================= # CI for binomial parameter p using standard normal approximation normal.ci <- function(pHat, n, alpha){ zStar <- qnorm(p=(1 - alpha/2)); e <- zStar * sqrt(pHat * (1-pHat) / n) pHat + c(-e, e)} # ================================= # example 6.49 cancer normal.ci(pHat=.040, n=10000, alpha=.05) # example 6.51 cancer # exact method for CI for binomial parameter p p1 <- seq(0,1,by=.01) p2 <- p1 b1 <- 1 - pbinom(q=1, size=20, prob=p1) b2 <- pbinom(q=2, size=20, prob=p2) cbind(p1, b1, p2, b2) # choose p1 == 0.01 and p2 == 0.32 ci <- c(0.01, 0.32); ci # example 6.52 health promotion # test for normal approximation --- fails test n <- 10 pHat <- .6 qHat <- 1 - pHat n * pHat * qHat # exact method, 99% CI p1 <- seq(0,1,by=.01) p2 <- p1 b1 <- 1 - pbinom(q=5, size=10, prob=p1) b2 <- pbinom(q=6, size=10, prob=p2) cbind(p1, b1, p2, b2) # choose p1 == 0.19 and p2 == 0.92 ci <- c(0.19, 0.92); ci # using R's exact binomial test, binom.test binom.test(x=6, n=10, alternative="two.sided", conf.level=.99) # => # 99 percent confidence interval: # 0.1909163 0.9232318 # ================================= # example 6.57 exact CI for Poisson parameter # load package epitools library(epitools) pois.exact(x=12, pt=120000, conf.level=.95) # => # x pt rate lower upper conf.level # 1 12 120000 1e-04 5.16717e-05 0.0001746797 0.95 # values for lambda: # 5.16 and 17.47 events per 10^5 person-years is 95% CI for lambda # corresponding values for mu == lambda * T: # 6.20 and 20.962 # these are exactly the values for x == 12 alpha == .95 in table 8, p.827 # there is also a poisson.test() in the stats package # ================================= # upper one-sided CI for binomial parameter p upper.onesided.normal.ci <- function(pHat, n, alpha){ pHat - qnorm(alpha) * sqrt(pHat * (1-pHat) / n) } # ================================= # example 6.57 exact 95% CI for Poisson parameter # check n <- 100 pHat <- .40 n * pHat * (1-pHat) alpha <- .95 upper.onesided.normal.ci(pHat=.40, n=100, alpha=.95) # ================================= # lower one-sided CI for binomial parameter p lower.onesided.normal.ci <- function(pHat, n, alpha){ pHat + qnorm(alpha) * sqrt(pHat * (1-pHat) / n) } # ================================= # example 6.57 exact 95% CI for Poisson parameter lower.onesided.normal.ci(pHat=.6, n=100, alpha=.95)