To follow along, use the lab06 code.

library("tidyverse")

Set your working directory using Session > Set Working Directory > Choose Directory in RStudio or using the setwd() function. You can also save creativity.csv into your working directory if you want.

Normal data

Background

Suppose \(Y_i\stackrel{ind}{\sim} N(\mu,\sigma^2)\). Recall that a Bayesian analysis with the default prior \(p(\mu,\sigma^2) \propto 1/\sigma^2\) provides the same analysis for \(\mu\) as a non-Bayesian analysis. That is

  • MLE and the posterior mean/median for \(\mu\) are the same
  • confidence and credible intervals for \(\mu\) are exactly the same

This is because the posterior for \(\mu\) is \[ \mu|y \sim t_{n-1}(\overline{y}, s^2/n) \] which means that \[ \left.\frac{\mu-\overline{y}}{s/\sqrt{n}} \right|y \sim t_{n-1}(0,1). \] while the sampling distribution for \(\overline{y}\) is such that \[ T=\frac{\overline{Y}-\mu}{S/\sqrt{n}} \sim t_{n-1}(0,1). \]

Please note the difference in these two statements is what is considered random. In the first two \(\mu\) is considered random while the data which are used to calculate \(\overline{y}\) and \(s\) are fixed. This does not mean that a Bayesian considered the actual true value of \(\mu\) to be random. Instead it means that we are expressing our uncertainty in \(\mu\) through probability. In the last statement, the data are considered random which is why \(\overline{Y}\) and \(S\) are capitalized while \(\mu\) is considered fixed.

Manual analysis

Suppose you observe the following data

set.seed(20180219)
(y <- rnorm(10, mean = 3.2, sd = 1.1)) # 3.2 and 1.1 are the population parameters
##  [1] 1.974293 4.759967 2.812764 0.758789 2.256666 4.880966 3.732871 3.929907 4.613458
## [10] 1.649653

Then you can manually construct an MLE and posterior mean using mean().

(ybar <- mean(y))
## [1] 3.136933

and a 95% credible/confidence interval using

n <- length(y)
s <- sd(y)

ybar + c(-1,1)*qt(.975, df = n-1)*s/sqrt(n)
## [1] 2.099208 4.174659

Built-in analysis

You can use the t.test() function to perform this for you.

t.test(y)
## 
##  One Sample t-test
## 
## data:  y
## t = 6.8383, df = 9, p-value = 7.572e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  2.099208 4.174659
## sample estimates:
## mean of x 
##  3.136933

You can extract the estimates and confidence/credible intervals, but first you need to know how to access the appropriate objects within the t-test object.

tt <- t.test(y)
names(tt)
##  [1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"    "null.value" 
##  [7] "stderr"      "alternative" "method"      "data.name"
str(tt)
## List of 10
##  $ statistic  : Named num 6.84
##   ..- attr(*, "names")= chr "t"
##  $ parameter  : Named num 9
##   ..- attr(*, "names")= chr "df"
##  $ p.value    : num 7.57e-05
##  $ conf.int   : num [1:2] 2.1 4.17
##   ..- attr(*, "conf.level")= num 0.95
##  $ estimate   : Named num 3.14
##   ..- attr(*, "names")= chr "mean of x"
##  $ null.value : Named num 0
##   ..- attr(*, "names")= chr "mean"
##  $ stderr     : num 0.459
##  $ alternative: chr "two.sided"
##  $ method     : chr "One Sample t-test"
##  $ data.name  : chr "y"
##  - attr(*, "class")= chr "htest"

It isn’t always obvious what object we want and thus we often need to just look at all of them until we figure out which ones we need. In this case, conf.int seems like a good guess for the confidence/credible intervals. It turns out that estimate will gives us the MLE and posterior mean/median.

tt$estimate
## mean of x 
##  3.136933
tt$conf.int
## [1] 2.099208 4.174659
## attr(,"conf.level")
## [1] 0.95

Modifying arguments

We can also change the default argument values to

  • change the type of hypothesis test
  • change the null value
  • change the confidence level

Suppose we wanted to test \(H_0:\mu\ge 1\) vs \(H_a:\mu < 1\) at a signficance level of 0.9 and/or we wanted to construct a 90% one-sided lower confidence interval.

t.test(y, alternative = "less",    mu = 1,  conf.level = 0.9)
## 
##  One Sample t-test
## 
## data:  y
## t = 4.6583, df = 9, p-value = 0.9994
## alternative hypothesis: true mean is less than 1
## 90 percent confidence interval:
##      -Inf 3.771374
## sample estimates:
## mean of x 
##  3.136933

Suppose we wanted to test \(H_0:\mu\le -1\) vs \(H_a:\mu > -1\) at a signficance level of 0.99 and/or we wanted to construct a 99% one-sided upper confidence interval.

t.test(y, alternative = "greater", mu = -1, conf.level = 0.99)
## 
##  One Sample t-test
## 
## data:  y
## t = 9.0182, df = 9, p-value = 4.199e-06
## alternative hypothesis: true mean is greater than -1
## 99 percent confidence interval:
##  1.842647      Inf
## sample estimates:
## mean of x 
##  3.136933

Activity

Using the following data, compare the point estimate and confidence/credible intervals obtained using the t.test() function to estimates and intervals you create yourself.

If you have more time, try to play around with the function arguments to create one-sided confidence intervals, change the null hypothesis value, and change the confidence/significance level.

set.seed(1)
y <- rnorm(1001, mean = 256, sd = 34.6)
Click for solution 
t.test(y)
## 
##  One Sample t-test
## 
## data:  y
## t = 225.84, df = 1000, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  253.4154 257.8578
## sample estimates:
## mean of x 
##  255.6366
mean(y)
## [1] 255.6366
mean(y) + c(-1,1)*qt(.975, df = length(y)-1)*sd(y)/sqrt(length(y))
## [1] 253.4154 257.8578

Binomial data

Bayesian analysis

Let \(Y\sim Bin(n,\theta)\). Recall that our default prior for \(\theta\) is \(\theta \sim Be(1,1)\), which is equivalent to a \(Unif(0,1)\) prior. The posterior under this prior is \(\theta|y \sim Be(1+y,1+n-y)\). In order to perform this analysis, we simply use the beta distribution in R.

Suppose you observe 9 successes out of 13 attempts.

a <- b <- 1
n <- 13
y <- 9

The posterior is

ggplot(data.frame(x=c(0,1)), aes(x=x)) + 
  stat_function(fun = dbeta, 
                args = list(shape1 = a+y,
                            shape2 = b+n-y)) + 
  labs(x = expression(theta),
       y = paste(expression("p(",theta,"|y)")),
       title = "Posterior distribution for probability of success") +
  theme_bw()

The posterior expectation is

(a+y)/(a+b+n)
## [1] 0.6666667

A 95% equal-tail credible interval is

qbeta(c(.025,.975), a+y, b+n-y)
## [1] 0.4189647 0.8724016

The probability that \(\theta\) is greater than 0.5, i.e. \(P(\theta>0.5|y)\) is

1-pbeta(0.5, a+y, b+n-y)
## [1] 0.9102173

Jeffreys prior activity

An alternative prior is called Jeffreys prior and it corresponds to a Be(0.5,0.5) prior. Suppose you observed 17 successes out of 20 attempts and you are willing to assume independence and a common probability of success. Use Jeffreys prior on this probability of success to do the following

  • Plot the posterior for the true probability of success
  • Calculate the posterior median for the true probability of success
  • Calculate a one-sided upper 95% credible interval for the true probability of success
  • Calculate the probability that the true probability of success is greater than 0.9.
Click for solution 
a <- b <- 0.5
n <- 20
y <- 17
curve(dbeta(x, a+y, b+n-y))

qbeta(0.5, a+y, b+n-y)
## [1] 0.8439957
c(qbeta(0.05, a+y, b+n-y), 1)
## [1] 0.6863196 1.0000000
1-pbeta(0.9, a+y, b+n-y)
## [1] 0.2137305

Non-Bayesian analyses when n is small

To perform non-Bayesian analyses when n is small, use the binom.test function. This will calculate pvalues and confidence intervals (based on inverting hypothesis tests).

Suppose you observe 9 successes out of 13 attempts.

n <- 13
y <- 9
binom.test(y, n)
## 
##  Exact binomial test
## 
## data:  y and n
## number of successes = 9, number of trials = 13, p-value = 0.2668
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.3857383 0.9090796
## sample estimates:
## probability of success 
##              0.6923077

This is equivalent to calculating the probability of observing y equal to 0, 1, 2, 3, 4, 9, 10, 11, 12, 13 if \(\theta\) is 0.5.

sum(dbinom(c(0:4,9:13), size = n, prob = 0.5))
## [1] 0.2668457

Recall that there is a one-to-one correspondence between pvalues and confidence intervals. This confidence interval is constructed by finding those values for \(\theta\) such that the pvalue is half of the confidence level (since it is a two-sided interval).

ci <- binom.test(y,n)$conf.int
binom.test(y, n, p=ci[2])$p.value # This one matches exactly
## [1] 0.025
binom.test(y, n, p=ci[1])$p.value # This one is close
## [1] 0.04124282

binom.test Activity

Suppose you observe 11 successes out of 12 attempts. Calculate a pvalue for the two-sided test that \(\theta=0.5\) and construct a 95% confidence interval.

Click for solution 
binom.test(11,12)
## 
##  Exact binomial test
## 
## data:  11 and 12
## number of successes = 11, number of trials = 12, p-value = 0.006348
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.6152038 0.9978924
## sample estimates:
## probability of success 
##              0.9166667

Non-Bayesian analyses when n is large (and p is not close to 0 or 1)

If you observe 78 successes out of 100 attempts, then you can use the prop.test() function to generate a number of statistics automatically based on the CLT.

n <- 100
y <- 78
prop.test(y,n)
## 
##  1-sample proportions test with continuity correction
## 
## data:  y out of n, null probability 0.5
## X-squared = 30.25, df = 1, p-value = 3.798e-08
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.6838686 0.8541734
## sample estimates:
##    p 
## 0.78

The estimate is

pt <- prop.test(y,n)
pt$estimate
##    p 
## 0.78

An approximate 95% confidence interval is

pt$conf.int
## [1] 0.6838686 0.8541734
## attr(,"conf.level")
## [1] 0.95

We can construct this ourself using the following formula

p <- y/n
p + c(-1,1)*qnorm(.975)*sqrt(p*(1-p)/n)
## [1] 0.6988092 0.8611908

The results don’t quite match because prop.test uses the continuity correction (which is the appropriate thing to do).

prop.test(y, n, correct = FALSE)
## 
##  1-sample proportions test without continuity correction
## 
## data:  y out of n, null probability 0.5
## X-squared = 31.36, df = 1, p-value = 2.144e-08
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.6892965 0.8499872
## sample estimates:
##    p 
## 0.78

Unfortunately, the results still don’t quite match. It turns out there is a large literature on how to do this and the suggestions are to either 1) use Wilson’s interval (which is what prop.test does) or to use 2) Jeffreys, i.e. a Beta(1/2,1/2) prior in a Bayesian analysis.

Activity

Suppose you observed 694 success out of 934 attempts. Compute an approximate 95% equal-tail confidence interval using prop.test and compare this to a 95% equal-tail confidence interval you construct using the Central Limit Theorem.

Click for solution 
y <- 694
n <- 934
prop.test(y, n)$conf.int
## [1] 0.7135099 0.7705425
## attr(,"conf.level")
## [1] 0.95
p <- y/n
p + c(-1,1)*qnorm(.975)*sqrt(p*(1-p)/n)
## [1] 0.7150178 0.7710636
# If you turn off the continuity correction, you will get closer
prop.test(y, n, correct = FALSE)$conf.int
## [1] 0.7140620 0.7700283
## attr(,"conf.level")
## [1] 0.95
If you really want to figure out what the function is doing, you can look at the function by just typing prop.test and hitting enter.

Activity

Suppose you observe 0 success out of 77 attempts. Compare 95% confidence intervals given by prop.test and binom.test to an interval you construct based on the CLT and to 95% credible intervals.

Click for solution

So this is a bit of a trick question since there were 0 successes. When you run prop.test and binom.test you are given one-sided confidence intervals. The CLT interval doesn’t exist since the standard error is zero. The appropriate credible interval to use is a one-sided interval.

y <- 0
n <- 77
prop.test(y, n)$conf.int
## [1] 0.00000000 0.05921002
## attr(,"conf.level")
## [1] 0.95
binom.test(y, n)$conf.int
## [1] 0.00000000 0.04677807
## attr(,"conf.level")
## [1] 0.95
qbeta(c(0,.95), 1+y, 1+n-y)
## [1] 0.00000000 0.03767863

Reading data from files

First, let’s write some data to a file.

set.seed(20180220)
# Generate some simulated data
n <- 100
d <- data.frame(rep = 1:n,
                response = sample(c("Yes","No"), n, replace=TRUE, prob = c(.2,.8)),
                measurement = rnorm(n, mean = 55, sd = 12))

# Write it to a file
# make sure you have set your working directory to someplace where you want this
# file to be written
write.csv(d, 
          file = "data.csv",
          row.names = FALSE)

Alternatively, you could use the write_csv() function in the readr package.

install.packages("readr")
library("readr") 
write_csv(d, path = "data.csv")

Now let’s read this data back in.

my_data <- read.csv("data.csv")

If you want to delete the file, you can run the following

if (file.exists("data.csv")) file.remove("data.csv")
## [1] TRUE

Take a look at the data to make sure it looks correct:

head(my_data)
##   rep response measurement
## 1   1       No    49.16525
## 2   2       No    66.62369
## 3   3       No    35.14872
## 4   4       No    61.15329
## 5   5       No    61.34038
## 6   6       No    53.79730
str(my_data)
## 'data.frame':    100 obs. of  3 variables:
##  $ rep        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ response   : chr  "No" "No" "No" "No" ...
##  $ measurement: num  49.2 66.6 35.1 61.2 61.3 ...

Binomial data

To use prop.test() and binom.test(), you need to calculate the number of successes and the number of attempts.

y <- sum(my_data$response == "Yes")
n <- length(my_data$response)
prop.test(y, n)
## 
##  1-sample proportions test with continuity correction
## 
## data:  y out of n, null probability 0.5
## X-squared = 32.49, df = 1, p-value = 1.198e-08
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.1375032 0.3052597
## sample estimates:
##    p 
## 0.21
binom.test(y, n)
## 
##  Exact binomial test
## 
## data:  y and n
## number of successes = 21, number of trials = 100, p-value = 4.337e-09
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.1349437 0.3029154
## sample estimates:
## probability of success 
##                   0.21

Normal data

To analyze the normal data, you can just use t.test() directly.

t.test(my_data$measurement)
## 
##  One Sample t-test
## 
## data:  my_data$measurement
## t = 42.956, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  50.91322 55.84453
## sample estimates:
## mean of x 
##  53.37888

Online activity

Read in the data at creativity.csv and then construct confidence/credible intervals for mean creativity score for both the Intrinsic and Extrinsic groups.

Click for solution

There are a variety of ways to do this. I will construct two new data frames to contain the Intrinsic and Extrinsic data and then construct the intervals.

creativity <- read.csv("http://www.jarad.me/courses/stat587Eng/labs/lab06/creativity.csv")

intrinsic_score <- creativity$Score[creativity$Treatment == "Intrinsic"]
extrinsic_score <- creativity$Score[creativity$Treatment == "Extrinsic"]
t.test(intrinsic_score)
## 
##  One Sample t-test
## 
## data:  intrinsic_score
## t = 21.941, df = 23, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  18.00869 21.75798
## sample estimates:
## mean of x 
##  19.88333
t.test(extrinsic_score)
## 
##  One Sample t-test
## 
## data:  extrinsic_score
## t = 14.37, df = 22, p-value = 1.16e-12
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  13.46774 18.01052
## sample estimates:
## mean of x 
##  15.73913
# Using the `subset` command which subsets the data.frame
intrinsic <- subset(creativity, Treatment == "Intrinsic")
extrinsic <- subset(creativity, Treatment == "Extrinsic")

t.test(intrinsic$Score)
## 
##  One Sample t-test
## 
## data:  intrinsic$Score
## t = 21.941, df = 23, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  18.00869 21.75798
## sample estimates:
## mean of x 
##  19.88333
t.test(extrinsic$Score)
## 
##  One Sample t-test
## 
## data:  extrinsic$Score
## t = 14.37, df = 22, p-value = 1.16e-12
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  13.46774 18.01052
## sample estimates:
## mean of x 
##  15.73913
# Another way to subset uses the dplyr package
# This will be prettier when you have a long sequence of commands
library(dplyr)
intrinsic <- creativity %>% filter(Treatment == "Intrinsic")
extrinsic <- creativity %>% filter(Treatment == "Extrinsic")

If you want to find out more about these data, take a look at the help file for case0101 in the Sleuth3 package.

install.packages("Sleuth3")
library("Sleuth3")
?case0101