El mgcv
paquete R
tiene dos funciones para ajustar las interacciones del producto tensorial: te()
y ti()
. Entiendo la división básica del trabajo entre los dos (ajustar una interacción no lineal versus descomponer esta interacción en efectos principales y una interacción). Lo que no entiendo es por qué te(x1, x2)
y ti(x1) + ti(x2) + ti(x1, x2)
puede producir resultados (ligeramente) 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))
Las diferencias no son muy grandes en este ejemplo, pero me pregunto por qué debería haber diferencias.
Información de la sesión:
> 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
fuente
fuente
Respuestas:
Estos son superficialmente el mismo modelo, pero en la práctica, cuando se ajusta, existen algunas diferencias sutiles. Una diferencia importante es que el modelo con
ti()
términos está estimando más parámetros de suavidad en comparación con elte()
modelo:y esto se debe a que hay más matrices de penalización asociadas con los dos modelos; en el
ti()
modelo tenemos uno por "término" en comparación con solo dos en elte()
modelo, uno por marginal liso.Veo modelosy^=β0+s(x,y) y^=β0+s(x)+s(y) s(x,y) s(x,y)
ti()
que se usan para decidir si quiero o . No puedo comparar esos modelos si uso términos, así que los uso . Una vez que he determinado si necesito puedo reajustar el modelo con si lo necesito o con separado para cada efecto marginal si no necesito .te()
ti()
te()
s()
Tenga en cuenta que puede acercar los modelos un poco entre sí mediante el uso de
method = "ML"
(o"REML"
, pero no debería comparar los efectos "fijos" a"REML"
menos que todos los términos estén totalmente penalizados, lo que por defecto no lo son, pero se diría conselect = TRUE
).fuente