A Function for Meng's Test of the Heterogeneity of Correlated Correlation Coefficients

Here's a first draft of a simple function for the test of heterogeneity of correlated correlation coefficients in Meng, Rosenthal and Rubin (1992), which I blogged about way back when. There are no tests included (e.g. tests that input is numeric) and the code is undoubtedly a bit flabby, but it works perfectly well.

It's a potentially useful test for many, but doesn't appear to be provided in any of the mainstream statistical software. 

Importantly, the function works for any number of correlated correlation coefficients and also allows for contrast tests. 

As input, the function takes a vector of correlations, the median correlation between the variables, sample size, and contrast weights (if required). 

	meng_test <- function(corrs, mdn_r, n, contrast_weights = rep(0, length(corrs))){
	library(psych)
	cors <- corrs # put correlated correlation coefficients here
	zs <- fisherz(cors) #fisherz transforms
	rs <- fisherz2r(zs) #transformations back
	mz <- mean(zs) #mean fisher z
	zmz <- zs-mz #fisher z minus mean fisher z
	zmz2 <- zmz*zmz #fisher z minus mean fisher z squared
	cw <- contrast_weights #put contrast weights in here, MUST SUM TO ZERO
	cwzs <- cw*zs #contrast weights multiplied by fisher z
	cw2 <-cw^2 #contrast weight squared
	r <- corrs # put correlation coefficients (Xs -> Y) in here
	rc <- matrix(r.con(r,n),ncol=2) #confidence intervals of r
	t <- r*sqrt(n-2)/sqrt(1-r^2) #t
	p <- (1-pt(t,n-2))/2 #p
	square <-r*r #r squared
	r.rc <- data.frame(r=r,z=fisherz(r),lower=rc[,1],upper=rc[,2],t=t,p=p,rsquare=square,zmz=zmz, zmz2 = zmz2, cw=cw, cwzs=cwzs, cw2=cw2)
	#average r square
	avrsqur <- mean(square)
	#sum of fishers z- mean fishers z
	sumzmz <- sum(zmz2)
	#enter median intercorrelation of Xs
	mi <- mdn_r
	f <- (1-mi)/(2*(1-avrsqur)) #find f
	h <- (1-(f*avrsqur))/(1-avrsqur) #find h
	#find chi-square (k-1) # THIS CHI-SQUARE IS A TEST OF THE HETEROGENEITY OF THE CORRELATED COEFFICIENTS
	chisqu1 <- (n-3)*sumzmz
	chisqu2 <- (1-mi)*h
	chisqu <- chisqu1/chisqu2
	df <- length(cors) - 1
	p <- pchisq(chisqu, df, lower.tail=FALSE)
	scwzs <- sum(cwzs) #sum of contrast weights fisher z transform (first part of test) #2.025 in example
	denom1 <- sum(cw2) #sum of contrast weights squared
	denom2 <- 1 - mi #1 - median intercorrelation of X's
	ct1 <- (n-3)/ (denom1*denom2*h)
	ct2 <- sqrt(ct1)
	contrastz <- scwzs*ct2
	p2 <- pnorm(-abs(contrastz))
	return(list(csquare = chisqu, df = df, pval = p, zcontrast = contrastz, zpval = p2))
	}
	# confirm it works with the example from Meng et al. (1992)
	meng_test(corrs = c(.63, .53, .54, -.03), mdn_r = .37, n = 15)

view raw gistfile1.txt hosted with ❤ by GitHub