partial r hypothesis testing does not account for uncertainty in residualizing about correlation HOT 18 OPEN

any thoughts

No. My only understanding of correlation is: the higher, the higher - or the higher, the lower.

from correlation.

I mean, residualizing is a thing people do, without accounting for the df of the residualizing process...
I don't know...

This might be common, but it's not correct. The distribution of r/z′ is pretty well captured by df = N - 2 - np (number of variables partialed). This applies well to Pearson, Spearman, Kendall (though they each have their own SE formula). Not sure about other methdos.

from correlation.

res <- correlation::correlation(mtcars, partial = TRUE, p_adjust = "none")
as.data.frame(res[res$Parameter1=="mpg", c("Parameter2", "r", "t", "p")])
#>    Parameter2           r          t          p
#> 1         cyl -0.02326429 -0.1274582 0.89942825
#> 2        disp  0.16083460  0.8925471 0.37920361
#> 3          hp -0.21052027 -1.1795002 0.24746855
#> 4        drat  0.10445452  0.5752679 0.56940012
#> 5          wt -0.39344938 -2.3440688 0.02589012
#> 6        qsec  0.23809863  1.3427357 0.18942961
#> 7          vs  0.03293117  0.1804693 0.85799776
#> 8          am  0.25832849  1.4646374 0.15342166
#> 9        gear  0.09534261  0.5246028 0.60371458
#> 10       carb -0.05243662 -0.2876029 0.77562818

res <- ppcor::pcor(mtcars)
data.frame(r = res$estimate[-1,1],
           t = res$statistic[-1,1],
           p = res$p.value[-1,1])
#>                r          t          p
#> cyl  -0.02326429 -0.1066392 0.91608738
#> disp  0.16083460  0.7467585 0.46348865
#> hp   -0.21052027 -0.9868407 0.33495531
#> drat  0.10445452  0.4813036 0.63527790
#> wt   -0.39344938 -1.9611887 0.06325215
#> qsec  0.23809863  1.1234133 0.27394127
#> vs    0.03293117  0.1509915 0.88142347
#> am    0.25832849  1.2254035 0.23398971
#> gear  0.09534261  0.4389142 0.66520643
#> carb -0.05243662 -0.2406258 0.81217871

^{Created on 2020-04-06 by the reprex package (v0.3.0)}

good point. Should we recompute the t a posteriori based on the partial r?

Otherwise, if that's too complicated we should just remove them from the output as it doesn't make sense to add in indices that are "false"?

from correlation.

Should we recompute the t a posteriori based on the partial r?

How would you do this?

Otherwise, if that's too complicated we should just remove them from the output as it doesn't make sense to add in indices that are "false"?

I think if we cannot provide the correct stats, this would be best for the frequentist r.

However for the Bayesian r, not sure what to do as the whole posterior distribution is biased in some way.... So any index (pd, BF, pRope) would biased as well...

from correlation.

How would you do this

effectsize::r_to_t? We have the be r and the df, so in theory we should be able to retrieve the t?

But wait hold on, now that I think of it, our partial correlations are not just computed by transforming matrices, but by really computing the new correlations based on data partialization. So in fact the t and everything should be correctly computed? I.e., we do not compute the partial r from the raw r (as we could do through correlation::cor_to_pcor). So I wonder whether it is not ppcor that has issues?

from correlation.

But you don't have the correct df - you have the df of the simple correlation between two vectors. The functions doesn't know if the vectors are raw scores, or transformed scores, or residualized scores.

Which is exactly what you're saying here:

our partial correlations are not just computed by transforming matrices, but by really computing the new correlations based on data partialization.

Which is why the t and everything aren't computed correctly.
ppcor does take all of this into account, and so produces the correct t etc, which should be the same as those produced by lm() (the frequentist test for a partial correlation is equal to the frequentist for a slope in multiple regression):

res <- ppcor::pcor(mtcars)
data.frame(t = res$statistic[-1,1],
           p = res$p.value[-1,1])
#>               t          p
#> cyl  -0.1066392 0.91608738
#> disp  0.7467585 0.46348865
#> hp   -0.9868407 0.33495531
#> drat  0.4813036 0.63527790
#> wt   -1.9611887 0.06325215
#> qsec  1.1234133 0.27394127
#> vs    0.1509915 0.88142347
#> am    1.2254035 0.23398971
#> gear  0.4389142 0.66520643
#> carb -0.2406258 0.81217871

summary(lm(mpg ~ ., mtcars))$coef[-1,c(3,4)]
#>         t value   Pr(>|t|)
#> cyl  -0.1066392 0.91608738
#> disp  0.7467585 0.46348865
#> hp   -0.9868407 0.33495531
#> drat  0.4813036 0.63527790
#> wt   -1.9611887 0.06325215
#> qsec  1.1234133 0.27394127
#> vs    0.1509915 0.88142347
#> am    1.2254035 0.23398971
#> gear  0.4389142 0.66520643
#> carb -0.2406258 0.81217871

^{Created on 2020-04-06 by the reprex package (v0.3.0)}

from correlation.

All of this is also true for cor_to_ci / cor_to_p - the results of which are only true for simple correlations.

This is also somewhat true for correlation(..., multilevel = TRUE) (but I've seen multilevel correlations tested without any concern for this issue).

from correlation.

But you don't have the correct df

You're right mmmh that's tricky due to the way partial correlations are created. We could maybe retrieve the right df by counting how many variables there is in the data that underwent partialization?

from correlation.

Hmmm... that would work!

(Recall that effectsize::r_to_t was removed a while back... So you'd have to build it as an internal function here.)

res <- ppcor::pcor(mtcars)

r_to_t <- function(r, df_error) {
  r * sqrt(df_error / (1 - r ^ 2))
}

df <- res$n - nrow(res$estimate)

r_to_t(res$estimate[-1,1], df)
#>        cyl       disp         hp       drat         wt       qsec         vs 
#> -0.1066392  0.7467585 -0.9868407  0.4813036 -1.9611887  1.1234133  0.1509915 
#>         am       gear       carb 
#>  1.2254035  0.4389142 -0.2406258
res$statistic[-1,1]
#>        cyl       disp         hp       drat         wt       qsec         vs 
#> -0.1066392  0.7467585 -0.9868407  0.4813036 -1.9611887  1.1234133  0.1509915 
#>         am       gear       carb 
#>  1.2254035  0.4389142 -0.2406258

^{Created on 2020-04-06 by the reprex package (v0.3.0)}

from correlation.

Continuing 464544d here:

With r and df you can build back the t value, the p value, and the CI - so essentially I think you wouldn't need to even use stats::cor.test:

• Use cor to get the r.
• Then with n and n_parm get the df ( = n - n_parm - 1)
• Then with r and df get the t ( = r * sqrt(df_error / (1 - r ^ 2)))
• Then with t and df
• get the p.value ( = 2 * stats::pt(abs(t), df, lower.tail = FALSE) )
• get the CI with effectsize::t_to_r(t, df).

This won't work for Kendall's Tau..., but will for Pearson and Spearman.
Is patrial = TRUE supported by any of the other methods? Not sure if / how it is solvable for them....

from correlation.

Is patrial = TRUE supported by any of the other methods? Not sure if / how it is solvable for them....

Since initialization is pretty much-taking place here, it's "compatible" with all methods...

But we can start with re-writing .cor_test_freq first so that it takes n_param as an argument and remove the call to cor.test

from correlation.

Hmmmm... I think it would perhaps be better to not support this at all and thoroughly document that partialling is done by residualizing and exactly how the residualizing is done,
than to partially support it...

I mean, residualizing is a thing people do, without accounting for the df of the residualizing process...

I don't know...

from correlation.

This is also true... I mean here the df, t etc. are not wrong per se, it's more like they consider that the partialization is an independent step and that the data obtained from it is the data... (which isn't completely wrong depending on how you conceptualize it)

So yeah its quite tricky 😕

from correlation.

True... any thoughts @IndrajeetPatil @strengejacke ?

from correlation.

Is my explanation okay?

from correlation.

💥😲💥

from correlation.

@DominiqueMakowski I would change

This is why small some discrepancies are to be expected...

As this can be quite a large discrepancy, depending on:

1. How many covariates are partialled out.
2. The strength of the correlation between X, Y and the partialled out covariates (the tolerance).

from correlation.

The df are wrong here. When partial = TRUE, df needs to subtract out the number of partialed variables. Compare the psych and ppcor results:

cr <- correlation::correlation(iris[,1:4], partial = TRUE, p_adjust = "none")
crs <- summary(cr, redundant = TRUE)

pp <- ppcor::pcor(iris[,1:4])

# 2 is number of partialed variables (pp$gp)
ps <- psych::corr.p(psych::partial.r(iris[,1:4]), 150 - 2, adjust = "none")

# r
as.matrix(cr)
pp$estimate
ps$r

# t
attributes(crs)$t
pp$statistic
ps$t

# p
attributes(crs)$p
pp$p
ps$p

from correlation.