R/mgcv: Why do te() and ti() tensor products produce different surfaces? The mgcv package for R has two functions for fitting tensor product interactions: te() and ti(). I understand the basic division of labour between the two (fitting a non-linear interaction vs. decomposing this interaction into main effects and an interaction). What I don't understand is why te(x1, x2) and ti(x1) + ti(x2) + ti(x1, x2) may produce (slightly) different results.
MWE (adapted from ?ti):
require(mgcv)
test1 <- function(x,z,sx=0.3,sz=0.4) { 
  x <- x*20
 (pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+
             0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))
}
n <- 500

x <- runif(n)/20;z <- runif(n);
xs <- seq(0,1,length=30)/20;zs <- seq(0,1,length=30)
pr <- data.frame(x=rep(xs,30),z=rep(zs,rep(30,30)))
truth <- matrix(test1(pr$x,pr$z),30,30)
f <- test1(x,z)
y <- f + rnorm(n)*0.2

par(mfrow = c(2,2))

# Model with te()
b2 <- gam(y~te(x,z))
vis.gam(b2, plot.type = "contour", color = "terrain", main = "tensor product")

# Model with ti(a) + ti(b) + ti(a,b)
b3 <- gam(y~ ti(x) + ti(z) + ti(x,z))
vis.gam(b3, plot.type = "contour", color = "terrain", main = "tensor anova")

# Scatterplot of prediction b2/b3
plot(predict(b2), predict(b3))

The differences aren't very large in this example, but I'm just wondering why there should be differences at all.
Session info:
 > devtools::session_info("mgcv")
 Session info
 -----------------------------------------------------------------------------------
 setting  value                       
 version  R version 3.3.1 (2016-06-21)
 system   x86_64, linux-gnu           
 ui       RStudio (0.99.491)          
 language en_US                       
 collate  en_US.UTF-8                 
 tz       <NA>                        
 date     2016-09-13                  

 Packages      ---------------------------------------------------------------------------------------
 package * version date       source        
 lattice   0.20-33 2015-07-14 CRAN (R 3.2.1)
 Matrix    1.2-6   2016-05-02 CRAN (R 3.3.0)
 mgcv    * 1.8-12  2016-03-03 CRAN (R 3.2.3)
 nlme    * 3.1-128 2016-05-10 CRAN (R 3.3.1)

 A: These are superficially the same model but in practice when fitting there are some subtle differences. One important difference is that the model with ti() terms is estimating more smoothness parameters compared with the te() model:
> b2$sp
te(x,z)1 te(x,z)2 
3.479997 5.884272 
> b3$sp
    ti(x)     ti(z)  ti(x,z)1  ti(x,z)2 
 8.168742 60.456559  2.370604  2.761823

and this is because there are more penalty matrices associated with the two models; in the ti() model we have one per "term" compared with just two in the te() model, one per marginal smooth.
I see models with ti() as being used to decide whether I want $\hat{y} = \beta_0 + s(x, y)$ or $\hat{y} = \beta_0 + s(x) + s(y)$. I can't compare those models if I use te() terms so I use ti(). Once I've determined if I need $s(x,y)$ I can refit the model with te() if I need it or with separate s() for each marginal effect if I don't need $s(x,y)$.
Note that you can get the models somewhat closer to each other by fitting using method = "ML" (or "REML", but you shouldn't be comparing "fixed" effects with "REML" unless all the terms are fully penalized, which by default they aren't, but would be say with select = TRUE).
