R / mgcv: Pourquoi les produits tenseurs te () et ti () produisent-ils des surfaces différentes?

11

Le mgcvpackage pour Ra deux fonctions pour ajuster les interactions des produits tensoriels: te()et ti(). Je comprends la division de base du travail entre les deux (ajustement d'une interaction non linéaire vs décomposition de cette interaction en effets principaux et interaction). Ce que je ne comprends pas, c'est pourquoi te(x1, x2)et ti(x1) + ti(x2) + ti(x1, x2)peut produire des résultats (légèrement) différents.

MWE (adapté de ?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))

Les différences ne sont pas très importantes dans cet exemple, mais je me demande simplement pourquoi il devrait y avoir des différences.

Informations sur la session:

 > 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
la source
4
Sérieusement les gens!? Alors que l'implémentation est clairement une chose spécifique au mgcv (je ne connais aucun autre logiciel standard pour les GAM qui permette cette décomposition de lissage bivariée de type ANOVA), le problème et la réponse sont clairement statistiques; les modèles ajustés ne sont pas les mêmes sous le capot en raison des matrices de pénalités supplémentaires qui surviennent lors de la décomposition des termes marginaux de la composante "interaction". Ce n'est pas spécifique au mgcv.
Gavin Simpson
1
@GavinSimpson J'ai posé une question sur Meta au sujet de l'actualité, ou autrement, de cette question
Silverfish

Réponses:

15

Ce sont superficiellement le même modèle, mais dans la pratique lors du montage, il y a quelques différences subtiles. Une différence importante est que le modèle avec ti()termes estime plus de paramètres de lissage par rapport au te()modèle:

> 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

et c'est parce qu'il y a plus de matrices de pénalités associées aux deux modèles; dans le ti()modèle, nous en avons un par «terme», contre seulement deux dans le te()modèle, un par lisse marginale.

Je vois des modèles avec ti()comme étant utilisés pour décider si je veux ou . Je ne peux pas comparer ces modèles si j'utilise des termes donc j'utilise . Une fois que j'ai déterminé si j'ai besoin de je peux refit le modèle avec si j'en ai besoin ou avec séparé pour chaque effet marginal si je n'ai pas besoin de .y^=β0+s(x,y)y^=β0+s(x)+s(y)te()ti()s(x,y)te()s()s(x,y)

Notez que vous pouvez rapprocher les modèles les uns des autres en ajustant à l'aide de method = "ML"(ou "REML", mais vous ne devriez pas comparer les effets "fixes" avec à "REML"moins que tous les termes ne soient entièrement pénalisés, ce qui par défaut ne l'est pas, mais dirait avec select = TRUE).

Gavin Simpson
la source