如何使用参数自举计算协方差矩阵？

问题描述

使用以下数据集“数据”：load(url("https://www.math.ntnu.no/emner/TMA4315/2020h/hoge-veluwe.Rdata")) 我已经安装了Poisson GLM model = glm(y~t + I(t^2),family = poisson,data)。我现在想通过使用带有1000个模拟的参数自举估计从GLM回归获得的β系数的协方差矩阵。到目前为止，我的代码是：


ysim = simulate(mod_quad,1)
betahat = matrix(0,nrow = 1000,ncol = 3)
for (i in 1:1000){  
  sim_data = cbind(ysim,data$t)
  betahat[i,]= glm(ysim ~ data$t + I(data$t^2),data = sim_data )$coefficients
  ysim = simulate(glm(ysim ~ data$t + I(data$t^2)family = poisson,data = sim_data ),1)
}

var(betahat[1000,])

每次都将betahat矩阵变成零矩阵，所以我不确定我的方法中缺少什么？

解决方法

这实际上是编码而不是统计问题。精简，我会这样做：

## set vals to NA (not 0) to make detection of problems easier
betahat <- matrix(NA,nrow = 1000,ncol = 3)
for (i in 1:1000) {
    ## replace response with parametric simulation
    sim_data <- transform(data,y=simulate(mod_quad,1)[[1]])
    ## refit model with new data
    newfit <- update(mod_quad,data=sim_data)
    ## store new coefficients
    betahat[i,] <- coef(newfit) 
}
## compute variance
var(betahat)

r r variance