Manipulación del modelo de regresión logística.

Me gustaría entender qué está haciendo el siguiente código. La persona que escribió el código ya no funciona aquí y está casi completamente indocumentado. Me pidieron que lo investigara alguien que piensa " es un modelo de regresión logística bayesiano "

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

Puedo ver que se ajusta a un modelo logístico, toma la transposición de la factorización de Cholseky de la matriz de covarianza estimada, post-multiplica esto por un vector de sorteos de y luego se agrega a las estimaciones del modelo. Esto luego se multiplica por la matriz de diseño, se toma el logit inverso de esto, en comparación con un vector de sorteos de y se devuelve el vector binario resultante. Pero lo que hace todo esto media estadísticamente?$N(0,1)$$U(0,1)$

Qué hace la función:
en esencia, la función genera nuevos datos de respuesta pseudoaleatoria (es decir, ) a partir de un modelo de sus datos. El modelo que se utiliza es un modelo frecuentista estándar. Como es habitual, se supone que sus datos * son constantes conocidas, no se muestrean de ninguna manera. Lo que veo como la característica importante de esta función es que está incorporando incertidumbre sobre los parámetros estimados. $Y$$X$

* Tenga en cuenta que debe agregar manualmente un vector de como la columna más a la izquierda de su matriz antes de ingresarlo a la función, a menos que desee suprimir la intersección (que generalmente no es una buena idea).$1$$X$

¿Cuál era el punto de esta función
? Sinceramente, no lo sé. Podría haber sido parte de una rutina Bayesian MCMC, pero solo habría sido una pieza: necesitaría más código en otro lugar para ejecutar un análisis bayesiano. No me siento lo suficientemente experto en métodos bayesianos para comentar definitivamente sobre esto, pero la función no me "parece" como lo que normalmente se usaría.

También podría haberse utilizado en análisis de potencia basados ​​en simulación. (Vea mi respuesta aquí: Simulación de análisis de potencia de regresión logística: experimentos diseñados , para obtener información sobre este tipo de cosas). Vale la pena señalar que los análisis de potencia basados ​​en datos anteriores que no tienen en cuenta la incertidumbre de las estimaciones de los parámetros son a menudo optimista. (Discuto ese punto aquí: tamaño del efecto deseado versus tamaño del efecto esperado ).

Si desea utilizar esta función:
como señala @whuber en los comentarios, esta función será ineficiente. Si desea utilizar esto para hacer (por ejemplo) análisis de potencia, dividiría la función en dos nuevas funciones. El primero leería en sus datos y generaría los parámetros y las incertidumbres. La segunda nueva función generaría los nuevos datos pseudoaleatorios . El siguiente es un ejemplo (aunque puede ser posible mejorarlo aún más): $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))
}