[Rd] potential bug in simulate.lm when gaussian(link != "identity")

[Alexandre Courtiol]

Dear all,

I think that there is a bug in the function simulate.lm() when called upon
a glm fitted with gaussian family with a link other than identity. The
variance of the simulated response is clearly off and an inspection to the
code of simulate.lm reveals that it is because the variance is divided by
the model weights (precisely, not the prior ones), which is not documented.

Can somebody file this bug for me (if you agree that this is a bug)?
Many thanks.
Alex

Simple demonstration:

set.seed(1L)
y <- 10 + rnorm(n = 100)
mean(y) ##  10.10889
var(y)  ##   0.8067621

mod_glm1  <- glm(y ~ 1, family = gaussian())
new.y1 <- simulate(mod_glm1)[, 1]
mean(new.y1) ## 10.07493
var(new.y1)  ##  0.7402303

mod_glm2  <- glm(y ~ 1, family = gaussian(link = "log"))
new.y2 <- simulate(mod_glm2)[, 1]
mean(new.y2) ## 10.11152
var(new.y2)  ##  0.008445062  ##### WRONG #####

mod_glm3 <- mod_glm2
mod_glm3$weights <- NULL
new.y3 <- simulate(mod_glm3)[, 1]
mean(new.y3) ## 10.15524
var(new.y3)  ##  0.7933856  ##### OK(?) #####


### my session ###

> sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-apple-darwin16.5.0 (64-bit)
Running under: macOS Sierra 10.12.4

Matrix products: default
BLAS:
/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK:
/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libLAPACK.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

loaded via a namespace (and not attached):
[1] compiler_3.4.0 tools_3.4.0

-- 
Alexandre Courtiol

http://sites.google.com/site/alexandrecourtiol/home

*"Science is the belief in the ignorance of experts"*, R. Feynman

	[[alternative HTML version deleted]]