Work on PR for gratia

[nicholasjclark]

Would likely require methods for:

• vcov() with same arguments as vcov.gam()
• update get_vcov() to return the Bayesian covariance matrix; shouldn't impact marginaleffects functionality
• predict(object, type = 'terms')
• predict(object, type = 'lpmatrix')
• residuals(object, type = 'working')
• weights(object)
• smoothCon() and PredictMat() functions for gp() and dynamic() terms

... actually will be easier to

1. use expand.grid() to get pred values for terms of interest. Accept an argument that allows multidimensional effects to be drawn either as heatmaps or in plot_predictions() style with faceting
2. fix all other vars to representative row idxs by replicating a single entry in original data (only filling in the vars of interest). This will work with df or list types and the values for non-focal vars won't matter as we will use 'terms' predictions
3. use predict(type = 'terms') to get partial contribution from effect of interest
4. plot by sending to gratia::draw_smooth_estimates() and with modified multidimensional plot functions, i.e. plot_smooth.bivariate_smooth_facet() perhaps