I have some raw data on which I compute percent changes, and I calculate rolling averages of two different lengths on the percent changes. I want to use a tensor interaction of the two rolling averages as an input to a GAM. The percent changes have some bias (edit/update: if you remove the bias from the response this problem goes away, as expected -- but I don't think this makes the question any less valid), and I have fit models with and without intercept terms to the data (Note: I recognize that these data are meaningless and it's probably not a good idea to force a zero intercept in this case, but as you'll see below, this evokes some behavior that I do not understand).
What I find confusing is that models with a univariate smooth term have just about exactly the same coefficient values whether or not the intercept is included, but models with a two-dimensional tensor interaction yield very different coefficients with and without intercept. For example:
Single univariate smooth (case 1) --
n <- 10000
mean <- 0
sd <- 1
raw.data <- (1:n)/1000 + rnorm(n=n, mean=mean, sd=sd) + 100
y <- diff(raw.data) / raw.data[-n] +.003
x1 <- rollmean(x=y, k=n/100, fill=NA, align='right')
x2 <- rollmean(x=y, k=n/500, fill=NA, align='right')
fit.s.1 <- gam(formula=as.formula("y~s(x1, bs='cr', k=5)"))
fit.s.2 <- gam(formula=as.formula("y~s(x1, bs='cr', k=5)-1"))
data.frame(fit.s.1$coefficients, c(NA, fit.s.2$coefficients))
Output:
with.intercept without.intercept
(Intercept) 0.003098624 NA
s(x1).1 0.013182549 0.013168485
s(x1).2 -0.001523374 -0.001526166
s(x1).3 0.022703465 0.022685008
s(x1).4 0.027567431 0.027645856
Tensor interaction (case 2) --
fit.ti.1 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5)"))
fit.ti.2 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5)-1"))
data.frame(fit.ti.1$coefficients, c(NA, fit.ti.2$coefficients))
Output:
with.intercept without.intercept
(Intercept) 0.00317231024 NA
ti(x1,x2).1 -0.00163077279 -0.00601973430
ti(x1,x2).2 -0.00007483268 0.00191561176
ti(x1,x2).3 -0.01574417645 -0.01731746447
ti(x1,x2).4 0.00279508904 -0.03562415825
ti(x1,x2).5 -0.00579243943 -0.00447007353
ti(x1,x2).6 0.00001615630 -0.00029669803
ti(x1,x2).7 -0.00795505768 -0.00706923863
ti(x1,x2).8 -0.00189071125 -0.00009736479
ti(x1,x2).9 -0.00915478215 -0.01452895684
ti(x1,x2).10 -0.00386687694 -0.00206597434
ti(x1,x2).11 -0.02346438804 -0.02221899031
ti(x1,x2).12 0.02186800209 -0.01123017941
ti(x1,x2).13 -0.01277916564 0.01313048335
ti(x1,x2).14 -0.00707688210 -0.00096508269
ti(x1,x2).15 -0.00427769077 0.02633769973
ti(x1,x2).16 0.03576640248 0.03967865100
I thought this might have something to do with the centering constraints placed on ti() marginals, so I disabled them and reran the model. This makes the coefficients much closer but not nearly identical like in the case of the univariate smooth.
Tensor interaction with centering constraints disabled (case 3) --
fit.ti.3 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5, mc=c(F,F))"))
fit.ti.4 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5, mc=c(F,F))-1"))
data.frame(format(fit.ti.3$coefficients, scientific=F), format(c(NA, fit.ti.4$coefficients), scientific=F))
Output:
with.intercept without.intercept
(Intercept) 0.00332043150 NA
ti(x1,x2).1 -0.03286230398 -0.02954187248
ti(x1,x2).2 -0.02003935751 -0.01671892601
ti(x1,x2).3 -0.01707730042 -0.01375686892
ti(x1,x2).4 -0.01400703252 -0.01068660102
ti(x1,x2).5 -0.00139113625 0.00192929525
ti(x1,x2).6 -0.01945937590 -0.01613894440
ti(x1,x2).7 -0.00627836343 -0.00295793193
ti(x1,x2).8 -0.00323360264 0.00008682886
ti(x1,x2).9 -0.00007762914 0.00324280236
ti(x1,x2).10 0.01289038386 0.01621081536
ti(x1,x2).11 -0.01660715724 -0.01328672574
ti(x1,x2).12 -0.00334994895 -0.00002951745
ti(x1,x2).13 -0.00028759225 0.00303283924
ti(x1,x2).14 0.00288662690 0.00620705840
ti(x1,x2).15 0.01592956097 0.01924999246
ti(x1,x2).16 -0.01376515649 -0.01044472499
ti(x1,x2).17 -0.00043203373 0.00288839777
ti(x1,x2).18 0.00264785857 0.00596829006
ti(x1,x2).19 0.00584025237 0.00916068387
ti(x1,x2).20 0.01895784225 0.02227827375
ti(x1,x2).21 0.00168433249 0.00500476399
ti(x1,x2).22 0.01543005061 0.01875048211
ti(x1,x2).23 0.01860523453 0.02192566603
ti(x1,x2).24 0.02189639869 0.02521683019
ti(x1,x2).25 0.03541988069 0.03874031219
Now, if I add a univariate smooth to the tensor interaction term (but use the ti() function to construct it), the interaction term coefficients are all exactly the same but the coefficients of the univariate smooth are very different.
Case 4 --
fit.ti.ti.1 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5, np=F, mc=c(F,T))+ti(x1, bs='cr', k=5, np=T, mc=F)"))
fit.ti.ti.2 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5, np=F, mc=c(F,T))+ti(x1, bs='cr', k=5, np=T, mc=F)-1"))
Output:
with.intercept without.intercept
(Intercept) 2.362528e-03 NA
ti(x1,x2).1 7.183172e-03 7.183172e-03
ti(x1,x2).2 -2.343559e-05 -2.343559e-05
ti(x1,x2).3 1.253984e-02 1.253984e-02
ti(x1,x2).4 1.488033e-02 1.488033e-02
ti(x1,x2).5 7.941893e-03 7.941893e-03
ti(x1,x2).6 -2.590619e-05 -2.590619e-05
ti(x1,x2).7 1.386434e-02 1.386434e-02
ti(x1,x2).8 1.645201e-02 1.645201e-02
ti(x1,x2).9 8.103348e-03 8.103348e-03
ti(x1,x2).10 -2.643482e-05 -2.643482e-05
ti(x1,x2).11 1.414620e-02 1.414620e-02
ti(x1,x2).12 1.678646e-02 1.678646e-02
ti(x1,x2).13 8.264221e-03 8.264221e-03
ti(x1,x2).14 -2.695747e-05 -2.695747e-05
ti(x1,x2).15 1.442704e-02 1.442704e-02
ti(x1,x2).16 1.711969e-02 1.711969e-02
ti(x1,x2).17 9.138680e-03 9.138680e-03
ti(x1,x2).18 -2.980949e-05 -2.980949e-05
ti(x1,x2).19 1.595358e-02 1.595358e-02
ti(x1,x2).20 1.893113e-02 1.893113e-02
ti(x1).1 -1.691305e-02 -1.455052e-02
ti(x1).2 -2.274243e-03 8.828505e-05
ti(x1).3 7.960668e-04 3.158595e-03
ti(x1).4 3.733641e-03 6.096169e-03
ti(x1).5 1.702011e-02 1.938264e-02
Interestingly, a full tensor model behaves similarly to a univariate smooth model in that its coefficients are almost identical whether or not the intercept term is present:
Case 5 --
fit.te.1 <- gam(formula=as.formula("y~te(x1, x2, bs='cr', k=5)"))
fit.te.2 <- gam(formula=as.formula("y~te(x1, x2, bs='cr', k=5)-1"))
data.frame(fit.te.1$coefficients, c(NA, fit.te.2$coefficients))
Output:
with.intercept without.intercept
(Intercept) 0.0030986243 NA
te(x1,x2).1 -0.0184806514 -0.0184806500
te(x1,x2).2 -0.0172327917 -0.0172327903
te(x1,x2).3 -0.0136081242 -0.0136081233
te(x1,x2).4 -0.0011718971 -0.0011719087
te(x1,x2).5 -0.0176757491 -0.0176757526
te(x1,x2).6 0.0147984258 0.0147984322
te(x1,x2).7 -0.0054231322 -0.0054231330
te(x1,x2).8 0.0076246608 0.0076246630
te(x1,x2).9 0.0133582270 0.0133582238
te(x1,x2).10 -0.0168794331 -0.0168794363
te(x1,x2).11 -0.0064646710 -0.0064646717
te(x1,x2).12 0.0049684688 0.0049684702
te(x1,x2).13 -0.0001598714 -0.0001598725
te(x1,x2).14 0.0156425216 0.0156425199
te(x1,x2).15 -0.0133135304 -0.0133135324
te(x1,x2).16 0.0075429791 0.0075429814
te(x1,x2).17 0.0008931768 0.0008931761
te(x1,x2).18 0.0270142207 0.0270142264
te(x1,x2).19 0.0208405029 0.0208405035
te(x1,x2).20 0.0019047894 0.0019047918
te(x1,x2).21 0.0158133845 0.0158133838
te(x1,x2).22 0.0185985651 0.0185985650
te(x1,x2).23 0.0231006405 0.0231006421
te(x1,x2).24 0.0357532266 0.0357532341
I thought that perhaps the models which did not seem to "care" about the intercept were set up such that some linear combination of their bases was equal (at least approximately) to a constant, so I tested this by running
X <- model.matrix(fit.te.2)
X %*% solve(t(X) %*% X) %*% t(X) %*% rep(mean(y), nrow(X))
with the expectation that the result would be mean(y) repeated for all observations. However, this is not what happens (note: it does happen in case 4, in which there is a tensor interaction plus a univariate smooth, but no other cases).
So, I apologize for the extremely long post, but I would really appreciate some insight into why this behavior occurs. Thanks in advance for your help and guidance.
The following is the entirety of the script I am working with, for your reference:
library(zoo)
n <- 10000
mean <- 0
sd <- 1
raw.data <- (1:n)/1000 + rnorm(n=n, mean=mean, sd=sd) + 100
y <- diff(raw.data) / raw.data[-n] +.003
x1 <- rollmean(x=y, k=n/100, fill=NA, align='right')
x2 <- rollmean(x=y, k=n/500, fill=NA, align='right')
fit.s.1 <- gam(formula=as.formula("y~s(x1, bs='cr', k=5)"))
fit.s.2 <- gam(formula=as.formula("y~s(x1, bs='cr', k=5)-1"))
data.frame(fit.s.1$coefficients, c(NA, fit.s.2$coefficients))
X <- model.matrix(fit.s.2)
X %*% solve(t(X) %*% X) %*% t(X) %*% rep(mean(y), nrow(X))
fit.s.3 <- gam(formula=as.formula("y~s(x1, x2, bs='tp')"))
fit.s.4 <- gam(formula=as.formula("y~s(x1, x2, bs='tp')-1"))
data.frame(fit.s.3$coefficients, c(NA, fit.s.4$coefficients))
X <- model.matrix(fit.s.4)
X %*% solve(t(X) %*% X) %*% t(X) %*% rep(mean(y), nrow(X))
fit.ti.1 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5)"))
fit.ti.2 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5)-1"))
data.frame(with.intercept=format(fit.ti.1$coefficients, scientific=F), without.intercept=format(c(NA, fit.ti.2$coefficients), scientific=F))
X <- model.matrix(fit.ti.2)
X %*% solve(t(X) %*% X) %*% t(X) %*% rep(mean(y), nrow(X))
fit.ti.3 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5, mc=c(F,F))"))
fit.ti.4 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5, mc=c(F,F))-1"))
data.frame(with.intercept=format(fit.ti.3$coefficients, scientific=F), without.intercept=format(c(NA, fit.ti.4$coefficients), scientific=F))
X <- model.matrix(fit.ti.4)
X %*% solve(t(X) %*% X) %*% t(X) %*% rep(mean(y), nrow(X))
fit.te.1 <- gam(formula=as.formula("y~te(x1, x2, bs='cr', k=5)"))
fit.te.2 <- gam(formula=as.formula("y~te(x1, x2, bs='cr', k=5)-1"))
data.frame(fit.te.1$coefficients, c(NA, fit.te.2$coefficients))
X <- model.matrix(fit.te.2)
X %*% solve(t(X) %*% X) %*% t(X) %*% rep(mean(y), nrow(X))
fit.ti.ti.1 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5, np=F, mc=c(F,T))+ti(x1, bs='cr', k=5, np=T, mc=F)"))
fit.ti.ti.2 <- gam(formula=as.formula("y~ti(x1, x2, bs='cr', k=5, np=F, mc=c(F,T))+ti(x1, bs='cr', k=5, np=T, mc=F)-1"))
data.frame(fit.ti.ti.1$coefficients, c(NA, fit.ti.ti.2$coefficients))
X <- model.matrix(fit.ti.ti.2)
X %*% solve(t(X) %*% X) %*% t(X) %*% rep(mean(y), nrow(X))
ti()smooths??tisuggests that they assume that you have done this. I would expect your models to haveti(x1) + ti(x2) + ti(x1, x2). If you don't need this functional ANOVA-like decomposition, just usete(x1, x2, ....)where you haveti(x1, x2, ....). Also, it would be nicer to read if you just wrote the formula outside"and didn't coerce it to aformulaobject; R knows what to do if you just doy ~ ti(x1) + ti(x2) + ti(x1, x2)without quoting and redundant coercion. $\endgroup$y ~ ti(x1) + ti(x2) + ti(x1, x2)andy ~ ti(x1) + ti(x2) + ti(x1, x2)-1? This behaves just like theti()interaction by itself (the coefficients are very different with and without the intercept). On another note, I had been under the impression thaty ~ ti(x1) + ti(x2) + ti(x1, x2)was equivalent toy ~ te(x1, x2), but this does not seem to be the case. $\endgroup$