Манипулирование моделью логистической регрессии

Я хотел бы понять, что делает следующий код. Человек, который написал код, больше не работает здесь, и он почти полностью без документов. Меня попросили исследовать это кто-то, кто думает, что это "модель байесовской логистической регрессии "

bglm <- function(Y,X) {
    # Y is a vector of binary responses
    # X is a design matrix

    fit <- glm.fit(X,Y, family = binomial(link = logit))
    beta <- coef(fit)
    fs <- summary.glm(fit)
    M <- t(chol(fs$cov.unscaled))
    betastar <- beta + M %*% rnorm(ncol(M))
    p <- 1/(1 + exp(-(X %*% betastar)))
    return(runif(length(p)) <= p)
}

Я могу видеть, что она соответствует логистической модели, принимает транспонирование факторизации Чолсеки оцененной ковариационной матрицы, умножает это на вектор ничьих из и затем добавляет к оценкам модели. Затем это предварительно умножается на матрицу проекта, берется обратный логит, сравниваемый с вектором отрисовок из и возвращаемый результирующий двоичный вектор. Но что все это значит статистически?$N(0,1)$$U(0,1)$

Что делает функция:
По сути, функция генерирует новые данные псевдослучайного ответа (т. ) из модели ваших данных. Используемая модель является стандартной моделью для частых исследований. Как обычно, предполагается, что ваши данные * являются известными константами - они никак не выбираются. Что я считаю важной особенностью этой функции, так это то, что она включает неопределенность в отношении оценочных параметров. $Y$$X$

* Обратите внимание, что вы должны вручную добавить вектор из в качестве крайнего левого столбца вашей матрицы перед вводом его в функцию, если вы не хотите подавить перехват (что, как правило, не очень хорошая идея).$1$$X$

Какой смысл этой функции:
я не знаю честно. Это могло бы быть частью процедуры байесовского MCMC, но это был бы только один кусок - вам понадобилось бы больше кода в другом месте, чтобы фактически выполнить байесовский анализ. Я не чувствую себя достаточно опытным в методах Байеса, чтобы комментировать это окончательно, но функция не «чувствует» меня как то, что обычно используется.

Он также мог бы использоваться при анализе мощности на основе моделирования. (См. Мой ответ здесь: моделирование анализа мощности логистической регрессии - эксперименты , спроектированные для получения информации об этом типе вещей.) Стоит отметить, что анализ мощности, основанный на предыдущих данных, которые не учитывают неопределенность оценок параметров, часто оптимистичный. (Я обсуждаю этот момент здесь: желаемый размер эффекта против ожидаемого размера эффекта .)

$Y$

simulationParameters <- function(Y,X) {
                        # Y is a vector of binary responses
                        # X is a design matrix, you don't have to add a vector of 1's 
                        #   for the intercept

                        X    <- cbind(1, X)  # this adds the intercept for you
                        fit  <- glm.fit(X,Y, family = binomial(link = logit))
                        beta <- coef(fit)
                        fs   <- summary.glm(fit)
                        M    <- t(chol(fs$cov.unscaled))

                        return(list(betas=beta, uncertainties=M))
}

simulateY <- function(X, betas, uncertainties, ncolM, N){

             # X      <- cbind(1, X)  # it will be slightly faster if you input w/ 1's
             # ncolM  <- ncol(uncertainties) # faster if you input this
             betastar <- betas + uncertainties %*% rnorm(ncolM)
             p        <- 1/(1 + exp(-(X %*% betastar)))

             return(rbinom(N, size=1, prob=p))
}