Before we run any analyses, we upload and inspect our dataset:
data = read.table("http://matt.colorado.edu/teaching/stats/labs/lab8-1.txt", header = TRUE) #load data
data #print out data
## Correct Group
## 1 587 0
## 2 342 1
## 3 326 0
## 4 362 1
## 5 602 1
## 6 648 0
## 7 494 0
## 8 647 1
## 9 336 0
## 10 579 0
## 11 696 1
## 12 346 1
## 13 696 1
## 14 350 0
## 15 681 0
## 16 352 1
## 17 527 1
## 18 672 0
## 19 353 1
## 20 685 1
## 21 699 0
## 22 688 0
## 23 691 0
## 24 344 1
## 25 472 1
## 26 615 0
## 27 606 0
## 28 634 1
## 29 633 1
## 30 688 0
## 31 694 0
## 32 671 0
## 33 339 1
## 34 692 1
## 35 605 0
## 36 640 0
## 37 340 1
## 38 692 1
## 39 642 0
## 40 676 0
## 41 495 1
## 42 601 1
## 43 370 0
## 44 687 1
## 45 688 1
## 46 356 0
## 47 358 0
## 48 692 0
## 49 386 1
## 50 438 1
## 51 679 1
## 52 359 0
## 53 692 0
## 54 336 1
## 55 579 0
## 56 653 0
## 57 486 1
## 58 585 1
## 59 333 0
## 60 356 0
## 61 700 0
## 62 613 1
## 63 649 1
## 64 637 0
## 65 609 0
## 66 608 1
## 67 385 0
## 68 652 0
## 69 691 0
## 70 349 1
## 71 689 0
## 72 698 0
## 73 338 1
## 74 338 1
## 75 673 1
## 76 347 1
## 77 364 0
## 78 344 1
## 79 693 1
## 80 538 1
## 81 361 0
## 82 555 0
## 83 673 0
## 84 341 1
## 85 646 0
## 86 337 1
## 87 342 1
## 88 602 0
## 89 690 0
## 90 359 1
## 91 578 1
## 92 548 1
## 93 340 0
## 94 691 0
## 95 361 1
## 96 471 1
## 97 681 0
## 98 694 0
## 99 353 1
## 100 649 0
## 101 477 0
## 102 592 1
## 103 560 0
## 104 660 1
## 105 679 0
## 106 672 0
## 107 352 1
## 108 620 1
## 109 371 0
## 110 693 0
## 111 355 1
## 112 675 1
## 113 669 0
## 114 528 1
## 115 688 1
## 116 349 1
## 117 344 0
## 118 656 0
## 119 663 1
## 120 683 1
## 121 418 0
## 122 685 0
## 123 403 1
## 124 679 1
## 125 377 0
## 126 346 1
## 127 694 1
## 128 365 0
## 129 629 0
## 130 360 1
## 131 686 1
## 132 356 0
## 133 698 0
## 134 364 1
## 135 688 1
## 136 356 0
## 137 510 0
## 138 341 1
## 139 345 1
## 140 336 0
## 141 669 1
## 142 625 1
## 143 345 0
## 144 684 0
## 145 686 0
## 146 342 1
## 147 336 1
## 148 577 0
## 149 665 0
## 150 352 1
## 151 686 1
## 152 673 0
## 153 364 1
## 154 342 1
summary(data) #print basic descriptive statistics
## Correct Group
## Min. :326.0 Min. :0.0000
## 1st Qu.:356.5 1st Qu.:0.0000
## Median :589.5 Median :1.0000
## Mean :532.0 Mean :0.5065
## 3rd Qu.:675.8 3rd Qu.:1.0000
## Max. :700.0 Max. :1.0000
For our purposes here, we can think of effect sizes in two ways: the distance between the mean and some theoretically important value in standard deviation units and the difference between two subgroup means also in (pooled) standard deviation units. We compute the first of these two types of effect sizes below:
M = mean(data$Correct) #mean of all subjects
s = sd(data$Correct) #standard deviation of all subjects
d1 = (M-350)/s #effect size as # of standard deviations away from 350 (chance level)
d1
## [1] 1.235511
We have \(d = 1.24\), which is understandably large given that most subjects probably should do better than chance on this task.
The second type of effect size, which considers the difference between means for some subset of the data, is computed below:
groupA = data$Correct[data$Group==0]
groupB = data$Correct[data$Group==1]
MA = mean(groupA)
MA
## [1] 563.4737
MB = mean(groupB)
MB
## [1] 501.3077
SA = sd(groupA)
SB = sd(groupB)
s = (SA+SB)/2 #average standard deviations for each group; note: this is not exactly how the pooled standard deviation would be calculated typically
d2 = (MA-MB)/s #compute standardized effect size
d2
## [1] 0.4307027
In some contexts it may make sense to compute power after one knows the type of test, fixed sample size, significance level, plausible effect size, and standard deviation. For instance, if we plan on running a two-sample t-test with \(n = 30\), \(\alpha = .05\), \(d = 5\), and \(s = 10\), we can easily compute power with power.t.test:
power.t.test(type="two.sample",n=30,sig.level=.05,
delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 30
## delta = 5
## sd = 10
## sig.level = 0.05
## power = 0.477841
## alternative = two.sided
##
## NOTE: n is number in *each* group
Here our power is around 48%, meaning that we have a 48% chance of rejecting the null hypothesis if our hypothesized effect exists.
It’s useful to consider what happens to power when the parameters in power.t.test change. For instance, here’s what happens when sample size increases:
power.t.test(type="two.sample",n=5,sig.level=.05,
delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 5
## delta = 5
## sd = 10
## sig.level = 0.05
## power = 0.1038399
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=30,sig.level=.05,
delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 30
## delta = 5
## sd = 10
## sig.level = 0.05
## power = 0.477841
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=60,sig.level=.05,
delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 60
## delta = 5
## sd = 10
## sig.level = 0.05
## power = 0.7752644
## alternative = two.sided
##
## NOTE: n is number in *each* group
…and here’s what happens when the expected effect size increases:
power.t.test(type="two.sample",n=40,sig.level=.05,
delta=1,sd=10)
##
## Two-sample t test power calculation
##
## n = 40
## delta = 1
## sd = 10
## sig.level = 0.05
## power = 0.06447851
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=40,sig.level=.05,
delta=3,sd=10)
##
## Two-sample t test power calculation
##
## n = 40
## delta = 3
## sd = 10
## sig.level = 0.05
## power = 0.2627742
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=40,sig.level=.05,
delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 40
## delta = 5
## sd = 10
## sig.level = 0.05
## power = 0.5981316
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=40,sig.level=.05,
delta=10,sd=10)
##
## Two-sample t test power calculation
##
## n = 40
## delta = 10
## sd = 10
## sig.level = 0.05
## power = 0.9929848
## alternative = two.sided
##
## NOTE: n is number in *each* group
…and here’s what happens when we change \(\alpha\), the significance level:
power.t.test(type="two.sample",n=40,sig.level=.1,delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 40
## delta = 5
## sd = 10
## sig.level = 0.1
## power = 0.7162549
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=40,sig.level=.05,delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 40
## delta = 5
## sd = 10
## sig.level = 0.05
## power = 0.5981316
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=40,sig.level=.01,delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 40
## delta = 5
## sd = 10
## sig.level = 0.01
## power = 0.3493075
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=40,sig.level=.001,delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 40
## delta = 5
## sd = 10
## sig.level = 0.001
## power = 0.1289781
## alternative = two.sided
##
## NOTE: n is number in *each* group
…but power remains unchanged when the standardized effect size remains constant (even if both \(d\) and \(s\) change):
power.t.test(type="two.sample",n=30,sig.level=.05,
delta=0.5,sd=1)
##
## Two-sample t test power calculation
##
## n = 30
## delta = 0.5
## sd = 1
## sig.level = 0.05
## power = 0.477841
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=30,sig.level=.05,
delta=5,sd=10)
##
## Two-sample t test power calculation
##
## n = 30
## delta = 5
## sd = 10
## sig.level = 0.05
## power = 0.477841
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",n=30,sig.level=.05,
delta=50,sd=100)
##
## Two-sample t test power calculation
##
## n = 30
## delta = 50
## sd = 100
## sig.level = 0.05
## power = 0.477841
## alternative = two.sided
##
## NOTE: n is number in *each* group
One of the most important uses for power.t.test is to compute the required sample size to achieve a given level of power (usually 80% or higher). This can be done by swapping “n” in power.t.test for “power” as parameters:
power.t.test(type="two.sample",power=.80,sig.level=.05,
delta=2,sd=10)
##
## Two-sample t test power calculation
##
## n = 393.4067
## delta = 2
## sd = 10
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group
power.t.test(type="two.sample",power=.80,sig.level=.05,
delta=8,sd=10)
##
## Two-sample t test power calculation
##
## n = 25.52463
## delta = 8
## sd = 10
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group