Manipulação do modelo de regressão logística

Gostaria de entender o que o código a seguir está fazendo. A pessoa que escreveu o código não trabalha mais aqui e é quase completamente indocumentada. Me pediram para investigá-lo por alguém que pensasse " é um modelo de regressão logística bayesiana "

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)
}

Eu posso ver que ele se encaixa em um modelo logístico, aceita a transposição da fatoração de Cholseky da matriz de covariância estimada, multiplica isso por um vetor de empates de e é então adicionado às estimativas do modelo. Isso é pré-multiplicado pela matriz de design, o logit inverso disso é obtido, comparado com um vetor de draws de e o vetor binário resultante retornado. Mas o que tudo isso significa estatisticamente?$N(0,1)$$U(0,1)$

O que a função faz:
Em essência, a função gera novos dados de resposta pseudo-aleatória (ou seja, ) a partir de um modelo dos seus dados. O modelo usado é um modelo freqüentador padrão. Como é habitual, supõe-se que seus dados * sejam constantes conhecidas - eles não são amostrados de forma alguma. O que vejo como característica importante dessa função é que ela está incorporando incerteza sobre os parâmetros estimados. $Y$$X$

* Observe que você deve adicionar manualmente um vetor como a coluna mais à esquerda da sua matriz antes de inseri-lo na função, a menos que queira suprimir a interceptação (o que geralmente não é uma boa ideia).$1$$X$

Qual era o objetivo dessa função:
não sei honestamente. Poderia ter sido parte de uma rotina bayesiana do MCMC, mas seria apenas uma peça - você precisaria de mais código em outro local para realmente executar uma análise bayesiana. Não me sinto suficientemente especialista em métodos bayesianos para comentar definitivamente sobre isso, mas a função não me parece o que normalmente seria usado.

Também poderia ter sido usado em análises de potência baseadas em simulação. (Veja minha resposta aqui: Simulação de análise de poder de regressão logística - experimentos projetados , para obter informações sobre esse tipo de coisa.) Vale a pena notar que as análises de poder com base em dados anteriores que não levam em consideração a incerteza das estimativas de parâmetros são frequentemente otimista. (Discuto esse ponto aqui: tamanho do efeito desejado vs. tamanho do efeito esperado .)

Se você quiser usar esta função:
Como o @whuber observa nos comentários, essa função será ineficiente. Se você quiser usar isso para fazer (por exemplo) análises de energia, eu dividiria a função em duas novas funções. O primeiro leria seus dados e produziria os parâmetros e as incertezas. A segunda nova função geraria os novos dados pseudo-aleatórios . A seguir, é apresentado um exemplo (embora seja possível melhorá-lo ainda mais): $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))
}