#R Code behind Figure 1 in DataColada[39] "Power Naps"
#www.datacolada.org
#Written by:Uri Simonsohn
#uws@wharton.upenn.edu
#Last update: 2016 11 05
###################################################################################################
#Corrected on 2016 11 05: Aurelien Allard identifed two errors. In the original post i had used
# d.within=d.true/sqrt(1-r^2)
#I should have used
# d.within=d.true/sqrt(1-r)
#And the power calculation had to multiply n by 2 and did not.
#The correction makes within subject manipulations have a stronger benefit, especially for more highly correlated variables.
#The blog post was updated on 2016 11 05 as well.
###################################################################################################
#Function that finds between subject equivalent to within-subject manipulation for a given correlation across measures
library(pwr)
nequiv3=function(n,r,power)
{
#Sytax;
#n: actual sample size (# of participants)
#r: test-retest reliability for r(baseline,dv)
#power: power of test when done between subject
df=2*n-2 #d.f. from between subject est
d.true=pwr.t.test(n=n,power=power)$d #effect size that gives between subject design desired power
tc=qt(.975,df=4*n-2) #critical t-vale
#d.within=d.true/sqrt(1-r^2) #WRONG: but used in origional post
d.within=d.true/sqrt(2*(1-r)) #RIGHT: corrected by Aurelien Allard in November 2016
ncp=sqrt(n/2)*d.true #noncentrality parameter for between subject
ncp.within=d.within*sqrt(n)
p.within=1-pt(tc,df=4*n-2,ncp=ncp.within)+pt(-tc,df=4*n-2,ncp=ncp.within) #power with covariate
n.within=round(2*pwr.t.test(power=p.within,d=d.true)$n) #find n that would give that effect size the obtained within subject power
return(c(n,power,n.within,round(p.within,3))) #return between n, betwen power, within n equivalent, within power.
}
#Example 1. Verify is equivalent to twice sample size if r=0
nequiv3(n=20,r=0,power=.5)
nequiv3(n=33,r=0,power=.7)
nequiv3(n=41,r=0,power=.33)
#Example 2. if r=.5 sample is n=79 regardless of effect size
nequiv3(n=20,r=.5,power=.5)
nequiv3(n=20,r=.5,power=.25)
nequiv3(n=20,r=.5,power=.75)
nequiv3(n=20,r=.9,power=.25)
#Now generate numbers behind figure in blogpost, all n equivalents for r=.0 to .095
#Set of test-retest correlations to consider
rs=seq(from=0,to=.9,by=.05)
#Setup loop
i=1 #counter
res20=matrix(nrow=19,ncol=4) #store results when power=20%
res70=matrix(nrow=19,ncol=4) #store results when power=70%
#Loop over possible r() values
for (r in rs)
{
res20[i,]=nequiv3(n=20,r=r,power=.2) #Save nequivalent() for power=20%
res70[i,]=nequiv3(n=20,r=r,power=.7) #now 70%
i=i+1
}
#Ugly chart, done for 20%, but same chart for 70%
plot(rs,res20[,3],type='b',ylim=c(0,1000))
text(rs,res20[,3]+10,res20[,3])
abline(h=20)
#Copy paste onto excel to make nicer chart easily