// Apr 2 2012 14:03:51 // Calculate the correlation between two distance matrices // Correlation is between the lower half (including diagonal) values mata: real vector rankorder(real vector dv) { real vector l1, l2 real matrix workspace l1 = (1::rows(dv)) workspace = (dv, l1) l2 = workspace[order(workspace,(1,2)),(2,1)] return(l1[order(l2,1),]) } // invorder(order(x,1)) == rankorder(x) real corrsqm(string v1, string v2, scalar sp) { real matrix m1, m2, cpair m1 = st_matrix(v1) m2 = st_matrix(v2) if (!issymmetric(m1) | !issymmetric(m2)) { return(-1) } else { m1 = vech(m1) m2 = vech(m2) if (sp) { m1 = rankorder(m1) m2 = rankorder(m2) m1 == m2 } cpair = correlation((m1, m2), 1) return(cpair[1,2]) } } real corrsqmnodiag(string v1, string v2, scalar sp) { real matrix m1, m2, cpair m1 = st_matrix(v1) m2 = st_matrix(v2) if (!issymmetric(m1) | !issymmetric(m2)) { return(-1) } else { nr = rows(m1) m1 = m1[2..nr,1..(nr-1)] m2 = m2[2..nr,1..(nr-1)] m1 = vech(m1) m2 = vech(m2) // m1 // m2 if (sp) { m1 = rankorder(m1) m2 = rankorder(m2) m1 == m2 } cpair = correlation((m1, m2), 1) return(cpair[1,2]) } } end capture program drop corrsqm program define corrsqm, rclass syntax namelist, [SPEarman NODiag] local mat1 : word 1 of `namelist' local mat2 : word 2 of `namelist' local sp 0 if ("`spearman'" != "") local sp 1 if (rowsof(`mat1') != rowsof(`mat2')) { di in red "Matrices are not the same size" exit } else { tempname retval if ("`nodiag'" == "") { mata: st_numscalar("`retval'",corrsqm("`mat1'", "`mat2'",`sp')) } else { mata: st_numscalar("`retval'",corrsqmnodiag("`mat1'", "`mat2'",`sp')) } if (`retval' == -1) { di in red "At least one matrix not symmetric" } else { di "VECH correlation between `mat1' and `mat2': " %6.4f `retval' } if ("`nodiag'" != "") { di "Diagonal suppressed" } return scalar rho = `retval' } end