Significance of correlation with range restriction in x

Question

Background

Range restriction can happen in various study designs and it affects the results of bivariate correlation or in multiple regression, for example. In some cases the range restriction may be done for practical reasons. For example, in medicine only patients with age >18 and <50 may be included because of insurance costs. More importantly, in some cases the range restriction is not obvious and may remain unnoticed. To list some examples, this can happen in prediction of job performance, biology or in the evaluation of the University Clinical Aptitude Test (UCAT). Surely, there are many other situations where range restriction can affect the results.

Consequences and correction of range restriction

The range restriction causes the standard deviation of the variable x to shrink which reduces the correlation of x to some variable y. Analogue, range restriction can also affect the coefficients in multiple regression. There are suggestions how to adjust for range restriction in multiple regression or in bivariate correlation.

My question

Using the approach from this paper indeedly gives me correlation coefficients that are closer to the expected values. But the adjustment for the p value seems not to work because using this adjustment gives me p values that are even more different from the expected p values in the unrestricted population compared to the p values I would get using no adjustment for range restriction. The difference is small but I still don't see a point to use the p value for the adjusted corelation if the result seems to be even more biased. Why is this?

My R code

I used to following R code to apply the adjustments for range restriction as described in the source.

# defining the functions
##########
# random variable with fixed mean and sd. code as decribed here:
# https://stackoverflow.com/questions/18919091/generate-random-numbers-with-fixed-mean-and-sd
rnorm2 <- function(n,mean,sd) { mean+sd*scale(rnorm(n)) }

# functions based on literature provided
# determining true correlation (page 53, equation 1 from source provided)
rxy_adjust <- function(k, rxy){
  rxx <- 1
  ryy <- 1
  rxy_adjusted <- (k* rxy)/ sqrt(rxx*ryy - rxy**2 + k**2*rxy**2)
  return(rxy_adjusted)}

# eastimating error variance for fixed reliabilities in x and y (page 54,    equation 10 from source provided)
errorvar_adjust <- function(n, rxy, k){
  errv <- (k**2*(1 - rxy**2)**2)/
    (n - 1)*(1 - rxy**2 + k**2*rxy**2)**3
  return(errv)}
##########

# choosing parameters
##########
# choosing the unrestricted sd of x
sdx_unrestricted_choosen <- 30

# choosing the unrestricted correlation between x and y
rxy_unrestricted_choosen <- 0

# choosing number of "simulations" being run
nsim <- 1000000

# choosing saplz size
n <- 100
##########

# run simulations
##########
results <- sapply(1:nsim, function(x){
  # generate unrestricted variables
  m_unrestricted <-   matrix(c(rnorm2(n, 50, sdx_unrestricted_choosen), rnorm(n)), ncol= 2)

  # name the matrix columns
  colnames(m_unrestricted) <- c("x", "y")

  # choosing only objects that lay inside a range >0 and < 80 of x for the restricted variables
  m_restricted <- m_unrestricted[m_unrestricted[ ,"x"] > 0 & m_unrestricted[ ,"x"] < 80, ]

  # correlations
  # calculate correlation for unrestricted data (always as set before)
  rxy_unrestricted <- cor(m_unrestricted)[1,2]
  # calculate correlation for restricted data
  rxy_restricted <- cor(m_restricted)[1,2]

  # calculate correlation which adjusts for the range restriction in x
  # 1. what is the restricted sd of x? (following 4 steps, 1-4)
  sd_restricted <- sd(m_restricted[ , "x"])
  # 2. determine u (the ratio between restricted and unrestricted sd of x)
  # as described in text on page 53 from source provided
  u <- sd_restricted/ sdx_unrestricted_choosen
  # 3. k, defined as 1/u
  # as described in text on page 53 from source provided
  k <- 1/u
  # 4. apply function for adjustment in range restriction in x to estimate adjusted
  # correlation between x and y
  rxy_adjusted <- rxy_adjust(k= k, rxy= rxy_restricted)


  # p values
  # p value for the unrestricted correlation
  p_unrestricted <- cor.test(m_unrestricted[ ,1], m_unrestricted[ ,2])$p.value
  # p value for the restricted correlation
  p_restricted <- cor.test(m_restricted[ ,1], m_restricted[ ,2])$p.value
  # p value for the adjusted correlation (following 3 steps, 1-3)
  # 1. apply function to estimate error varaince
  errorvar_adjusted <- errorvar_adjust(n= nrow(m_restricted),
                                       rxy= rxy_restricted,
                                       k= k)
  # 2. calculate the z value
  zval <- rxy_adjusted/ sqrt(errorvar_adjusted)
  # 3. determine p value for z
  p_adjusted <- (1 - pnorm(abs(zval))) * 2


  # make a vector with all values
  results_temp <- c(rxy_unrestricted, rxy_restricted, rxy_adjusted,
                    p_unrestricted, p_restricted, p_adjusted)
  # return vector with all values
  return(results_temp)})
##########

# results
##########
# rotating the matrix of results
results <- t(results)

# naming the columns
colnames(results) <- c("rxy_unrestricted", "rxy_restricted", "rxy_adjusted",
                       "p_unrestricted", "p_restricted", "p_adjusted")


# comparing correlations
# absolute difference between the unrestricted and the restricted correlations
diff_rxy_unrestricted_restricted <- abs(results[ , "rxy_unrestricted"] - results[ , "rxy_restricted"])

# absolute difference between the unrestricted and the adjusted correlations
diff_rxy_unrestricted_adjusted <- abs(results[ , "rxy_unrestricted"] - results[ , "rxy_adjusted"])

# median deviations from the true (unrestricted) correlations
median(diff_rxy_unrestricted_restricted); median(diff_rxy_unrestricted_adjusted)
# 0.06474827 and 0.08936159    
# result: the unadjusted correlation coefficients are closer to the unrestricted coefficients!


# p values
# type 1 error of p_unrestricted
sum(results[ , "p_unrestricted"] < .05)/ nrow(results)
# 0.049903

# type 1 error of p_restricted
sum(results[ , "p_restricted"] < .05)/ nrow(results)
# 0.050148

# type 1 error of p_adjusted
sum(results[ , "p_adjusted"] < .05)/ nrow(results)
# 0.033725
##########

Some notes

To my understanding the issue is not the code and that is why I ask this question here and not on stackoverflow. For those who don't want to follow the link here is the paper I try to realize in R: Paper: Raju, N. S., & Brand, P. A. (2003). Determining the significance of correlations corrected for unreliability and range restriction. Applied Psychological Measurement, 27(1), 52-71.

EDIT

I have changed the R code. Before, I set the unrestricted correlation to .5 and the results indicated that the adjustment worked in regard to the correlation coefficient but not to the p value (as described above). After Reading the answer from Frans Rodenburg I adapted the code; now two random variables are generated where expecter correlation is zero. Now the correlation coefficient of the unadjusted correlation is closer to the unrestricted data than the coefficient of the adjusted analysis! Further, the type 1 error is worse for the adjusted analysis (0.033725) than for the unadjusted analysis (0.050148) or the unrestricted data (0.049903). This would mean that running a usual pearson correlation would give better results than the adjustment, hence, I am wondering even more about the use of the adjustment in the cited paper.

Answer 1

There are two issues:

• Since you simulate data, you have the true correlation! So use it, don't compare your method's performance to the estimated correlation of the unrestricted data;
• There is no such thing as a true $p$-value. In the population, the difference either is or isn't. When sampling, the $p$-value is a random variable that depends on your sample size. In your restricted correlation test, you are using a proper subset of the data. Hence you will have fewer degrees of freedom and thus less compelling evidence against the null. This is reflected by the higher $p$-value, so nothing is wrong with your implementation there. In fact, it still comes to the same conclusion at commonly used levels of significance.

The bias and variance of your method can be estimated solely from comparisons with the true correlation (e.g. using bootstrapping). The only reason you'd want to compute a $p$-value at all, is to compare the type I and type II error rates of the adjusted and unadjusted approach. The actual $p$-values of individual tests don't matter, as long as the test comes to the right conclusion more often.

Then again, even if the adjusted method has a higher type II error, it might still be preferred for inference if the goal is to obtain more accurate estimates of correlation, so I would do away with the 'significance' of the correlation entirely. (If you don't, then you could simulate with true correlations of $\rho=0$ to determine the type I error rate.)