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.