library("tidyverse")
## ── Attaching core tidyverse packages ─────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4.9000     ✔ readr     2.1.5     
## ✔ forcats   1.0.0          ✔ stringr   1.5.1     
## ✔ ggplot2   3.5.1          ✔ tibble    3.2.1     
## ✔ lubridate 1.9.3          ✔ tidyr     1.3.1     
## ✔ purrr     1.0.2          
## ── Conflicts ───────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Normal model

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 frequentist 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

# Normal 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
##  [9] 4.613458 1.649653

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

# Normal mean MLE and posterior expectation
(ybar <- mean(y))
## [1] 3.136933

and a 95% credible/confidence interval using

# Manual confidence/credible interval
n <- length(y)
s <- sd(y)
a <- 0.05

ybar + c(-1,1) * qt(1 - a/2, 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 (one-sample normal model)
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.

# T-test objects
tt <- t.test(y)
names(tt)
##  [1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"   
##  [6] "null.value"  "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.

# Normal estimators for the mean (mu)
tt$estimate # MLE and Posterior expectation
## mean of x 
##  3.136933
tt$conf.int # Credible and confidence interval
## [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.

# Normal arguments 
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.

# Normal activity data
set.seed(1)
y <- rnorm(1001, mean = 256, sd = 34.6)
Click for solution 
# Normal activity 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 model

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.

# Data
n <- 13
y <- 9

# Prior, Be(a,b)
a <- b <- 1

The posterior is

# Posterior
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

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

A 95% equal-tail credible interval is

# 95% equal-tail credible interval
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

# Binomial probabilities
1 - pbeta(0.5, a + y, b + n - y) # P(theta > 0.5 | 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, i.e. (L, 1), for the true probability of success \(\theta\)
  • Calculate the probability that the true probability of success \(\theta\) is greater than 0.9.
Click for solution 
# Binomial activity solutions
a <- b <- 0.5 # Jeffreys prior, Be(a,b)
n <- 20 # attempts
y <- 17 # successes

curve(dbeta(x, a + y, b + n - y))   # Posterior density

qbeta(0.5, a + y, b + n - y)        # Posterior median
## [1] 0.8439957
c(qbeta(0.05, a + y, b + n - y), 1) # One-sided credible interval
## [1] 0.6863196 1.0000000
1 - pbeta(0.9, a + y, b + n - y)    # P(theta > 0.9 | y)
## [1] 0.2137305

Frequentist analyses when n is small

To perform a frequentist 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.

# Small n data
n <- 13
y <- 9

(bt <- 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

The p-value here is 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. This range of integers from 0 to 4 and 9 to 13 is the as or more extreme regions.

# Binomial manual p value calculation
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 one minus the confidence level (since it is a two-sided interval).

# Binomial confidence interval
(ci <- bt$conf.int)
## [1] 0.3857383 0.9090796
## attr(,"conf.level")
## [1] 0.95
# Using the end points as the hypothesized probability should result in 
# a p-value that is half of 1 - confidence level
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
# Find the "correct" endpoint
# Create a function that should be zero when we have the correct probability
f <- function(p) {
  binom.test(y, n, p = p)$p.value - 0.025 # 0.025 is what the p-value should be
}

# Use unitroot to find value for p such that its p-value is 0.025
(u <- uniroot(f, lower = 0, upper = y/n))
## $root
## [1] 0.3772149
## 
## $f.root
## [1] -0.001481535
## 
## $iter
## [1] 19
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 0.000116585
# Check the p-value for this probability
binom.test(y, n, p = u$root)$p.value      
## [1] 0.02351846

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 
# Binomial binom.test activity solution
(bt <- 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
# Extract just the relevant results
bt$p.value
## [1] 0.006347656
bt$conf.int
## [1] 0.6152038 0.9978924
## attr(,"conf.level")
## [1] 0.95

Frequentist 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.

# Data with large n (and p not too close to 0 or 1)
n <- 100
y <- 78

When n is large, you can use the prop.test function.

# prop.test for large n (and p not too close to 0 or 1)
(pt <- 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

# Binomial MLE
pt$estimate
##    p 
## 0.78

An approximate 95% confidence interval is

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

We can construct this yourself using the following formula

# Manual CI
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).

# You should always use the continuity correct
# this is just for illustrative purposes
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 actually does) or to
2) perform a Bayesian analysis using Jeffreys prior, i.e. a Beta(1/2, 1/2).

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 
# Binomial activity data
y <- 694
n <- 934

# Use prop.test since n is large (and p is not too close to 0 or 1)
prop.test(y, n)$conf.int
## [1] 0.7135099 0.7705425
## attr(,"conf.level")
## [1] 0.95
# MLE
p <- y/n

# CI from CLT
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 to the 
# pure CLT interval
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.

# Binomial activity 3 data
y <- 0
n <- 77

# Large n (but p is close to 0)
prop.test(y, n)$conf.int
## [1] 0.00000000 0.05921002
## attr(,"conf.level")
## [1] 0.95
# Small n (but y/n is 0)
binom.test(y, n)$conf.int
## [1] 0.00000000 0.04677807
## attr(,"conf.level")
## [1] 0.95
# Bayesian one-sided interval with uniform prior
qbeta(c(0, .95), 1 + y, 1 + n - y)
## [1] 0.00000000 0.03767863

Reading data from files

Set working directory

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.

Write/read data from file

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/stat5870Eng/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