Manipulacja modelem regresji logistycznej

Chciałbym zrozumieć, co robi następujący kod. Osoba, która napisała kod, już tu nie pracuje i jest prawie całkowicie nieudokumentowana. Zostałem poproszony o zbadanie go przez kogoś, kto myśli „ to bayesowski model regresji logistycznej ”

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

Widzę, że pasuje do modelu logistycznego, przyjmuje transpozycję faktoryzacji Cholseky'ego szacowanej macierzy kowariancji, pomnoża ją przez wektor losowań z a następnie dodaje do oszacowań modelu. Jest to następnie mnożone przez macierz projektową, pobierana jest odwrotna logit tego, w porównaniu z wektorem losowań z i zwracany wynikowy wektor binarny. Ale co to wszystko oznacza statystycznie?$N(0,1)$$U(0,1)$

Działanie funkcji:
Zasadniczo funkcja generuje nowe dane odpowiedzi pseudolosowej (tj. ) z modelu danych. Używany model jest standardowym modelem dla osób często podróżujących. Jak zwykle, zakłada się, że twoje dane X * są znanymi stałymi - nie są w żaden sposób próbkowane. Ważną cechą tej funkcji jest to, że obejmuje ona niepewność co do szacowanych parametrów. $Y$$X$

* Pamiętaj, że musisz ręcznie dodać wektor jako lewą kolumnę macierzy X przed wprowadzeniem go do funkcji, chyba że chcesz stłumić przecięcie (co na ogół nie jest dobrym pomysłem).$1$$X$

Jaki był sens tej funkcji:
nie wiem szczerze. Mogłoby to być częścią rutyny MCMC Bayesa, ale byłby to tylko jeden kawałek - potrzebowałbyś więcej kodu gdzie indziej, aby przeprowadzić analizę bayesowską. Nie czuję się wystarczająco ekspertem w zakresie metod bayesowskich, aby definitywnie wypowiedzieć się na ten temat, ale funkcja nie „wydaje mi się”, jak by to było zwykle.

Mógłby być również wykorzystany w symulacyjnych analizach mocy. (Zobacz moją odpowiedź tutaj: Symulacja analizy mocy regresji logistycznej - zaprojektowane eksperymenty , aby uzyskać informacje na temat tego rodzaju rzeczy.) Warto zauważyć, że analizy mocy oparte na wcześniejszych danych, które nie uwzględniają niepewności oszacowań parametrów, są często optymistyczny. (Omawiam ten punkt tutaj: pożądany rozmiar efektu vs. oczekiwany rozmiar efektu ).

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