search

refunders / refund

Regression with functional data
39 stars 23 forks source link

pffr: ti and discrete = TRUE don't play well together #70

Closed fabian-s closed 7 years ago

fabian-s commented 8 years ago 
library(refund)
library(mgcv)

set.seed(122212)
dat <- pffrSim(scenario = c("int", "smoo", "lin"), n=200, SNR=1)
m_discrete <- pffr(Y ~ s(xsmoo) + s(xlin), data = dat, yind = attr(dat, "yindex"),
  method="fREML", algorithm = "bam", discrete  = TRUE, chunk.size = 1000)
m_nodiscrete <- pffr(Y ~ s(xsmoo) + s(xlin), data = dat, yind = attr(dat, "yindex"),
  method="fREML", algorithm = "bam", discrete  = FALSE, chunk.size = 1000)
m_t2 <- pffr(Y ~ s(xsmoo) + s(xlin), data = dat, yind = attr(dat, "yindex"),
  method="fREML", algorithm = "bam", tensortype  = "t2", chunk.size = 1000)
# gam-call used internally by pffr:
m_discrete$call
# bam(formula = Y ~ s(x = yindex.vec, bs = "ps", k = 20, m = c(2,
#   1)) + ti(xsmoo, d = c(1, 1), yindex.vec, mc = c(TRUE, FALSE),
#     bs = c("tp", "ps"), m = c(2, 1), k = c(8, 5)) + ti(xlin,
#       d = c(1, 1), yindex.vec, mc = c(TRUE, FALSE), bs = c("tp",
#         "ps"), m = c(2, 1), k = c(8, 5)), data = pffrdata, method = "fREML",
#   chunk.size = 1000, discrete = TRUE)

# refit pffr model without custom constraints: use ti without mc-arg
#       and the prepped data from the original pffr-fit
m_nomc_discrete <- bam(formula = Y ~ s(x = yindex.vec, bs = "ps", k = 20, m = c(2,
     1)) + ti(xsmoo, d = c(1, 1), yindex.vec,
       bs = c("tp", "ps"), m = c(2, 1), k = c(8, 5)) + ti(xlin,
         d = c(1, 1), yindex.vec, bs = c("tp",
           "ps"), m = c(2, 1), k = c(8, 5)), data = m_discrete$model, method = "fREML",
  chunk.size = 1000, discrete = TRUE)

m_nomc_nodiscrete <- bam(formula = Y ~ s(x = yindex.vec, bs = "ps", k = 20, m = c(2,
  1)) + ti(xsmoo, d = c(1, 1), yindex.vec,
    bs = c("tp", "ps"), m = c(2, 1), k = c(8, 5)) + ti(xlin,
      d = c(1, 1), yindex.vec, bs = c("tp",
        "ps"), m = c(2, 1), k = c(8, 5)), data = m_discrete$model, method = "fREML",
  chunk.size = 1000, discrete = FALSE)

layout(matrix(1:15, 3, 5))
plot(m_discrete, pers=TRUE, pages = 0, ticktype="detailed")
plot(m_nomc_discrete, pers=TRUE, pages = 0, ticktype="detailed")
plot(m_nodiscrete, pers=TRUE, pages = 0, ticktype="detailed")
plot(m_nomc_nodiscrete, pers=TRUE, pages = 0, ticktype="detailed")
plot(m_t2, pers=TRUE, pages = 0, ticktype="detailed")

image

?!? why are axis orders different between plots ?!?)
true functions are (see pffrSim code):

so discrete makes a huge difference if mc is set (columns 1 & 2 vs 3 & 4; it should not!)
and t2 does weird things as well, at least with the intercept function...

all.equal(predict(m_discrete), predict(m_nodiscrete))
all.equal(predict(m_discrete), predict(m_t2))
all.equal(as.vector(predict(m_nomc_discrete)), as.vector(predict.gam(m_discrete)))
all.equal(as.vector(predict(m_nomc_discrete)), as.vector(predict(m_nomc_nodiscrete)))
# rel. diffs of .034 (!!), 0.035, .005 (!!), .033 (!!)
## --> discrete seems the culprit... ?!?
fabian-s commented 8 years ago

sessionInfo() R version 3.3.1 (2016-06-21) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 14.04.5 LTS

locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=de_DE.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=de_DE.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] mgcv_1.8-15 nlme_3.1-128 refund_0.1-16

loaded via a namespace (and not attached): [1] Rcpp_0.12.7 knitr_1.14 splines_3.3.1 magic_1.5-6
[5] MASS_7.3-45 munsell_0.4.3 colorspace_1.2-6 lattice_0.20-33 [9] grpreg_3.0-2 fda_2.4.4 minqa_1.2.4 plyr_1.8.4
[13] tools_3.3.1 parallel_3.3.1 grid_3.3.1 gtable_0.2.0
[17] htmltools_0.3.5 digest_0.6.10 gamm4_0.2-4 lme4_1.1-12
[21] Matrix_1.2-6 nloptr_1.0.4 ggplot2_2.1.0 RLRsim_3.1-2
[25] rmarkdown_0.9.6 pbs_1.1 scales_0.4.0 boot_1.3-18

fabian-s commented 8 years ago

Understood so far:

fabian-s commented 7 years ago

fixed in mgcv 1.8-20