Issue of reproducibility with gam and lm.wfit in different versions of R

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

Issue of reproducibility with gam and lm.wfit in different versions of R

Sebastien Bihorel
Dear R users,

I recently stumbled upon problems of reproducibility while running GAM analyses in different R and gam package versions. In the example below, a small dataset is created in which the y and x1 variables are 100% correlated. The intents of this example were primarily for regression testing and, secondarily, to evaluate how the gam algorithm behaves under extreme/limit conditions.

I ran this little snippet in 5 different environments and got 100% consistent results until I switched to R 3.3.2.
* Comparing results from environments 1, 2, and 3 shows that changing the version of the gam package did not change the results under R 3.3.0.
* Comparing results from environments 3 and 4 shows that changing the version of R altered the values of the AIC and the output of the step.gam call (changed to a NULL object)
* Comparing results from environments 4 and 5 shows that reverting to an older version of the gam package in R 3.3.2 still produced altered AIC values and the NULL output from step.gam call

Further investigations into these differences seem to show that the lm.wfit call in the gam.fit function (called from within gam and step.gam) may result in different values in R 3.3.0 and 3.3.2.

So my questions are 2-fold:
1- Would you have any information about why lm.wfit would produce different outcomes in R 3.3.0 and 3.3.2? The source code does not appear significantly different.
2- Is it expected that the step.gam function returns a NULL object in R 3.3.2 while a valid model has been identified? Looking at the source of step.gam, the line 157 (if(is.null(form.list)) break) seems to be the reason the function breaks out and returns a NULL value.

I thank you in advance for your time.

Sebastien Bihorel






library(gam)

dat <- data.frame(
y = c(57,57,98,83,122,69,108,86,80,87,75,76,97,101,121,111,105,84,65,54,61,71,125,60,50,112,102,110,77,45,93,62,120,115,70,113,117,85,46,123,89,95,116,55,110,92,109,100,72,88,105,119,94,45,67,58,60,45,107,73,100,79,47,99,51,53,68,125,90,48,82,85,65,52,70,59,125,49,118,103,91,124,78,81,63,63),
x1 = c(52,52,93,78,117,64,103,81,75,82,70,71,92,96,116,106,100,79,60,49,56,66,120,55,45,107,97,105,72,40,88,57,115,110,65,108,112,80,41,118,84,90,111,50,105,87,104,95,67,83,100,114,89,40,62,53,55,40,102,68,95,74,42,94,46,48,63,120,85,43,77,80,60,47,65,54,120,44,113,98,86,119,73,76,58,58),
x2 = c(0.0001,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0001,0)
)

summary(dat)

scope <- list(
x1=c('1','x1','ns(x1, df=2)'),
x2=c('1','x2')
)

gam.object <- gam(y~1, data=dat)
step.object <- step.gam(gam.object, scope=scope, trace=2)
is.null(step.object)

# Run this if you want to evaluate one lm.wfit example call made from within gam functions
x <- cbind(rep(1,nrow(dat)), dat$x1)
names(x) <- c('Intercept', 'x1')
z <- dat$y
w <- rep(1,nrow(dat))
str(eval(expression(lm.wfit(x, z, w, method = "qr", singular.ok = TRUE))))





######################################################################
#1 R 3.1.2 (x86_64-redhat-linux-gnu (64-bit)) / gam 1.09.1
######################################################################
> step.object <- step.gam(gam.object, scope=scope, trace=2)
Start: y ~ 1; AIC= 797.0034
Trial: y ~ x1 + 1 ; AIC= -5121.796
Trial: y ~ 1 + x2 ; AIC= 796.915
Step:1 y ~ x1 ; AIC= -5121.796
Trial: y ~ ns(x1, df=2) + 1 ; AIC= -5123.816
Trial: y ~ x1 + x2 ; AIC= -5174.408
Step:2 y ~ x1 + x2 ; AIC= -5174.408
Trial: y ~ ns(x1, df=2) + x2 ; AIC= -5164.829
> is.null(step.object)
[1] FALSE

######################################################################
#2: R 3.3.0 (x86_64-redhat-linux-gnu (64-bit)) / gam 1.12
######################################################################
> step.object <- step.gam(gam.object, scope=scope, trace=2)
Start: y ~ 1; AIC= 797.0034
Trial: y ~ x1 + 1 ; AIC= -5121.796
Trial: y ~ 1 + x2 ; AIC= 796.915
Step:1 y ~ x1 ; AIC= -5121.796
Trial: y ~ ns(x1, df=2) + 1 ; AIC= -5123.816
Trial: y ~ x1 + x2 ; AIC= -5174.408
Step:2 y ~ x1 + x2 ; AIC= -5174.408
Trial: y ~ ns(x1, df=2) + x2 ; AIC= -5164.829
> is.null(step.object)
[1] FALSE

######################################################################
#3: R 3.3.0 (x86_64-redhat-linux-gnu (64-bit)) / gam 1.14-4
######################################################################
> step.object <- step.gam(gam.object, scope=scope, trace=2)
Start: y ~ 1; AIC= 797.0034
Trial: y ~ x1 + 1 ; AIC= -5121.796
Trial: y ~ 1 + x2 ; AIC= 796.915
Step:1 y ~ x1 ; AIC= -5121.796
Trial: y ~ ns(x1, df=2) + 1 ; AIC= -5123.816
Trial: y ~ x1 + x2 ; AIC= -5174.408
Step:2 y ~ x1 + x2 ; AIC= -5174.408
Trial: y ~ ns(x1, df=2) + x2 ; AIC= -5164.829
> is.null(step.object)
[1] FALSE

######################################################################
#4: R 3.3.2 (x86_64-redhat-linux-gnu (64-bit)) / gam 1.14-4
######################################################################
> step.object <- step.gam(gam.object, scope=scope, trace=2)
Start: y ~ 1; AIC= 797.0034
Trial: y ~ x1 + 1 ; AIC= -5122.886
Trial: y ~ 1 + x2 ; AIC= 796.915
Step:1 y ~ x1 ; AIC= -5122.886
Trial: y ~ ns(x1, df=2) + 1 ; AIC= -5177.531
Trial: y ~ x1 + x2 ; AIC= -5177.531
Step:2 y ~ ns(x1, df = 2) ; AIC= -5177.531
Trial: y ~ ns(x1, df=2) + x2 ; AIC= -5222.703
Step:3 y ~ ns(x1, df = 2) + x2 ; AIC= -5222.703
> is.null(step.object)
[1] TRUE

######################################################################
#5: R 3.3.2 (x86_64-redhat-linux-gnu (64-bit)) / gam 1.09.1
######################################################################
> step.object <- step.gam(gam.object, scope=scope, trace=2)
Start: y ~ 1; AIC= 797.0034
Trial: y ~ x1 + 1 ; AIC= -5122.886
Trial: y ~ 1 + x2 ; AIC= 796.915
Step:1 y ~ x1 ; AIC= -5122.886
Trial: y ~ ns(x1, df=2) + 1 ; AIC= -5177.531
Trial: y ~ x1 + x2 ; AIC= -5177.531
Step:2 y ~ ns(x1, df = 2) ; AIC= -5177.531
Trial: y ~ ns(x1, df=2) + x2 ; AIC= -5222.703
Step:3 y ~ ns(x1, df = 2) + x2 ; AIC= -5222.703
> is.null(step.object)
[1] TRUE

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.