4
$\begingroup$

This is a follow-up to a previous question of mine, explaining the problem in more detail in the hopes of getting more precise advice.

Consider the following structured additive regression model or GAM, commonly fit by the mgcv or BayesX packages in R. I will write it a bit differently than usual, to emphasize the structure of the problem I'm considering:

$$\mu = h^{-1}(\eta)$$ $$\eta = \eta_0(\mathbf{w}) + \Delta \eta(\mathbf{x},\mathbf{z})$$ $$\Delta \eta = f_1(\mathbf{z}) + . . . + f_p(\mathbf{z}) + \mathbf{x}^T \mathbf{\gamma}$$

Here we have

  • $\mu$ is the response variable
  • $h$ is the link function
  • $\eta$ is the linear scale predictor
  • $\eta_0(\mathbf{w})$ is an offset depending on some covariates $\mathbf{w}$. We can think of $\eta_0$ as a "respected baseline model".
  • $\Delta \eta(\mathbf{x},\mathbf{z})$ is a deviation from the baseline model depending on some other covariates ($\mathbf{w}, \mathbf{x}, \mathbf{z}$ are distinct but overlapping, they are all derived from the same underlying dataset).
  • $\mathbf{z}$ are the covariates subject to semiparametric smooth regression via the unknown functions $f_1, \ldots, f_p$, which will be determined by fitting to data. (It's basically just fancy linear modeling with penalized smoothing splines).
  • $\mathbf{x}$ are the covariates for the parametric terms with unknown coefficients $\gamma$, which will be determined by fitting to data.

Suppose I elicit from experts a prior belief that the deviations $\Delta \eta$ should not be too large. I would therefore like to impose a simple Gaussian prior of known variance:

$$\Delta \eta \sim N(0, \sigma^2), \quad \sigma = 0.05 $$

How might I go about adding this prior into either mgcv or BayesX? Or should I use another package entirely? I thought it should not be terribly difficult to add this penalty into mgcv, but I can't find any way to do it. I feel like I should be able to obtain the model matrix $M$ for $\Delta \eta$ or at least $G = M^T M$, then I could just add that as a quadratic penalty to my model. But I don't see any kind of slot that this would fit into in either mgcv or BayesX. The H slot in mgcv::gam looks the closest but there is no guidance on how to use it.

Or maybe I'm thinking about this the wrong way? People normally add priors to the variables that appear "prior" to the covariates in the model DAG. Parameters that depend on the covariates, like $\eta$ and $\Delta \eta$, are typically not assigned priors. Is there something inappropriate about the way I'm conceptualizing this?

$\endgroup$
0
$\begingroup$

An interesting question. This not an answer, but a few rambling thoughts at the moment.

It sounds like you want to penalize a functional of the parameters, ie, your penalty on $\Delta \eta$ is implicitly a function $g(\beta_{11}, \beta_{12}, \dotsc, \beta_{21}, \beta_{22}, \dotsc, )$ of the parameters in the basis expansion of the smoothed fits $f_1(z) = \beta_1B_1(z) + \beta_2B_2(z) + \dotsc$.

Though mgcv has a L2-penalty (smooth.construct.re.smooth.spec) as a smoother, this applies to columns of the (possibly-augmented) design matrix, whereas you would need it to apply to the current basis expansion of the covariates.

I'm wonder if there might not be an iterative procedure in which you somehow approximate the total variation smoothness penalty $\int f''(x)^2dx$ through some combination of rotating the response and the design matrix, and weighting the L2 penalty after doing your own basis expansion? EG,you want to solve $$ argmin_{\beta} \lVert Y - B(X) \beta \rVert_2^2 + \lambda p(\beta) $$ which you approximate with $$ argmin_{\beta} \lVert Q(Y - B(X)) \beta \rVert_2^2 + \lambda \lVert D \beta \rVert_2^2 $$ for some arbitrary matrix $Q$ and diagonal matrix $D$?

$\endgroup$
  • $\begingroup$ I guess I'm not perceiving the motivation for this idea. You want to use a total variation penalty? $\endgroup$ – Paul Oct 16 '16 at 13:42
  • $\begingroup$ Errm, yeah, my bad. See edit. $\endgroup$ – Andrew M Oct 17 '16 at 21:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.