Yes, I think it definitely makes sense to fit a linear regression to both plots and compare the intercepts and slopes (considering a linear fit is appropriate, be sure to check "residual by predicted" plot). You could then run a formal statistical test (perhaps an ANCOVA) to check whether the difference between the two intercepts is significant, and also whether the difference between the two slopes is significant (interaction effect).

There are many references about this, try searching "comparing two regression lines" in Google. Below are some examples

http://r-eco-evo.blogspot.com/2011/08/comparing-two-regression-slopes-by.html

https://blog.minitab.com/en/adventures-in-statistics-2/how-to-compare-regression-lines-between-different-models

One stupid pragmatic solution would be to create multiple pdps with the same grid size (not sure how to define this properly). Then substract the predictions of the pdp functions for the same grid from each other, square them and divide by the grid size. In simplicity, it is similar as the mean squared error between (MSE) two pdp functions. The the larger the MSE the larger the difference between the pdp predictions (see example). Yet, some times the grid locations (how to call this stuff ???) slightly deviate. One option would be to write an own function where you can "fix" the positions of the grid as exactly similar. You could do this by "pirating" the source code of an package or code it yourself.

As an RF model is non-linear, fitting a linear model to the pdps would be improper, unless it has a nice linear function. Yet, if your focus is directed to the intercept and slope, and the pdps show a nice linear trend, an RF model is not necessary. In this case you could better directly resort to the use of a linear model. The strength of an RF model lies in the capture of these "non-linear" patterns (and many more), but this does not necessarily mean you have to use a linear model for a linear pattern. The strength of linear model is the return of a clear mathematical equation, which represents the model in a very compressed, direct and understandable manner. At this moment, it seems that using an statistical "test" on the pdp predictions would be biased. If you would increase the grid size representing the number of samples. Then setting the grid size to 100000 would always result in p < .05.

Another pragmatic and (very) stupid, not really thought through solution I came up with at the moment would be to calculate the residuals of the two RF models Fitted1-Predicted1 and Fitted2-Predicted2 and the difference between Predicted1-Predicted2 (somehow). Then you have the variability of the two models and the difference between the two models and you could somehow compare them (not sure how though). This could also be performed for a single pdp function.

Not sure if any of these solution are viable though.

Thank you Jad and Wim for your interest and answers. Here are some more thoughts and clarifications about this from me: Of course my motivation for analyzing the data set with random forest or other tree-based model its ability to capture complex nonlinear dependences and interactions between many variables. However, say one variable of interest has a pdp plot that looks close to linear or maybe somewhat more complicated (e.g. 2nd order polynomial). Does it make sense to fit such functions to the pdp for this variable? I agree with Wim that comparing these regressions statistically does not make sense because the number of points is arbitrary due to grid size during pdp generation. Also, how reliable is the pdp plot in case some other predictor variables are substantially correlated with the variable of interest?

What do you want to do with the polynomial function?

Depends what you mean with reliable. The pdp does what it suppose to do. Best thing is to check pdps with some additional scatter or boxplots (etc.). For interference it does matter, since you do not know if causality exist, yet this is with all predictive models (unless carefully designed experiment or known causal relation).

I think the party package in R tries to compensate for all these things via permutations and also has some pdp plot functions for correlate variables (???). However, I am not sure if it actually solves the problem. I am a bit warry of the claims and neither do I notice major differences with the randomForest package but clearly highlights some important points to take in account. Hence, it is worth exploring the party package.

Otherwise, you could carefully select your samples reducing co-correlated variables (if you have a large dataset). If x1 and x2 are positively correlated select sample were x1 and x2 show the "opposite" of this trend. Otherwise, I once tried to substract the residuals of y~x1 from y and then perform a model on on y(y-resid(y~x1)) ~ x2. It does not solve the problem at all, but gives some insight on the co-correlation. Of course, the best solution would be to perform a factorial experiment.

Thank you Wim, your suggestions sound very useful and I will look into the party package. For the random forest analysis I use the ranger package and caret package to split the data into training and testing parts, tune it with cross validation, etc. Regarding the polynomial function, here is some more details: Suppose I am interested in visualizing and quantifying the dependence of the outcome variable on some variable x, and plotting the pdp for this dependence looks roughly linear, but with some "bumps" that deviate somewhat away from a linear trend. Fitting a straight line or second order polynomial to this pdp will "smooth" it and provide numerical parameters. Can these parameters then be compared with analogous parameters for a pdp for variable x from another model on the same data? The data set also contains several other predictor variables whose pdp does not look linear.

To be honest, I dont know. I do not feel qualified enough to answer this. However, I can give some food for thought.

The pdps are dependent predictions. RF models are very good in recognizing patterns in the data. The "bumps" are the results of this (e.g. gabs in the data or otherwise). You could try to tweak the nodesize of the models to smooth the pdps. Otherwise, you can also tweak the maximum number of nodes in each tree. For the latter option I have read that people often do not prefer this, although they never gave an explanation why (perhaps someone can highlight this?). The caret package gives the option to let these two tuning parameters vary between a specific range and selects the "optimal" final mode.

Additionally, if you have model 1 (x = 1.3 and 1.4) and 2, the RF model 1 will only create the split points on the data presented, it cannot select/see a split-point between 1.3 and 1.4 (e.g. 1.35) if not in the data. If, model 2 only has values 1.3 and 2.0 it can only make a split-point at 1.3 or 2.0. The slope of a model fitting between 1.4-1.3 and 2.0-1.3 might be very different. Although the range of the predictor might be similar, the split-points and thus pdp, will never be exactly similar. It seems to me that applying a polynomial is an attempt to generalize the pattern observed in the pdp (correct me if I am wrong)? The strength of the RF model lies in the bagging/bootstrapping generalizing/averaging these patterns in the data. Yet, the data is still the most important aspect. Personally I am not a person that complains about a ~5-10% difference, but I think this is inherently dependent on the topic, personal and socially accepted standards. I think it is well worth exploring (as the idea you propose seems exploratory) fitting a polynomial, but cannot guess to the meaning of it, but make sure this is clear.

If the pdps do not look linear you could try transforming the data (e.g. log, sqrt, box-cox etc.). For the RF model it matters very little, but it might show a more linear pattern. Additionally, if the pdps show a wave like pattern, this does not directly mean there the pattern is non-linear (as a pdp is dependent). It would be worth while to fit a linear model and check the residuals, and QQ-plot. The DHARMa package in R is very useful in this regard and also has a nice vignette. Perhaps Generalized Additive Models GAMs or GAMMs might fit your purpose better. Yet, I am quite unfamiliar with them as they have a lot of buttons to tweak and I never really explored them.

As always, take these points with a grain of salt.

Thank you again Wim for detailed reply and interesting thoughts. To clarify the issue with random forest tuning: I tend to use the caret package to tune the following - mtry, splitrule, and min.node.size, using repeatedcv. This is done on a training part of the data (e.g. random 75%), and performance is tested on the withheld testing part. Regarding the pdps: its not a problem of course if the pdp of interest does not look linear in this case. I am simply trying to use the pdp to better interpret the random forest model behaviors. My attempt to fit a simple parametric function to the pdp is motivated by the same thing - having numerical parameters can help to interpret the result. Does this make sense? To take the same logic to an extreme, one could potentially make a completely parametric model (e.g. polynomial regression), motivated by the pdp plots for the top random forest predictors, and fit it to the data. This would likely perform worse than the forest, but would be easily understandable. Does this seem reasonable in principle, or misguided?

pdps are useful in understanding the marginal predictions of the model. It is the shape of the curve that is of interest, positive, negative, "linear", monotonous, abrupt, etc. Since, the pdps are dependent, direct quantifiable information is difficult to extract (to my knowledge). In the first figure you see the data for mtcars (mpg~disp) the pdp follows the the relation between mpg~disp quite nicely. However, when including another variable, creating the model mpg~disp+hp. The pdp does not show the same slope and intercept anymore, hence the term dependency. Of course you could do a two dimensional pdp (or even 3), but If you have a relative high number of variables it becomes difficult. I am not sure if you can build an "acceptable" model only on the pdps.

When a pdp plot of a multivariate model (one dependent variable, and several independent variables) is performed, the pdp plot shows the relationship between the selected independent variable (x-axis) and the dependent variable (y-axis), keeping the other independent variables static (ceteris-paribus).

The pdp value is the average of all ICE curves.

In the case of a linear regression, the pdp plot will show a line. The slope of that line is the beta coefficient of the independent variable (OLS model).

I agree with what Wim Kaijser says, "It is the shape of the curve that is of interest, positive, negative, "linear", monotonous, abrupt, etc. "

I usually take as a baseline a multiple linear regression, and compare the pdp curves with another algorithm (e.g. rf, xgb, lightgbm,...), to verify that these other algorithms are not linear, adapting better to non-linearities.