[R] Fisher's r to z' transformation - help needed

[Mike White]

I am trying to use Fisher's z' transformation of the Pearson's r but the
standard error does not appear to be correct.  I have simulated an example
using the R code below.  The z' data appears to have a reasonably normal
distribution but the standard error given by the formula 1/sqrt(N-3) (from
http://davidmlane.com/hyperstat/A98696.html) gives a different results than
sd(z).  Can anyone tell me where I am going wrong?

library(amap)

## SIMULATED DATA #########################################################
p<-10
n<-3000
means<-1000*c(1:p)
SDs<-rep(100,p)
set.seed(1)
dat<-mapply(rnorm, mean=means, sd=SDs, n=n)
colnames(dat)<-paste("V",1:p, sep="")
rownames(dat)<-1:n
# calculate centroid of simulated data
dat.mean<-apply(dat,2,mean)
# calculated Pearson's r to centroid
r<-apply(dat,1,cor, y=dat.mean)
plot(density(r))
# Fisher's z' transformation
z<-0.5*log((1+r)/(1-r))
plot(density(z))
sd(z)
# [1] 0.2661779

1/sqrt(p-3)
# [1] 0.3779645

## alternatively use comparisons for all possible pairs
## Centred Pearson's r on raw data
r<-1-Dist(dat,"corr")
plot(density(r))
z<-0.5*log((1+r)/(1-r))
plot(density(z))
sd(z)
# [1] 0.2669787
1/sqrt(p-3)
# [1] 0.3779645

Many thanks
Mike White