[R] paraPen in gam [mgcv 1.4-1.1] and centering constraints

Simon Wood s.wood at bath.ac.uk
Mon Feb 9 11:07:28 CET 2009

Daniel,

There's nothing built in to constraint the coefficients in this way, but it's
not too difficult to reparameterize to impose any linear constraint. The key
is to find an orthogonal basis for  null space of the constraint. Then it's
easy to re-parameterize to work in that space. Here's a simple linear model
based example ....

## simulated problem: linear model subject to constraint
## i.e y = X b + e subject to C b=0
X <- matrix(runif(40),10,4) ## a model matrix
y <- rnorm(10)              ## a response
C <- matrix(1,1,4)          ## a constraint matrix so C b = 0

## Get a basis for null space of constraint...
qrc <- qr(t(C))
Z <- qr.Q(qrc,complete=TRUE)[,(nrow(C)+1):ncol(C)]

## absorb constraint into basis...
XZ <- X%*%Z

## fit model....
b <- lm(y~XZ-1)

## back to original parameterization
beta <- Z%*%coef(b)
sum(beta) ## it works!

## If there had been a penalty matrix, S, as well,
## then it should have been transformed as follows....

S <- t(Z)%*%S%*%Z

... So provided that you have the model matrix for the term that you want to
penalize in `gam' then it's just a matter of transforming that model matrix
and corresponding penalty/ies, using a null space matrix like Z.

Note that the explicit formation of Z is not optimally efficient here, but
this won't have a noticeable impact on  the total computational cost in this
context anyway (given that mgcv is not able to make use of all that lovely
sparcity in the MRF, :-( ).

Hope this helps.

best,
Simon

On Saturday 07 February 2009 21:00, Daniel Sabanés Bové wrote:
> Dear Mr. Simon Wood, dear list members,
>
> I am trying to fit a similar model with gam from mgcv compared to what I
> did with BayesX, and have discovered the relatively new possibility of
> incorporating user-defined matrices for quadratic penalties on
> parametric terms using the "paraPen" argument. This was really a very
> good idea!
>
> However, I would like to constraint the coefficients belonging to one
> penalty matrix to sum to zero. So I would like to have the same
> centering constraint on user-defined penalized coefficient groups like
> it is implemented for the spline smoothing terms. The reason is that I
> have actually a factor coding different regions, and the penalty matrix
> results from the neighborhood structure in a Gaussian Markov Random
> Field (GMRF). So I can't choose one region as the reference category,
> because then the structure in the other regions would not contain the
> same information as before...
>
> Is there a way to constraint a group of coefficients to sum to zero?
>