# Chihara and Hesterberg, chap 7 # Chapter 7 Classical Inference: Confidence Intervals # 7.1 confidence intervals for means # 7.1.1 CI for normal with known sigma # =================================== # Example 7.1 weights # Xi ~ N(mu, sd^2) # XBar ~ N(mu, sd^2 / n) # Z == (XBar - mu) / (sd / sqrt(n)) ~ N(0,1) # .05 == P(z.025 <= Z <= z.975) # == P(z.025 <= (XBar - mu) / (sd / sqrt(n)) <= z.975) # == P(z.025 * SE <= (XBar - mu) <= z.975 * SE) and SE == sd / sqrt(n) # == P(xBar - z.975 * SE <= mu <= xBar - z.025 * SE) and SE == sd / sqrt(n) mean <- 101 sd <- 24.6 xBar <- 95 n <- 150 alpha <- .05 z.975 <- qnorm(1 - alpha/2) # => 1.959964 se <- sd / sqrt(n) (me <- z.975 * se) # margin of error # => 3.936748 (ci <- xBar + c(-me, me)) # => 91.06325 98.93675 # =================================== # simulation : draw random samples of size 30 from N(25, 4^2) # plot the CIs ... how many intervals contain mu? #set.seed(1) counter <- 0 # set counter to 0 plot(x = c(22, 28), y = c(1, 100), type = "n", xlab = "", ylab = "", main="Confidence intervals") abline(v = 25, col="red") # vertical line at mu for (i in 1:1000) { x <- rnorm(30, 25, 4) # draw a random sample of size 30 L <- mean(x) - 1.96*4/sqrt(30) # lower limit U <- mean(x) + 1.96*4/sqrt(30) # upper limit if (L < 25 && 25 < U) # check to see if 25 is in interval counter <- counter + 1 # increase counter by 1 if (i <= 100) #plot first 100 intervals segments(L, i, U, i, col="sandybrown") } abline(v = 25, col = "red") #vertical line at mu counter/1000 # proportion of times interval contains mu. # => 0.946 # =================================== # Example 7.3 90% CI for the mean # X ~ N(mu=??, sigma=2.5) sample <- c(3.4, 2.9, 2.8, 5.1, 6.3, 3.9) sigma <- 2.5 (n <- length(sample)) (xBar <- mean(sample)) # => 4.066667 alpha <- .10 (z.95 <- qnorm(1 - alpha/2)) # => 1.644854 se <- sigma / sqrt(n) (me <- z.95 * se) # margin of error # => 1.678772 (ci <- xBar + c(-me, me)) # => 2.387895 5.745438 # =================================== # Example 7.4 90% sample size # X ~ N(mu=??, sigma=24.6) sigma <- 24.6 alpha <- .05 (z.975 <- qnorm(1 - alpha/2)) # me == z.975 * se == z.975 * sigma / sqrt(n) <= 5 # => sqrt(n) >= z.975 * sigma / 5 (n <- ceiling((z.975 * sigma / 5)^2)) # => 93 # =================================== # Section 7.1.2 # simulate distribution of t statistic N <- 10^4 w <- numeric(N) n <- 15 # sample size for (i in 1:N) { x <- rnorm(n, 25, 7) #draw a size 15 sample from N(25, 7^2) xbar <- mean(x) s <- sd(x) w[i] <- (xbar-25) / (s/sqrt(n)) } hist(w, col="wheat3") dev.new() qqnorm(w, pch = ".") abline(0, 1, col = 2) # y = x line # pch = "." is point character. This option says to use . for the points. # =================================== # pdfs for four distributions curve(dnorm, from=-3, to=3, col="dark red", xlab="", ylab="probability", lty=1, main="PDFs for four distributions") curve(dt(x, df=4), lty=2, col="plum", add=TRUE) curve(dt(x, df=8), lty=3, col="plum2", add=TRUE) curve(dt(x, df=16), lty=4, col="plum4", add=TRUE) legend("topright", lty=1:4, legend=c("N(0,1)", "df=4", "df=8", "df=16"), col=c("dark red", "plum", "plum2", "plum4"), lwd=2) # =================================== # Example 7.5 90% CI # X ~ N(mu=??, sigma=??) xBar <- 110 sd <- 7.5 n <- 28 alpha <- .10 (z.95 <- qt(1 - alpha/2, df=n - 1)) # => 1.703288 se <- sd / sqrt(n) (me <- z.95 * se) # margin of error # => 2.414184 (ci <- xBar + c(-me, me)) # => 107.5858 112.4142 # =================================== # Example 7.6 99% CI # X ~ N(mu=??, sigma=??) xBar <- 3398.317 sd <- 485.691 n <- 521 alpha <- .01 (z.995 <- qt(1 - alpha/2, df=n - 1)) # => 2.585317 se <- sd / sqrt(n) (me <- z.995 * se) # margin of error # => 55.01169 (ci <- xBar + c(-me, me)) # => 3343.305 3453.329 # qqplot NCBirths2004 <- read.csv("NCBirths2004.csv") girls <- subset(NCBirths2004, select="Weight", subset=Gender == "Female", drop=TRUE) qqnorm(girls, col="mediumorchid", main="Weights of baby girls") qqline(girls, col = "rosybrown1") t.test(girls, conf.level = 1 - alpha)$conf # => [1] 3343.305 3453.328 # => attr(,"conf.level") # => [1] 0.99 # =================================== # Example 7.7 Simulation 95% confidence interval from # skewed gamma distribution # set.seed(0) curve(dgamma(x, shape=5, rate=2), from=0, to=6, col="dark red", xlab="", ylab="probability", lty=1, main="Skewed gamma distribution") tooLow <- 0 #set counter to 0 tooHigh <- 0 #sest counter to 0 n <- 20 # sample size N <- 10^5 for (i in 1:N) { x <- rgamma(n, shape=5, rate=2) xbar <- mean(x) s <- sd(x) lower <- xbar - abs(qt(.025, n-1))*s/sqrt(n) upper <- xbar + abs(qt(.025, n-1))*s/sqrt(n) if (upper < 5/2) tooLow <- tooLow + 1 if (lower > 5/2) tooHigh <- tooHigh + 1 } tooLow/N# => 0.0433 tooHigh/N # => 0.01359 # =================================== # qqplots # uniform, n == 10 n <- 10 my.Us<-numeric(10000) my.Ts<-numeric(10000) #set.seed(300) for (i in 1:10000) { x<-runif(n) my.Us[i]<- (mean(x) - .5) / (sd(x) / sqrt(n)) y<-rt(n, df=n-1) my.Ts[i]<- (mean(y) - 0) / (sd(y) / sqrt(n)) } qqplot(my.Ts, my.Us, pch=".", xlim=c(-3, 3), ylim=c(-3,3), col="yellowgreen", main="Uniform, n=10") qqline(my.Ts, col="dark red") # exponential, n == 10 n <- 10 my.Es<-numeric(10000) my.Ts<-numeric(10000) #set.seed(300) for (i in 1:10000) { x<-rexp(n) my.Es[i]<- (mean(x) - 1) / (sd(x) / sqrt(n)) y<-rt(n, df=n-1) my.Ts[i]<- (mean(y) - 0) / (sd(y) / sqrt(n)) } qqplot(my.Ts, my.Es, pch=".", xlim=c(-3, 3), ylim=c(-3,3), col="maroon2", main="Exponential, n=10") qqline(my.Ts, col="dark red") # exponential, n == 100 n <- 100 my.Es<-numeric(10000) my.Ts<-numeric(10000) #set.seed(300) for (i in 1:10000) { x<-rexp(n) my.Es[i]<- (mean(x) - 1) / (sd(x) / sqrt(n)) y<-rt(n, df=n-1) my.Ts[i]<- (mean(y) - 0) / (sd(y) / sqrt(n)) } qqplot(my.Ts, my.Es, pch=".", xlim=c(-3, 3), ylim=c(-3,3), col="tan4", main="Exponential, n=100") qqline(my.Ts, col="dark red") # =================================== # 7.1.3 confidence intervals for a difference in means # Example 7.8 reading scores xBar <- 51.48 s1 <- 11.01 n1 <- 21 yBar <- 41.52 s2 <- 17.14 n2 <- 23 se <- sqrt(s1^2 / n1 + s2^2 / n2) nu <- (s1^2 / n1 + s2^2 / n2)^2 / ((s1^2 / n1)^2 / (n1 - 1) + (s2^2 / n2)^2 / (n2 - 1)) alpha <- .05 q.975 <- qt(1 - alpha/2, df = floor(nu)) # 95% CI me <- q.975 * se (ci <- (xBar - yBar) + c(-me, me)) # => 1.234327 18.685673 # t.test # reading data is available at DASL # http://lib.stat.cmu.edu/DASL/DataArchive.html # http://lib.stat.cmu.edu/DASL/Datafiles/DRPScores.html # copy that data into the file "Reading.csv" (and add the commas to separate fields) Reading <- read.csv("Reading.csv") str(Reading) summary(Reading) treated <- subset(Reading, select="Response", subset= Treatment == "Treated", drop=TRUE) control <- subset(Reading, select="Response", subset= Treatment == "Control", drop=TRUE) t.test(treated, control)$conf # => [1] 1.23302 18.67588 # => attr(,"conf.level") # => [1] 0.95 plot.ecdf(treated, xlab="x", ylab="Fn(x)", xlim=c(0, 90), col="olivedrab3", main="Reading Scores") plot.ecdf(control, col="palegreen2", add=TRUE) legend(5, .9, legend=c("Control", "Treated"), lwd=2, col=c("palegreen2", "olivedrab3")) # =================================== # qqplots # use df == nu == n1 - 1 for the theoretical t distribution # exponential, n1 == 10, n2 == 10 n1 <- 10 n2 <- 10 my.Ts<-numeric(10000) theoretical.Ts<-numeric(10000) #set.seed(300) for (i in 1:10000) { x <- rexp(n1) y <- rexp(n2) my.Ts[i]<- (mean(x) - mean(y)) / sqrt(sd(x)^2 / n1 + sd(y)^2 / n2) nu <- n1 - 1 y<-rt(n1, df=nu) theoretical.Ts[i]<- (mean(y) - 0) / (sd(y) / sqrt(n1)) } qqplot(theoretical.Ts, my.Ts, pch=".", xlim=c(-3, 3), ylim=c(-3,3), col="yellowgreen", main="Exponential: n1=10, n2 =10") qqline(theoretical.Ts, col="dark red") # exponential, n1 == 10, n2 == 20 n1 <- 10 n2 <- 20 my.Ts<-numeric(10000) theoretical.Ts<-numeric(10000) #set.seed(300) for (i in 1:10000) { x <- rexp(n1) y <- rexp(n2) my.Ts[i]<- (mean(x) - mean(y)) / sqrt(sd(x)^2 / n1 + sd(y)^2 / n2) nu <- n1 - 1 y<-rt(n1, df=nu) theoretical.Ts[i]<- (mean(y) - 0) / (sd(y) / sqrt(n1)) } qqplot(theoretical.Ts, my.Ts, pch=".", xlim=c(-3, 3), ylim=c(-3,3), col="maroon2", main="Exponential: n1=10, n2 =20") qqline(theoretical.Ts, col="dark red") # exponential, n == 100 n1 <- 10 n2 <- 100 my.Ts<-numeric(10000) theoretical.Ts<-numeric(10000) #set.seed(300) for (i in 1:10000) { x <- rexp(n1) y <- rexp(n2) my.Ts[i]<- (mean(x) - mean(y)) / sqrt(sd(x)^2 / n1 + sd(y)^2 / n2) nu <- n1 - 1 y<-rt(n1, df=nu) theoretical.Ts[i]<- (mean(y) - 0) / (sd(y) / sqrt(n1)) } qqplot(theoretical.Ts, my.Ts, pch=".", xlim=c(-3, 3), ylim=c(-3,3), col="tan4", main="Exponential: n1=10, n2=100") qqline(theoretical.Ts, col="dark red") # =================================== # Example 7.9 weights of boys and girls TXBirths2004 <- read.csv("TXBirths2004.csv") head(TXBirths2004) TXBirths2004.boys <- subset(TXBirths2004, select=Weight, subset=Gender=="Male", drop=T) TXBirths2004.girls <- subset(TXBirths2004, select=Weight, subset=Gender=="Female", drop=T) t.test(TXBirths2004.boys, TXBirths2004.girls) # => Welch Two Sample t-test # => data: TXBirths2004.boys and TXBirths2004.girls # => t = 4.193, df = 1552.906, p-value = 2.908e-05 # => alternative hypothesis: true difference in means is not equal to 0 # => 95 percent confidence interval: # => 61.68415 170.12395 # => sample estimates: # => mean of x mean of y # => 3336.843 3220.939 par(mfrow=c(2,2)) # set layout hist(TXBirths2004.boys, xlab="Grams", ylab="Frequency", col="royalblue2", main = "Weights of Boys") qqnorm(TXBirths2004.boys, col="blue", main = "Boys") qqline(TXBirths2004.boys, col="dark red") hist(TXBirths2004.girls, xlab="Grams", ylab="Frequency", col="lightcyan4", main = "Weights of Girls") qqnorm(TXBirths2004.girls, col="lightcyan4", main = "Girls") qqline(TXBirths2004.girls, col="dark red") par(mfrow=c(1,1)) # reset layout # =================================== # 7.2 confidence intervals in general # =================================== # Example 7.10 confidence intervals in general alpha <- .05 z.975 <- qnorm(1 - alpha/2) # L <- xBar - z.975 * sigma / sqrt(n) # U <- xBar + z.975 * sigma / sqrt(n) # note sqrt missing on p.183, but present in equation 7.6, p.168 # =================================== # Example 7.11 confidence intervals in general # Xi ~ exp(lambda == ??) # sum Xi ~ gamma(n, lambda) # lambda sum Xi ~ gamma(n, 1) n <- 5 alpha <- .05 q.025 <- qgamma(alpha / 2, n, 1) q.975 <- qgamma(1 - alpha / 2, n, 1) # 1 - alpha == P(q.025 <= lambda sum Xi <= q.975) # == P(q.025 / sum Xi <= lambda <= q.975 / sum Xi) # 95% CI # (q.025 / sum Xi, q.975 / sum Xi) # example data <- c(2, 2.5, 3, 4, 9) n <- length(data) sum <- sum(data) alpha <- .05 q.025 <- qgamma(alpha / 2, n, 1) q.975 <- qgamma(1 - alpha / 2, n, 1) (ci <- c(q.025 / sum, q.975 / sum)) # => 0.07919446 0.49958969 xBar <- mean(data) (lambdaHat <- 1 / xBar) # => 0.2439024 (a <- -q.025 / sum + lambdaHat) # => 0.164708 (b <- q.975 / sum - lambdaHat) # => 0.2556873 # 95% CI c(lambdaHat - a, lambdaHat + b) # gamma(5, 1) curve(dgamma(x, 5, 1), from=0, to=12, xlab="x", ylab="probability density", col="dark red", main="Gamma(5, 1)") poly.xs <- seq(q.025, q.975, by=.01) poly.ys <- dgamma(poly.xs, 5, 1) polygon(c(q.025, poly.xs, q.975), c(0, poly.ys, 0), col="orange") text(.5, .04, labels=expression(0.025)) text(11.2, .025, labels=expression(0.025)) # =================================== # Example 7.12 95% confidence interval for lambda # X ~ exp(lambda=??) # Y == lambda X has PDF given by e^(-y) for y > 0 # F_Y(y) == 1 - exp(-y) # solve F_Y(y) == .025 and F_Y(y) == .975 # 1 - exp(-y) == .025 => .975 == exp(-y) => y == -ln(.975) q.025 <- -log(.975) # => 0.02531781 # 1 - exp(-y) == .975 => .025 == exp(-y) => y == -ln(.025) q.975 <- -log(.025) # => 3.688879 # then, 0.95 == P(q.025 < Y < q.975) == P(q.025 < lambda X < q.975) # == P(q.025 / X < lambda < q.975 / X) # example x <- 0.03 (ci <- c(q.025 / x, q.975 / x)) # => 0.8439269 122.9626485 # =================================== # Example 7.13 tanks q1 <- .025^(1/400) q2 <- .975^(1/400) x <- 10000 (ci <- c(x / q2, x / q1)) # => 10000.63 10092.65 # =================================== # 7.2.1 location and scale parameters curve(dgamma(x + 2, 4, 1.3), from=-2, to=10, xlab="x", ylab="probability density", col="dark red", main="Sampling distributions for a location parameter") curve(dgamma(x, 4, 1.3), lty=2, col="magenta3", add=TRUE) curve(dgamma(x - 2, 4, 1.3), lty=3, col="mediumseagreen", add=TRUE) legend("topright", title="Sampling distribution when", lty=c(1,2,3), lwd=1, col=c("dark red", "magenta3", "mediumseagreen"), legend=c(expression(paste(theta," = 0")), expression(paste(theta," = 2")), expression(paste(theta," = 4")))) curve(dgamma(x, 3, 2/2), from=0, to=12, xlab="x", ylab="probability density", col="dark red", main="Sampling distributions for a scale parameter") curve(dgamma(x, 3, 2/4), lty=2, col="magenta3", add=TRUE) curve(dgamma(x, 3, 2/8), lty=3, col="mediumseagreen", add=TRUE) legend("topright", title="Sampling distribution when", lty=c(1,2,3), lwd=1, col=c("dark red", "magenta3", "mediumseagreen"), legend=c(expression(paste(theta," = 2")), expression(paste(theta," = 4")), expression(paste(theta," = 8")))) # =================================== # 7.3 one-sided confidence intervals # =================================== # Example 7.14 lead levels xBar <- 7 s <- 2 n <- 15 alpha <- .05 q.95 <- qt(1 - alpha, df=n-1) L <- xBar - q.95 * s / sqrt(n) # => 6.090463 # 95% CI for mean # ci <- c(L, \infty) # =================================== # one-sided confidence intervals using t.test NCBirths2004 <- read.csv("NCBirths2004.csv") t.test(NCBirths2004$Weight, alt="greater")$conf # => [1] 3422.98 Inf # => attr(,"conf.level") # => [1] 0.95 # =================================== # Example 7.15 95% lower confidence interval for lambda # from Example 7.12 # X ~ exp(lambda=??) # Y == lambda X has PDF given by e^(-y) for y > 0 # F_Y(y) == 1 - exp(-y) # solve F_Y(y) == .95 # 1 - exp(-y) == .95 => .05 == exp(-y) => y == -ln(.05) q.95 <- -log(.05) # => 2.995732 # then, 0.95 == P(Y < q.95) == P(lambda X < q.95) # == P(lambda < q.95 / X) # 95% CI for lambda # U <- q.95 / X # ci <- c(-\infty, U) # =================================== # 7.4 confidence intervals for proportions # =================================== # Example 7.17 confidence intervals for proportion prop.test(899, 1308, conf.level=.9, correct=FALSE)$conf # => [1] 0.6658563 0.7079882 # => attr(,"conf.level") # => [1] 0.9 prop.test(899, 1308, conf.level=.9, correct=FALSE, alt="greater")$conf # => [1] 0.6706553 1.0000000 # => attr(,"conf.level") # => [1] 0.9 # =================================== # Example 7.18 Agresti-Coull 95% confidence interval for proportion x <- 130 n <- 210 x.tilda <- x + 2 n.tilda <- n + 4 p.tilda <- x.tilda / n.tilda (ci <- p.tilda + 1.96 * sqrt(p.tilda * (1 - p.tilda) / n.tilda) * c(-1, 1)) # => 0.5516852 0.6819597 # =================================== # Example 7.19 sample size for confidence interval for proportion # me <- 1.96 * sqrt(p.tilda * (1 - p.tilda) / n.tilda) # eqn = me <= .04 # solve for n # use p.tilda == .5 to maximize the numerator # 1.96 * sqrt(.5 * (1 - .5) / n.tilda) == .04 # sqrt(n.tilda) == 1.96 * sqrt(.5 * (1 - .5) / .04^2) (n.tilda <- (1.96 * sqrt(.5 * (1 - .5) / .04^2))^2) (n <- n.tilda - 4) # => 596.25 # =================================== # Example 7.20 sample size for confidence interval for proportion prop.test(c(475, 424), c(631, 677), correct=FALSE)$conf # => [1] 0.07687197 0.17608985 # => attr(,"conf.level") # => [1] 0.95 # =================================== # 7.5 bootstrap confidence intervals # =================================== # Example 7.21 One sample bootstrap t confidence interval Bangladesh <- read.csv("Bangladesh.csv") Arsenic <- subset(Bangladesh, select = Arsenic, drop = T) xbar <- mean(Arsenic) N <- 10^4 n <- length(Arsenic) Tstar <- numeric(N) #set.seed(100) for (i in 1:N) { x <-sample(Arsenic, size = n, replace = T) Tstar[i] <- (mean(x)-xbar)/(sd(x)/sqrt(n)) } quantile(Tstar, c(0.025, 0.975)) hist(Tstar, xlab = "T*", main = "Bootstrap distribution of T*", col="tan3") dev.new() qqnorm(Tstar, col="salmon1") qqline(Tstar, col="slategray4") #------------------------------------------------------- # Example 7.22 Verizon # 2-Sample bootstrap t confidence interval Verizon <- read.csv("Verizon.csv") Time.ILEC <- subset(Verizon, select=Time, Group == "ILEC", drop=T) Time.CLEC <- subset(Verizon, select=Time, Group == "CLEC", drop=T) thetahat <- mean(Time.ILEC)-mean(Time.CLEC) nx <- length(Time.ILEC) #nx=1664 ny <- length(Time.CLEC) #ny=23 SE <- sqrt(var(Time.ILEC)/nx + var(Time.CLEC)/ny) N <- 10000 Tstar <- numeric(N) set.seed(0) for(i in 1:N) { bootx <- sample(Time.ILEC, nx, replace=TRUE) booty <- sample(Time.CLEC, ny, replace=TRUE) Tstar[i] <- (mean(bootx) - mean(booty) - thetahat) / sqrt(var(bootx)/nx + var(booty)/ny) } thetahat - quantile(Tstar, c(.975, .025)) * SE t.test(Time.ILEC, Time.CLEC)$conf #---------------------------------------------------------------- # Exercise 18 %%Simulation to compare pooled/unpooled t-confidence intervals pooled.count <- 0 unpooled.count <- 0 m <- 80 n <- 80 B <- 10000 for (i in 1:B) { x <- rnorm(m, 8,10) y <- rnorm(n, 3, 15) CI.pooled <- t.test(x,y,var.equal=T)$conf CI.unpooled <- t.test(x,y)$conf if (CI.pooled[1] < 5 & 5 < CI.pooled[2]) pooled.count <- pooled.count + 1 if (CI.unpooled[1] < 5 & 5 < CI.unpooled[2]) unpooled.count <- unpooled.count + 1 } pooled.count/B unpooled.count/B