R / mgcv: Dlaczego produkty tensorowe te () i ti () wytwarzają różne powierzchnie?

mgcvOpakowanie Rposiada dwie funkcje montowania interakcji produktów napinacz: te()i ti(). Rozumiem podstawowy podział pracy między nimi (dopasowanie interakcji nieliniowej vs. rozkładanie tej interakcji na główne efekty i interakcję). To, czego nie rozumiem, to dlaczego te(x1, x2)i ti(x1) + ti(x2) + ti(x1, x2)może powodować (nieznacznie) różne wyniki.

MWE (dostosowany z ?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))

Różnice nie są bardzo duże w tym przykładzie, ale zastanawiam się, dlaczego w ogóle powinny być takie różnice.

Informacje o sesji:

 > 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)

r gam mgcv conditional-probability mixed-model references bayesian estimation conditional-probability machine-learning optimization gradient-descent r hypothesis-testing wilcoxon-mann-whitney time-series bayesian inference change-point time-series anova repeated-measures statistical-significance bayesian contingency-tables regression prediction quantiles classification auc k-means scikit-learn regression spatial circular-statistics t-test effect-size cohens-d r cross-validation feature-selection caret machine-learning modeling python optimization frequentist correlation sample-size normalization group-differences heteroscedasticity independence generalized-least-squares lme4-nlme references mcmc metropolis-hastings optimization r logistic feature-selection separation clustering k-means normal-distribution gaussian-mixture kullback-leibler java spark-mllib data-visualization categorical-data barplot hypothesis-testing statistical-significance chi-squared type-i-and-ii-errors pca scikit-learn conditional-expectation statistical-significance meta-analysis intuition r time-series multivariate-analysis garch machine-learning classification data-mining missing-data cart regression cross-validation matrix-decomposition categorical-data repeated-measures chi-squared assumptions contingency-tables prediction binary-data trend test-for-trend matrix-inverse anova categorical-data regression-coefficients standard-error r distributions exponential interarrival-time copula log-likelihood time-series forecasting prediction-interval mean standard-error meta-analysis meta-regression network-meta-analysis systematic-review normal-distribution multiple-regression generalized-linear-model poisson-distribution poisson-regression r sas cohens-kappa

— jvh_ch
źródło

Poważnie ludzie !? Podczas gdy implementacja jest wyraźnie specyficzna dla mgcv (nie jestem świadomy żadnego innego gotowego oprogramowania dla GAM, które pozwala na podobny do ANOVA rozkład dwuwymiarowych wygładzeń), problem i odpowiedź są wyraźnie statystyczne; dopasowywane modele nie są takie same pod maską ze względu na dodatkowe macierze kar, które powstają podczas dekompozycji warunków krańcowych z komponentu „interakcji”. Nie dotyczy to mgcv.

— Gavin Simpson

@GavinSimpson Zadałem Meta pytanie na temat aktualności tego pytania lub inaczej

— Silverfish,

Są to pozornie ten sam model, ale w praktyce przy montażu występują pewne subtelne różnice. Jedną ważną różnicą jest to, że model z ti()terminami szacuje więcej parametrów gładkości w porównaniu z te()modelem:

> 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

a to dlatego, że istnieje więcej macierzy kar związanych z tymi dwoma modelami; w ti()modelu mamy jeden na „termin” w porównaniu z tylko dwoma w te()modelu, jeden na margines gładki.

Widzę modele ti()używane jako decydujące, czy chcę czy . Nie mogę porównać tych modeli, jeśli używam terminów, więc używam . Kiedy już ustalę, czy potrzebuję , mogę zmodernizować model, jeśli tego potrzebuję, lub osobno dla każdego efektu krańcowego, jeśli nie potrzebuję . $\hat{y} = \beta_0 + s(x, y)$ $\hat{y} = \beta_0 + s(x) + s(y)$ te()ti() $s(x,y)$ te()s() $s(x,y)$

Zauważ, że możesz zbliżyć modele nieco bliżej siebie, dopasowując je za pomocą method = "ML"(lub "REML", ale nie powinieneś porównywać „ustalonych” efektów, "REML"chyba że wszystkie warunki są w pełni ukarane, co domyślnie nie jest, ale byłoby powiedziane z select = TRUE).

— Gavin Simpson
źródło