Simon,
Many thanks as always for your help.
I see and appreciate the example that you cited, but I'm having a hard
time generalizing it to a multivariate case. A bit about my context --
my dataset is response ratios; the log of a treatment over a control.
One of my explanatory variables is treatment intensity. When intensity
goes to zero, the expectation of the response ratio should also go to
zero. Here is the model that I would like to fit:
model = (ResponseRatio~
+s(as.factor(study),bs="re",by=intensity)
+s(intensity)
+s(x1),by=intensity)
+s(x2),by=intensity)
+te(x1,x2),by=intensity)
+te(x1,intensity)
)
Here is the example that you gave:
library(mgcv)
set.seed(0)
n <- 100
x <- runif(n)*4-1;x <- sort(x);
f <- exp(4*x)/(1+exp(4*x));y <- f+rnorm(100)*0.1;plot(x,y)
dat <- data.frame(x=x,y=y)
## Create a spline basis and penalty, making sure there is a knot
## at the constraint point, (0 here, but could be anywhere)
knots <- data.frame(x=seq(-1,3,length=9)) ## create knots
## set up smoother...
sm <- smoothCon(s(x,k=9,bs="cr"),dat,knots=knots)[[1]]
## 3rd parameter is value of spline at knot location 0,
## set it to 0 by dropping...
X <- sm$X[,-3] ## spline basis
S <- sm$S[[1]][-3,-3] ## spline penalty
off <- y*0 + .6 ## offset term to force curve through (0, .6)
## fit spline constrained through (0, .6)...
b <- gam(y ~ X - 1 + offset(off),paraPen=list(X=list(S)))
lines(x,predict(b))
## compare to unconstrained fit...
b.u <- gam(y ~ s(x,k=9),data=dat,knots=knots)
lines(x,predict(b.u))
*My question*: how can I extend your example to more than one smooth
terms, and several smooth interactions? In the context of a model where
E[ResponseRatio] must equal 0 when intensity equals zero?
I am not sure that I understand exactly what is going on when you call
smoothCon to specify the `sm`. I've written my own penalized spline
code from scratch, but it is far less sophisticated than mgcv, and is
basically a ridge regression where I optimize to get a lambda after
specifying a model with a lot of knots. mgcv clearly has a lot more
going on, and is far preferable to my rudimentary code for its handling
of tensors and random effects.
(Also, for prediction, how can I do "by=dum" AND "by=intensity" at the
same time?)
Many thanks,
Andrew
PS: I am aware that interacting my model with an intensity variable
makes my model quite heteroskedastic. I am thinking of using a cluster
wild bootstrap to construct confidence intervals. If a better way
forward immediately comes to your mind -- especially if its
computationally cheaper -- I'd greatly appreciate if you could share it.
On 04/17/2013 02:16 AM, Simon Wood wrote:
> hi Andrew.
>
> gam does suppress the intercept, it's just that this doesn't force the
> smooth through the intercept in the way that you would like. Basically
> for the parameteric component of the model '-1' behaves exactly like
> it does in 'lm' (it's using the same code). The smooths are 'added on'
> to the parametric component of the model, with sum to zero constraints
> to force identifiability.
>
> There is a solution to forcing a spline through a particular point at
> http://r.789695.n4.nabble.com/Use-pcls-in-quot-mgcv-quot-package-to-achieve-constrained-cubic-spline-td4660966.html
>
>
> (i.e. the R help thread "Re: [R] Use pcls in "mgcv" package to achieve
> constrained cubic spline")
>
> best,
> Simon
>
> On 16/04/13 22:36, Andrew Crane-Droesch wrote:
>>> Dear List,
>>>
>>> I've just tried to specify a GAM without an intercept -- I've got one
>>> of the (rare) cases where it is appropriate for E(y) -> 0 as X ->0.
>>> Naively running a GAM with the "-1" appended to the formula and the
>>> calling "predict.gam", I see that the model isn't behaving as expected.
>>>
>>> I don't understand why this would be. Google turns up this old R help
>>> thread:
>>> http://r.789695.n4.nabble.com/GAM-without-intercept-td4645786.html
>>>
>>> Simon writes:
>>>
>>> *Smooth terms are constrained to sum to zero over the covariate
>>> values. **
>>> **This is an identifiability constraint designed to avoid
>>> confounding with **
>>> **the intercept (particularly important if you have more than one
>>> smooth). *
>>> If you remove the intercept from you model altogether (m2) then
>>> the
>>> smooth will still sum to zero over the covariate values, which in
>>> your
>>> case will mean that the smooth is quite a long way from the data.
>>> When
>>> you include the intercept (m1) then the intercept is effectively
>>> shifting the constrained curve up towards the data, and you get a
>>> nice fit.
>>>
>>> Why? I haven't read Simon's book in great detail, though I have read
>>> Ruppert et al.'s Semiparametric Regression. I don't see a reason why
>>> a penalized spline model shouldn't equal the intercept (or zero) when
>>> all of the regressors equals zero.
>>>
>>> Is anyone able to help with a bit of intuition? Or relevant passages
>>> from a good description of why this would be the case?
>>>
>>> Furthermore, why does the "-1" formula specification work if it
>>> doesn't work "as intended" by for example lm?
>>>
>>> Many thanks,
>>> Andrew
>>>
>>>
>>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
