Hello, All:
The default "simulate" method for lm and glm seems to ignore the
sampling variance of the parameter estimates; see the trivial lm and
glm examples below. Both these examples estimate a mean with formula =
x~1. In both cases, the variance of the estimated mean is 1.
* In the lm example with x0 = c(-1, 1), var(x0) = 2, and
var(unlist(simulate(lm(x0~1), 10000, 1))) is 2.0064. Shouldn't it be 3
= var(mean(x0)) + var(x0) = (2/2) + 2?
* In the glm example with x1=1,
var(unlist(simulate(glm(x1~1, poisson), 10000, 1))) = 1.006. Shouldn't
it be 2 = var(glm estimate of the mean) + var(simulated Poisson
distribution) = 1 + 1?
I'm asking, because I've recently written "simulate" methods for
objects of class stats::glm and BMA::bic.glm, where my primary interest
was simulating the predicted mean with "newdata". I'm doing this, so I
can get Monte Carlo prediction intervals. My current code for
"simulate.glm" and "simulate.bic.glm" are available in the development
version of the "Ecfun" package on GitHub:
https://github.com/sbgraves237/Ecfun
Comparing my new code with "stats:::simulate.lm" raises the
following questions in my mind regarding "simulate" of a fit object:
1. Shouldn't "simulate" start by simulating the random
variability in the estimated parameters? I need that for my current
application. If a generic "simulate" function should NOT include this,
what should we call something that does include this? And how does the
current stats:::simulate.lm behavior fit with this?
2. Shouldn't "simulate" of a model fit include an option
for "newdata"? I need that for my application.
3. By comparing with "predict.glm", I felt I needed an
argument 'type = c("link", "response")'. "predict.glm" has an argument
'type = c("link", "response", "terms")'. I didn't need "terms", so I
didn't take the time to code that. However, a general "simulate"
function should probably include that.
4. My application involves assumed Poisson counts. I need
to simulate those as well. If I combined those with "simulate.glm",
what would I call them? I can't use the word "response", because that's
already used with a different meaning. Might "observations" be the
appropriate term?
What do you think?
Thanks,
Spencer Graves
> x0 <- c(-1, 1)
> var(x0)
[1] 2
> fit0 <- lm(x0~1)
> vcov(fit0)
(Intercept)
(Intercept) 1
> sim0 <- simulate(fit0, 10000, 1)
> var(unlist(sim0))
[1] 2.006408
> x1 <- 1
> fit1 <- glm(x1~1, poisson)
> coef(fit1)
(Intercept)
4.676016e-11
> exp(coef(fit1))
(Intercept)
1
> vcov(fit1)
(Intercept)
(Intercept) 0.9999903
> sim1 <- simulate(fit1, 10000, 1)
> var(unlist(sim1))
[1] 1.00617
> sessionInfo()
R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.2
Matrix products: default
BLAS:
/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK:
/Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.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.6.2 tools_3.6.2
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
