Manipulation des logistischen Regressionsmodells

Ich möchte verstehen, was der folgende Code tut. Die Person, die den Code geschrieben hat, arbeitet hier nicht mehr und ist fast vollständig undokumentiert. Ich wurde gebeten, es von jemandem zu untersuchen, der denkt " es ist ein bayesianisches logistisches Regressionsmodell ".

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

Ich kann sehen, dass es zu einem logistischen Modell passt, die Transponierung der Cholseky-Faktorisierung der geschätzten Kovarianzmatrix übernimmt, diese mit einem Vektor von Ziehungen aus nachmultipliziert und dann zu den Modellschätzungen hinzugefügt wird. Dies wird dann durch die Entwurfsmatrix vormultipliziert, die inverse Protokollierung davon wird genommen, verglichen mit einem Vektor von Zügen von U ( 0 , 1 ) und dem resultierenden binären Vektor, der zurückgegeben wird. Aber was bedeutet das alles Mittelwert statistisch?$N(0,1)$$U(0,1)$

$Y$$X$

* Beachten Sie, dass Sie manuell einen Vektor von hinzufügen müssen$1$ als die Spalte ganz links von Ihnen$X$ Matrix vor der Eingabe in die Funktion, es sei denn, Sie möchten den Achsenabschnitt unterdrücken (was im Allgemeinen keine gute Idee ist).

Was war der Sinn dieser Funktion:
Ich weiß es nicht ehrlich. Es hätte Teil einer Bayes'schen MCMC-Routine sein können, aber es wäre nur ein Teil gewesen - Sie würden an anderer Stelle mehr Code benötigen, um tatsächlich eine Bayes'sche Analyse durchzuführen. Ich fühle mich nicht ausreichend mit Bayes'schen Methoden vertraut, um dies definitiv zu kommentieren, aber die Funktion fühlt sich für mich nicht so an, wie sie normalerweise verwendet wird.

Es könnte auch in simulationsbasierten Leistungsanalysen verwendet worden sein. (Siehe meine Antwort hier: Simulation der logistischen Regressionsleistungsanalyse - entworfene Experimente , um Informationen zu dieser Art von Dingen zu erhalten.) Es ist erwähnenswert, dass Leistungsanalysen, die auf früheren Daten basieren und die Unsicherheit der Parameterschätzungen nicht berücksichtigen, häufig sind optimistisch. (Ich diskutiere diesen Punkt hier: Gewünschte Effektgröße vs. erwartete Effektgröße .)

Wenn Sie diese Funktion verwenden möchten:
Wie @whuber in den Kommentaren feststellt, ist diese Funktion ineffizient. Wenn Sie dies (zum Beispiel) für Leistungsanalysen verwenden möchten, würde ich die Funktion in zwei neue Funktionen aufteilen. Der erste würde Ihre Daten einlesen und die Parameter und die Unsicherheiten ausgeben. Die zweite neue Funktion würde das neue Pseudozufall erzeugen$Y$Daten. Das Folgende ist ein Beispiel (obwohl es möglicherweise möglich ist, es weiter zu verbessern):

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