R / mgcv: Por que os produtos tensores te () e ti () produzem superfícies diferentes?

11

O mgcvpacote para Rpossui duas funções para ajustar as interações do produto tensorial: te()e ti(). Entendo a divisão básica do trabalho entre os dois (ajustando uma interação não linear versus decompondo essa interação em efeitos principais e uma interação). O que não entendo é o porquê te(x1, x2)e ti(x1) + ti(x2) + ti(x1, x2)pode produzir resultados (ligeiramente) diferentes.

MWE (adaptado 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))

As diferenças não são muito grandes neste exemplo, mas estou me perguntando por que deveria haver diferenças.

Informações da sessão:

 > 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
fonte
4
Sério gente !? Enquanto a implementação de forma clara uma coisa específica-mgcv (eu sou desconhecem qualquer outro software off-the-shelf para os modelos MAG que permite que este ANOVA-like decomposição de suaviza bivariadas), o problema e resposta são claramente estatística; os modelos que estão sendo ajustados não são os mesmos por baixo das capas devido às matrizes de penalidade extra que surgem ao decompor os termos marginais do componente "interação". Isto não é específico para mgcv.
Gavin Simpson
11
@GavinSimpson Criei uma pergunta no Meta sobre a actualidade em-, ou de outro modo, desta questão
Silverfish

Respostas:

15

Estes são superficialmente o mesmo modelo, mas na prática, quando ajustados, existem algumas diferenças sutis. Uma diferença importante é que o modelo com ti()termos está estimando mais parâmetros de suavidade em comparação com o te()modelo:

> 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

e isso ocorre porque existem mais matrizes de penalidade associadas aos dois modelos; no ti()modelo, temos um por "termo" em comparação com apenas dois no te()modelo, um por margem marginal.

ti()y^=β0 0+s(x,y)y^=β0 0+s(x)+s(y)te()ti()s(x,y)te()s()s(x,y)

Observe que você pode aproximar os modelos um do outro ajustando usando method = "ML"(ou "REML", mas você não deve comparar efeitos "fixos" com, a "REML"menos que todos os termos sejam totalmente penalizados, o que, por padrão, eles não são, mas diria com select = TRUE).

Gavin Simpson
fonte