Logistic regression model manipulation

Question

I would like to understand what the following code is doing. The person who wrote the code no longer works here and it is almost completely undocumented. I was asked to investigate it by someone who thinks "it's a bayesian logistic regression model"

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

I can see that it fits a logistic model, takes the transpose of the Cholseky factorisation of the estimated covariance matrix, post-multiplies this by a vector of draws from $N(0,1)$ and is then added to the model estimates. This is then premultiplied by the design matrix, the inverse logit of this is taken, compared with a vector of draws from $U(0,1)$ and the resulting binary vector returned. But what does all this mean statistically ?

Answer 1

What the function does:
In essence, the function generates new pseudorandom response (i.e., $Y$) data from a model of your data. The model being used is a standard frequentist model. As is customary, it is assuming that your $X$* data are known constants--they are not sampled in any way. What I see as the important feature of this function is that it is incorporating uncertainty about the estimated parameters.

* Note that you have to manually add a vector of $1$'s as the leftmost column of your $X$ matrix before inputting it to the function, unless you want to suppress the intercept (which is generally not a good idea).

What was the point of this function:
I don't honestly know. It could have been part of a Bayesian MCMC routine, but it would only have been one piece--you would need more code elsewhere to actually run a Bayesian analysis. I don't feel sufficiently expert on Bayesian methods to comment definitively on this, but the function doesn't 'feel' to me like what would typically be used.

It could also have been used in simulation-based power analyses. (See my answer here: Simulation of logistic regression power analysis - designed experiments, for information on this type of thing.) It is worth noting that power analyses based on prior data that do not take the uncertainty of the parameter estimates into account are often optimistic. (I discuss that point here: Desired effect size vs. expected effect size.)

If you want to use this function:
As @whuber notes in the comments, this function will be inefficient. If you want to use this to do (for example) power analyses, I would split the function into two new functions. The first would read in your data and output the parameters and the uncertainties. The second new function would generate the new pseudorandom $Y$ data. The following is an example (although it may be possible to improve it further):

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