Manipulacja modelem regresji logistycznej

12

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)

P Sellaz
źródło
Prawdopodobnie bardzo pomogłoby wiedzieć, w jakim polu jest to używane ...
naught101
2
Zasadniczo funkcja generuje dane z (częstego) modelu twoich danych, uwzględniając niepewność co do rzeczywistych parametrów. Może to być część rutyny MCMC Bayesa, ale może być również wykorzystane w opartej na symulacji analizie mocy (nb, analizy mocy oparte na wcześniejszych danych, które nie uwzględniają niepewności, są często optymistyczne ).
gung - Przywróć Monikę
Nie ma za co, @PSellaz. Ponieważ nikt inny nie odpowiedział, zamienię to w „oficjalną” odpowiedź.
gung - Przywróć Monikę

Odpowiedzi:

7

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. YX

* 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).1X

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))
}
gung - Przywróć Monikę
źródło
4
+1. Dla mnie dziwne jest to, że dopasowanie i symulowane przewidywania są wykonywane w obrębie pojedynczej funkcji. Zwykle takie operacje można wykonać najpierw obliczając dopasowanie (powrót betai M), a następnie tworząc liczne symulacje iid oparte na tym dopasowaniu. (Umieszczenie ich w tej samej funkcji niepotrzebnie spowodowałoby powtarzanie dopasowania za każdym razem, znacznie spowalniając obliczenia.) Z tych symulacji można było odzyskać ( między innymi ) przedziały prognozowania dla nieliniowych lub bardzo skomplikowanych kombinacji odpowiedzi.
whuber
@ whuber, zgadzam się. Właściwie zastanawiałem się nad zasugerowaniem, aby funkcję podzielić na 2 różne funkcje przed użyciem jej do symulacji, ale nie wyglądało na to, żeby była to część pytania.
gung - Przywróć Monikę