-
Notifications
You must be signed in to change notification settings - Fork 0
/
2.2_d_p_power_dist.R
201 lines (192 loc) · 10.7 KB
/
2.2_d_p_power_dist.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
if (!require(shiny)) {install.packages('shiny')}
library(shiny)
if (!require(shinythemes)) {install.packages('shinythemes')}
library(shinythemes)
ui <- fluidPage(theme=shinytheme("flatly"),
titlePanel("Distribution of Cohen's d, p-values, and power curves for an independent two-tailed t-test"),
sidebarLayout(
sidebarPanel(numericInput("N", "Participants per group:", 50, min = 1, max = 1000),
sliderInput("d", label = HTML("Cohen's d (δ)"), min = 0, max = 2, value = 0.5, step= 0.01),
sliderInput("p_upper", "alpha, or p-value (upper limit):", min = 0, max = 1, value = 0.05, step= 0.005),
uiOutput("p_low"),
h4(textOutput("pow0")),br(),
h4(textOutput("pow1")),br(),
h4(textOutput("pow2")),br(),
h4("The three other plots indicate power for a range of alpha levels (top right), sample sizes per group (bottom left), and effect sizes (bottom right). The bottom right figure illustrates the point that when the true effect size of a study is unknown, the power of a study is best thought of as a curve, not as a single value."),br()
),
mainPanel(
plotOutput("plot_d", width = "1004px", height = "500px"),
splitLayout(style = "border: 1px solid silver:", cellWidths = c(500,500),cellHeights = c(800,800),
plotOutput("pdf"),
plotOutput("cdf")
),
splitLayout(style = "border: 1px solid silver:", cellWidths = c(500,500),
plotOutput("power_plot"),
plotOutput("power_plot_d")
),
h4("Get the code at ", a("GitHub", href="https://github.com/Lakens/p-curves"))
)
)
)
server <- function(input, output) {
#surpress warnings
options(warn = -1)
output$pdf <- renderPlot({
N<-input$N
d<-input$d
p<-0.05
p_upper<-input$p_upper+0.00000000000001
p_lower<-input$p_lower+0.00000000000001
# if(p_lower==0){p_lower<-0.0000000000001}
ymax<-25 #Maximum value y-scale (only for p-curve)
#Calculations
se<-sqrt(2/N) #standard error
ncp<-(d*sqrt(N/2)) #Calculate non-centrality parameter d
#p-value function
pdf2_t <- function(p) 0.5 * dt(qt(p/2,2*N-2,0),2*N-2,ncp)/dt(qt(p/2,2*N-2,0),2*N-2,0) + dt(qt(1-p/2,2*N-2,0),2*N-2,ncp)/dt(qt(1-p/2,2*N-2,0),2*N-2,0)
par(bg = "aliceblue")
plot(-10,xlab="P-value", ylab="Density", axes=FALSE,
main=substitute(paste("P-value distribution for ", delta == d," and N =",N)), xlim=c(0,1), ylim=c(0, ymax))
abline(v = seq(0,1,0.1), h = seq(0,ymax,5), col = "lightgray", lty = 1)
axis(side=1, at=seq(0,1, 0.1), labels=seq(0,1,0.1))
axis(side=2)
cord.x <- c(p_lower,seq(p_lower,p_upper,0.001),p_upper)
cord.y <- c(0,pdf2_t(seq(p_lower, p_upper, 0.001)),0)
polygon(cord.x,cord.y,col=rgb(1, 0, 0,0.5))
curve(pdf2_t, 0, 1, n=1000, col="black", lwd=2, add=TRUE)
})
output$cdf <- renderPlot({
N<-input$N
d<-input$d
p_upper<-input$p_upper
p_lower<-input$p_lower
ymax<-25 #Maximum value y-scale (only for p-curve)
#Calculations
se<-sqrt(2/N) #standard error
ncp<-(input$d*sqrt(N/2)) #Calculate non-centrality parameter d
cdf2_t<-function(p) 1 + pt(qt(p/2,2*N-2,0),2*N-2,ncp) - pt(qt(1-p/2,2*N-2,0),2*N-2,ncp)
par(bg = "aliceblue")
plot(-10,xlab="Alpha", ylab="Power", axes=FALSE,
main=substitute(paste("Power for independent t-test with ", delta == d," and N =",N)), xlim=c(0,1), ylim=c(0, 1))
abline(v = seq(0,1,0.1), h = seq(0,1,0.1), col = "lightgray", lty = 1)
axis(side=1, at=seq(0,1, 0.1), labels=seq(0,1,0.1))
axis(side=2)
# cord.x <- c(p_lower,seq(p_lower,p_upper,0.001),p_upper)
# cord.y <- c(0,cdf2_t(seq(p_lower, p_upper, 0.001)),0)
# polygon(cord.x,cord.y,col=rgb(1, 0, 0,0.5))
curve(cdf2_t, 0, 1, n=1000, col="black", lwd=2, add=TRUE)
points(x=p_upper, y=(1 + pt(qt(input$p_upper/2,2*N-2,0),2*N-2,ncp) - pt(qt(1-input$p_upper/2,2*N-2,0),2*N-2,ncp)), cex=2, pch=19, col=rgb(1, 0, 0,0.5))
})
output$power_plot <- renderPlot({
N<-input$N
d<-input$d
p_upper<-input$p_upper
ncp<-(input$d*sqrt(N/2)) #Calculate non-centrality parameter d
plot_power <- (function(d, N, p_upper){
ncp <- d*(N*N/(N+N))^0.5 #formula to calculate t from d from Dunlap, Cortina, Vaslow, & Burke, 1996, Appendix B
t <- qt(1-(p_upper/2),df=(N*2)-2)
1-(pt(t,df=N*2-2,ncp=ncp)-pt(-t,df=N*2-2,ncp=ncp))
}
)
par(bg = "aliceblue")
plot(-10,xlab="sample size (per condition)", ylab="Power", axes=FALSE,
main=substitute(paste("Power for independent t-test with ", delta == d)), xlim=c(0,N*2), ylim=c(0, 1))
abline(v = seq(0,N*2, (2*N)/10), h = seq(0,1,0.1), col = "lightgray", lty = 1)
axis(side=1, at=seq(0,2*N, (2*N)/10), labels=seq(0,2*N,(2*N)/10))
axis(side=2, at=seq(0,1, 0.2), labels=seq(0,1,0.2))
curve(plot_power(d=d, N=x, p_upper=p_upper), 3, 2*N, type="l", lty=1, lwd=2, ylim=c(0,1), xlim=c(0,N), add=TRUE)
points(x=N, y=(1 + pt(qt(input$p_upper/2,2*N-2,0),2*N-2,ncp) - pt(qt(1-input$p_upper/2,2*N-2,0),2*N-2,ncp)), cex=2, pch=19, col=rgb(1, 0, 0,0.5))
})
output$power_plot_d <- renderPlot({
N<-input$N
d<-input$d
p_upper<-input$p_upper
ncp<-(input$d*sqrt(N/2)) #Calculate non-centrality parameter d
plot_power_d <- (function(d, N, p_upper)
{
ncp <- d*(N*N/(N+N))^0.5 #formula to calculate t from d from Dunlap, Cortina, Vaslow, & Burke, 1996, Appendix B
t <- qt(1-(p_upper/2),df=(N*2)-2)
1-(pt(t,df=N*2-2,ncp=ncp)-pt(-t,df=N*2-2,ncp=ncp))
}
)
par(bg = "aliceblue")
plot(-10,xlab=substitute(paste("Cohen's ", delta == d)), ylab="Power", axes=FALSE,
main=substitute(paste("Power for independent t-test with N =",N," per group")), xlim=c(0,2), ylim=c(0, 1))
abline(v = seq(0,2, 0.2), h = seq(0,1,0.1), col = "lightgray", lty = 1)
axis(side=1, at=seq(0,2, 0.2), labels=seq(0,2,0.2))
axis(side=2, at=seq(0,1, 0.2), labels=seq(0,1,0.2))
curve(plot_power_d(d=x, N=N, p_upper=p_upper), 0, 2, type="l", lty=1, lwd=2, ylim=c(0,1), xlim=c(0,2), add=TRUE)
points(x=d, y=(1 + pt(qt(input$p_upper/2,2*N-2,0),2*N-2,ncp) - pt(qt(1-input$p_upper/2,2*N-2,0),2*N-2,ncp)), cex=2, pch=19, col=rgb(1, 0, 0,0.5))
})
# make dynamic slider
output$p_low <- renderUI({
sliderInput("p_lower", "p-value (lower limit):", min = 0, max = input$p_upper, value = 0, step= 0.005)
})
output$pow0 <- renderText({
p_upper<-input$p_upper
N<-input$N
d<-input$d
VARd<-((N+N)/(N*N))+(d^2/(2*(N+N)))
SEd<-sqrt(VARd)
crit_d<-abs(qt(p_upper/2, (N*2)-2))/sqrt(N/2)
paste("On the top, you can see the distribution of Cohen's d assuming a true effect size of d =",d,"illustrated by the black line. Only observed effect sizes larger than d =",round(crit_d,2),"will be statisically significant with",N,"observations per group. The distribution assuming a Cohen's d of 0 is illustrated by the grey curve. Red areas illustrates Type 1 errors, the blue area illustrates the Type 2 error rate. The distribution has a standard error of",round(SEd,3))
})
output$pow1 <- renderText({
N<-input$N
d<-input$d
paste("On the right, you can see the p-value distribution for a two-sided independent t-test with",N,"participants in each group, and a true effect size of ", d)
})
output$pow2 <- renderText({
N<-input$N
d<-input$d
p_upper<-input$p_upper
p_lower<-input$p_lower
ncp<-(input$d*sqrt(N/2)) #Calculate non-centrality parameter d
p_u<-1 + pt(qt(p_upper/2,2*N-2,0),2*N-2,ncp) - pt(qt(1-p_upper/2,2*N-2,0),2*N-2,ncp) #two-tailed
p_l<-1 + pt(qt(p_lower/2,2*N-2,0),2*N-2,ncp) - pt(qt(1-p_lower/2,2*N-2,0),2*N-2,ncp) #two-tailed
paste("The statistical power based on an alpha of",p_upper,"and assuming the true effect size is d =",d,"is",100*round((1 + pt(qt(input$p_upper/2,2*N-2,0),2*N-2,ncp) - pt(qt(1-input$p_upper/2,2*N-2,0),2*N-2,ncp)),digits=4),"%. In the long run, you can expect ",100*round(p_u-p_l, 4),"% of p-values to fall in the selected area between p = ",p_lower,"and p = ",p_upper,".")
})
output$plot_d <- renderPlot({
N<-input$N
d<-input$d
p_upper<-input$p_upper
ncp<-(input$d*sqrt(N/2)) #Calculate non-centrality parameter d
low_x<--2
high_x<-3
#calc d-distribution
x=seq(low_x,high_x,length=10000) #create x values
d_dist<-dt(x*sqrt(N/2),df=(N*2)-2, ncp = ncp)*sqrt(N/2) #calculate distribution of d based on t-distribution
#Set max Y for graph
y_max<-max(d_dist)+1
#create plot
par(bg = "aliceblue")
d = round(d,2)
plot(-10,xlim=c(low_x,high_x), ylim=c(0,y_max), xlab=substitute(paste("Cohen's ", delta == d)), ylab="Density",main=substitute(paste("Distribution of Cohen's ", delta == d,", N = ",N)))
#abline(v = seq(low_x,high_x,0.1), h = seq(0,0.5,0.1), col = "lightgray", lty = 1)
axis(side = 1, at = seq(low_x,high_x,0.1), labels = FALSE)
lines(x,d_dist,col='black',type='l', lwd=2)
#add d = 0 line
d_dist<-dt(x*sqrt(N/2),df=(N*2)-2, ncp = 0)*sqrt(N/2)
lines(x,d_dist,col='grey',type='l', lwd=1)
#Add Type 1 error rate right
crit_d<-abs(qt(p_upper/2, (N*2)-2))/sqrt(N/2)
y=seq(crit_d,10,length=10000)
z<-(dt(y*sqrt(N/2),df=(N*2)-2)*sqrt(N/2)) #determine upperbounds polygon
polygon(c(crit_d,y,10),c(0,z,0),col=rgb(1, 0, 0,0.5))
#Add Type 1 error rate left
crit_d<--abs(qt(p_upper/2, (N*2)-2))/sqrt(N/2)
y=seq(-10, crit_d, length=10000)
z<-(dt(y*sqrt(N/2),df=(N*2)-2)*sqrt(N/2)) #determine upperbounds polygon
polygon(c(y,crit_d,crit_d),c(z,0,0),col=rgb(1, 0, 0,0.5))
#Add Type 2 error rate
crit_d<-abs(qt(p_upper/2, (N*2)-2))/sqrt(N/2)
y=seq(-10,crit_d,length=10000)
z<-(dt(y*sqrt(N/2),df=(N*2)-2, ncp=ncp)*sqrt(N/2)) #determine upperbounds polygon
polygon(c(y,crit_d,crit_d),c(0,z,0),col=rgb(0, 0, 1,0.5))
segments(crit_d, 0, crit_d, y_max-0.8, col= 'black', lwd=2)
text(crit_d, y_max-0.5, paste("Effects larger than d = ",round(crit_d,2),"will be statistically significant"), cex = 1)
})
}
shinyApp(ui = ui, server = server)
# Daniel Lakens, 2019.
# This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. https://creativecommons.org/licenses/by-nc-sa/4.0/