Я генерирую данные, используя двухкомпонентную модель препятствий через распределение Пуассона. Часть этой модели генерирует данные подсчета только из 0 и 1. Я использую метод начальной загрузки, чтобы сгенерировать достаточное количество наблюдений, чтобы найти ошибку и мощность в моей модели, чтобы увидеть, подходит ли эта модель
Мой вопрос просто: «Как мне посчитать количество 0 в каждой итерации?». Иногда R выдает ошибку на одной из итераций, которая имеет вид Error in solve.default(as.matrix(fit_zero$hessian)) : Lapack routine dgesv: system is exactly singular: U[1,1] = 0
. Это приводит меня к выводу, что матрица имеет либо все 0, либо все 1. Следующим моим шагом было попытаться посчитать количество нулей в каждой итерации. Однако я не уверен, как это сделать. Две проблемы, с которыми я столкнулся, заключались в том, как получить каждую итерацию, если я делаю 1000 итераций с 1000 бутстрапами, и как посчитать (и вывести) количество нулей в циклах. Моя попытка ниже:
errors <- 0
no.zeros <- 0
results <- foreach(iiii=icount(iterations), .combine = rbind) %do%{
message("Parameter set: ",parameters.index, "/", nrow(kParameters$matrix), '\t', "Iteration: ", iiii, "/", iterations)
data = data.frame(gen.hurdle(n, a, b, b, c, c, i, i, i))
data0 = data.frame(gen.hurdle(n, a, 0, 0, c, c, i, i, i))
p_0 =1-mean(data$z_p)
p_0_hat =1-mean(data$tstar)
p_0_b0 =1-mean(data0$z_p)
p_0_hat_b0 =1-mean(data0$tstar)
#power
fit1 = lm(m ~ x, data=data)
fit2 = hurdle(formula = y ~ m + x, data=data, dist = "poisson", zero.dist = "binomial")
a_hat = summary(fit1)$coef[2,1]
b1_hat = summary(fit2)[[1]]$count[2,1]
b2_hat = summary(fit2)[[1]]$zero[2,1]
ab1_hat = prod(a_hat,b1_hat)
ab2_hat = prod(a_hat,b2_hat)
#type I error
fit3 = lm(m ~ x, data=data0)
fit4 = hurdle(formula = y ~ m + x, data=data0, dist = "poisson", zero.dist = "binomial")
a_hat_b0 = summary(fit3)$coef[2,1]
b1_hat_b0 = summary(fit4)[[1]]$count[2,1]
b2_hat_b0 = summary(fit4)[[1]]$zero[2,1]
ab1_hat_b0 = prod(a_hat_b0,b1_hat_b0)
ab2_hat_b0 = prod(a_hat_b0,b2_hat_b0)
message(paste0("Bootstrapping..."))
# bootstrap
boot <- foreach(jjjj = icount(r), .combine = rbind, .errorhandling = "remove", .packages = c("pscl")) %dopar%{
#power
boot.data = data[sample(nrow(data), replace = TRUE), ]
has.zero <- prod(boot.data$y) > 0
if(!has.zero) {
no.zeros <- no.zeros + 1
boot.data$y[1] = 0
warning(paste0("Iteration #",iiii, " Bootstrap #",jjjj, " had no zeros!"), immediate. = TRUE, call. = FALSE)
}