I have a single dependent variable Y, K independent variables of interest (X1, X2... Xk), and a single independent variable of non-interest (Z).

I am interested in the relationship between Y and each X, while accounting for the effect of Z. In other words, I am interested in the partial correlations.

Further, I want to compare the slopes of each X to each other, both to evaluate whether some X's have a larger effect than others, and because I want to be able to talk about the presence/absence of a significant effect of each X without falling into the trap of "The Difference Between “Significant” and “Not Significant” is not Itself Statistically Significant".

Complicating factors:

• There is high co-linearity between many of the X variables. Because of this, and because I have a priori theoretical reasons to be interested in the non-unique variance explained by each X, I assess the effect of each X in a separate model (Y ~ X1 + Z, Y ~ X2 + Z, etc) rather than in a combined model (Y ~ X1 + X2.. + Xk + Z). My understanding is that a combined model would result in unstable (and obviously smaller) estimates of the effects of each X due to the model trying to split the shared variance.
• All of the variables have strong autocorrelation, leading to artificially inflated effect estimates (and thus significance) when assessed in a parametric framework. As such, I am using permutation testing (with autocorrelation-preserving shuffling) and the Freedman-Lane procedure (permute the residulas of the reduced model Y ~ Z) to assess the significance of the effect of X in each model.

I understand how to approach the comparison of slopes in the context of a single ANOVA with a categorical variable and multiple levels (i.e. check for an interaction between the effect of X and a dummy variable coding for category level), but I am having trouble figuring out how to compare slopes in this situation where I have:

• Multiple continuous, co-linear, autocorelated variables (X1, X2... Xk).
• Multiple separate models rather than a single combined model.
• Non-parametric permutation testing and violation of parametric assumptions about the null distribution of effect sizes.

I've looked into various methods for comparing dependent/overlapping correlations (e.g. the various tests implemented in the cocor package), but my (admittedly limited understanding) is that all of these approaches make parametric assumptions that would be violated by autocorrelation (e.g. Zou's method estimates and evaluates confidence intervals around the correlations by using Fisher's z transformation).

I feel like there is gap in my understanding that is causing me to miss an obvious solution. Any advice of solutions or even just related ideas to research would be appreciated.