# Function to calculate a truely pooled correlation
# matrix from two correlation matrices given the means,
# standard-deviations (as estimated population para-
# meters), numbers of (valid) cases (per correlation
# matrix, i.e. with listwise deletion of missing cases),
# and correlation matrices of the variables of both
# samples:

corr.tot = function(m1,s1,n1,rmat1,
                    m2,s2,n2,rmat2)
{
  n = n1+n2
  m = (m1*n1+m2*n2)/n
  ss1 = sqrt(s1**2*(n1-1)/n1)
  ss2 = sqrt(s2**2*(n2-1)/n2)
  s = sqrt(1/n * (n1*ss1**2 + n2*ss2**2) +
           1/n * (n1*(m1-m)**2 + n2*(m2-m)**2))
  rmat = (n1*rmat1*(ss1 %*% t(ss1)) +
          n2*rmat2*(ss2 %*% t(ss2)) +
          n1*(m1 %*% t(m1)) +
          n2*(m2 %*% t(m2)) -
          n*(m %*% t(m)))/(n*(s %*% t(s)))
  pooled = list('mean'=m,
                'sd'=sqrt(n/(n-1))*s,
                'n'=n,
                'rmat'=rmat)
  pooled
}

# example (de-comment the following lines):
#
# m1 = c(12.5,13.453,0.98)
# s1 = c(8.2,11.4,0.6034)
# n1 = 48
# rmat1 = matrix(c(1.000, .421, .123,
#                   .421,1.000,-.346,
#                   .123,-.346,1.000),nrow=3)
#
# m2 = c(9.2,14.541,1.02)
# s2 = c(8.2,9.3,0.5034)
# n2 = 72
# rmat2 = matrix(c(1.000, .395, .098,
#                   .395,1.000,-.288,
#                   .098,-.288,1.000),nrow=3)
#
# corr.tot(m1,s1,n1,rmat1,m2,s2,n2,rmat2)