#R Code behind Figure 1 in DataColada[39] "Power Naps"
#www.datacolada.org
#Written by:Uri Simonsohn
#uws@wharton.upenn.edu
#Last update: 2015 06 17
#
#Function that finds between subject equivalent to mixed design with n=x participants, subtracting baseline before treatment, with given correlation
#for test,retest reliability;
library(pwr) #load power library
#Create the function that find the equivalent
nequiv=function(n,r,power)
{
#Sytax;
#n: actual sample size (# of participants per cell)
#r: test-retest reliability for r(baseline,dv)
#power: power of test when done between subject, ignoring baseline
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=df) #critical t-vale for p=.05 two-sided, given sample size
ncp=sqrt(n/2)*d.true #noncentrality parameter for between subject
ncp.w=ncp*sqrt((1/(2-2*r))) #noncentrality parameter for within subject (computed by reducing SE of statistic based on r)
pw=1-pt(tc,df=2*n-2,ncp=ncp.w)+pt(-tc,df=2*n-2,ncp=ncp.w) #power subtracting baseline
nw=round(pwr.t.test(power=pw,d=d.true)$n) #find n that would give that effect size the obtained within subject power
cat("\nTest-retest: r=",r, "(entered by user)")
cat("\nActual sample size n=",n," (entered by user)")
cat("\nPower=",power," (entered by user)")
cat("\nAs if sample had been n=",nw," (calculated)")
cat("\nPower obtained with subtraction=",round(pw,3)," (calculated)\n\n")
return(c(n,power,nw,round(pw,3))) #return: actual sample size, power for that with betwen subject, within n equivalent, within power.
}
#Example 1. If r=.5 power and n is the same regardless of power
nequiv(n=20,r=.5,power=.06)
nequiv(n=20,r=.5,power=.76)
#Example 2. if r=.75 as if sample (nearly) doubles regardless of effect size or sample size
nequiv(n=20,r=.75,power=.5)
nequiv(n=20,r=.75,power=.25)
nequiv(n=20,r=.75,power=.75)
nequiv(n=50,r=.75,power=.5)
nequiv(n=100,r=.75,power=.25)
nequiv(n=200,r=.75,power=.75)
#LET'S MAKE THE FIGURE NOW.
#Set of test-retest correlations to consider, between 0 and .95 (post ended up plotting up to r=.9)
rs=seq(from=0,to=.95,by=.05)
#Setup loop i=1 #counter
res20=matrix(nrow=20,ncol=4) #store results when power=20%
res70=matrix(nrow=20,ncol=4) #store results when power=70%
#Loop over possible r() values
i=1 #counter
for (r in rs) #for every value of r in the set to consider
{
res20[i,]=nequiv(n=20,r=r,power=.2) #Save nequivalent() for power=20%
res70[i,]=nequiv(n=20,r=r,power=.7) #now 70%
i=i+1
}
#Ugly chart
plot(rs,res20[,3],type='b',ylim=c(0,200),
main='Subtracting baseline when n=20 is like running...', #this is for 20% power
xlab='test-retest r',
ylab='equivalent sample size')
lines(rs,res70[,3],col='red') #in red for 70%, same exact line.
text(rs,res20[,3]+10,res20[,3]) #add value labels
abline(h=20) #add horizonatl line at n=20
#Copy paste onto excel the matrix res[20] to make nicer chart easily