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

Question

Me gustaría entender lo que el código siguiente está haciendo. La persona que escribió el código ya no funciona aquí, y es casi completamente indocumentados. Me pidió investigar por alguien que piensa que "es un modelo de regresión 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)
}

Puedo ver que se ajuste a un modelo logístico, se lleva a la transpuesta de la Cholseky factorización de la estimación de la matriz de covarianza, post-multiplica esta cifra por un vector de sorteos de $N(0,1)$ y se añade a continuación las estimaciones del modelo. Esto es luego multiplicado por el diseño de la matriz, a la inversa logit de esto es tomado, en comparación con un vector de sorteos de $U(0,1)$ y el binario resultante del vector devuelto. Pero ¿qué tiene todo esto significa que estadísticamente ?

Answer 1

Lo que hace la función:
En esencia, la función genera nuevas pseudoaleatoria de respuesta (es decir, $Y$) los datos de un modelo de los datos. El modelo utilizado es un estándar de modelo frecuentista. Como es habitual, se asume que su $X$* los datos son conocidas las constantes--no se muestra en ninguna manera. Lo que yo veo como la característica importante de esta función es que es la incorporación de la incertidumbre acerca de los parámetros estimados.

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

¿Cuál fue el punto de esta función:
Yo honestamente no sé. Podría haber sido parte de un Bayesiano MCMC rutina, pero sólo habría sido una pieza--necesita más código en otra parte, para ejecutar un análisis Bayesiano. No me siento lo suficientemente experto en métodos Bayesianos para comentar definitivamente en esto, pero la función no "sentir" a mí me gusta lo que normalmente se usaría.

También podría haber sido utilizado en la simulación de poder basadas en el análisis. (Véase mi respuesta aquí: Simulación de la regresión logística, el análisis de poder - los experimentos diseñados, para obtener información sobre este tipo de cosas.) Vale la pena señalar que el poder de análisis basados en datos anteriores, que no toman la incertidumbre de las estimaciones de los parámetros en cuenta a menudo son optimistas. (Me discutir ese punto aquí: efecto Deseado tamaño vs espera que el tamaño del efecto.)

Si desea utilizar esta función:
Como @whuber notas en los comentarios, esta función será ineficaz. Si desea utilizar este para hacer (por ejemplo) el poder de los análisis, me gustaría dividir la función en dos nuevas funciones. El primero sería la lectura de los datos y la salida de los parámetros y de las incertidumbres. La segunda nueva función de generar el nuevo pseudoaleatoria $Y$ de los datos. El siguiente es un ejemplo (aunque es posible mejorar aún más):

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