R Under development (unstable) (2024-06-20 r86796 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> library(MuMIn)
> suppressPackageStartupMessages(library(mgcv))
> RNGkind("Mersenne")
> 
> 
> set.seed(0) ## simulate some data...
> dat <- gamSim(1, n = 200, dist = "binary", scale = 2)
Gu & Wahba 4 term additive model
> 
> 
> gamsmoothwrap <-
+ function(formula, k1 = NA, ...) { 
+     cl <- origCall <- match.call()
+     cl[[1]] <- as.name("gam")
+     cl$formula <- exprApply(formula, "s", function(e, k, x) {
+         i <- which(e[[2]] == x)[1]
+         if(!is.na(i) && !is.na(k[i])) e[["k"]] <- k[i]
+         e
+     }, k = c(k1), x = c("x0"))
+     cl$k1 <- NULL
+ 	fit <- eval(cl, parent.frame())
+     fit$call <- origCall # replace the stored call
+     fit
+ }
> 
> glo <- gamsmoothwrap(y ~ s(x0, fx = TRUE), k1 = 3, data = dat, family = binomial)
> 
> 
> print(ms <- model.sel(models <- lapply(dredge(glo, varying = list(k1 = 3:16), fixed = TRUE,
+ 	evaluate = FALSE), eval)))
Fixed terms are "s(x0, fx = TRUE, k = 3)" and "(Intercept)"
Model selection table 
   (Int) s(x0,T,3) s(x0,T,4) s(x0,T,5) s(x0,T,6) s(x0,T,7) s(x0,T,8) s(x0,T,9)
1  1.293         +                                                            
2  1.321                   +                                                  
3  1.321                             +                                        
4  1.329                                       +                              
5  1.348                                                 +                    
13 2.126                                                                      
10 1.554                                                                      
6  1.360                                                           +          
14 1.964                                                                      
11 1.575                                                                      
7  1.361                                                                     +
12 1.679                                                                      
9  1.424                                                                      
8  1.360                                                                      
   s(x0,T,10) s(x0,T,11) s(x0,T,12) s(x0,T,13) s(x0,T,14) s(x0,T,15) s(x0,T,16)
1                                                                              
2                                                                              
3                                                                              
4                                                                              
5                                                                              
13                                                                 +           
10                                +                                            
6                                                                              
14                                                                            +
11                                           +                                 
7                                                                              
12                                                      +                      
9                      +                                                       
8           +                                                                  
   k1 df   logLik  AICc delta weight
1   3  2 -102.341 210.8  0.00  0.373
2   4  4 -101.430 211.1  0.26  0.327
3   5  5 -101.391 213.1  2.29  0.119
4   6  6 -101.216 214.9  4.06  0.049
5   7  6 -100.616 215.8  5.01  0.030
13 15 15  -92.011 216.6  5.82  0.020
10 12 12  -95.543 216.8  5.95  0.019
6   8  8 -100.151 217.1  6.25  0.016
14 16 16  -91.066 217.1  6.30  0.016
11 13 13  -94.846 217.6  6.84  0.012
7   9  8 -100.142 219.2  8.43  0.006
12 14 13  -94.490 219.3  8.45  0.005
9  11 11  -98.147 219.7  8.89  0.004
8  10 10 -100.136 221.4 10.63  0.002
Models ranked by AICc(x) 
> 
> 
> # should throw a warning
> stopifnot(isTRUE(tryCatch(model.avg(models), warning = function(e) TRUE)))
> 
> #_______________________________________________________________________________
> 
> set.seed(0) ## simulate some data...
> dat <- gamSim(1, n = 400, dist = "binary", scale = 2)
Gu & Wahba 4 term additive model
> 
> testSmoothKConsistency <- MuMIn:::testSmoothKConsistency
> 
> 
> glm1 <- glm(y ~ x1 * x2, data = dat, family = binomial)
> gam53 <- gam(y ~ s(x0, k = 5) + s(x1, k = 3), data = dat, family = binomial)
> gam37 <- update(gam53, . ~ s(x0, k = 3) + s(x1, k = 7))
> gam3e44 <- update(gam53, . ~ s(x0, k = 3) + te(x2, x1, k = 4))
> gam3e54 <- update(gam53, . ~ s(x0, k = 3) + te(x2, x1, k = c(5, 4)))
> gam3i54 <- update(gam53, . ~ s(x0, k = 3) + ti(x2, x1, k = c(5, 4)))
> 
> gam3e54 <- update(gam53, . ~ s(x0, k = 3) + te(x2, x1, k = c(5, 4)))
> gam3i34 <- update(gam53, . ~ s(x0, k = 3) + ti(x2, x1, k = c(3, 4)))
> 
> 
> ms <- model.sel(gam53, gam37, gam3e44, gam3e54, gam3i54, glm1)
> testSmoothKConsistency(ms)
Warning message:
In testSmoothKConsistency(ms) :
  smooth term dimensions differ between models for variables 'x0', 'x1' and 'x1,x2'. Related coefficients are incomparable.
> testSmoothKConsistency(list(gam3e54, gam3i34))
> 
> 
> 
> #_______________________________________________________________________________
> 
> 
> 
> proc.time()
   user  system elapsed 
   1.54    0.12    1.67