LogLik von Hand aus einer logistischen Regression berechnen

Question

Ich lief eine gemischte Modell logistische Regression Anpassung meines Modells mit genetischen Beziehung Matrix mit einem R-Paket als GMMAT (Funktion: glmmkin()).

aus dem Modell My Ausgabe enthält (aus der Bedienungsanleitung entnommen):

• theta: die Dispersion Parameterschätzung [1] und der Varianz Komponentenparameterschätzwert [2]
• coefficients: fester Effektparameter Schätzungen (einschließlich des Abschnitts).
• linear.predictors: die linearen Prädiktoren.
• fitted.values: angepasste Mittelwerte auf der ursprünglichen Skala.
• Y: ein Vektor der Länge entspricht der Stichprobengröße für den endgültigen Arbeitsvektor.
• P: die Projektionsmatrix mit Abmessungen gleich der Stichprobengröße.
• residuals: Reste im Originalmaßstab. NICHT durch den Dispersionsparameter neu skaliert.
• cov: Kovarianzmatrix für die festen Effekte (einschließlich der Achsenabschnitt).
• converged: ein logischer Indikator für die Konvergenz.

Ich versuche, die Log-Likelihood zu erhalten, um die erklärte Varianz zu berechnen. Mein erster Instinkt war, die logLik.glm-Funktion auseinander zu ziehen, um das "von Hand" zu berechnen, und ich blieb stecken, als ich versuchte, AIC zu berechnen. Ich habe die Antwort von here verwendet.

Ich habe eine Plausibilitätsprüfung mit einer logistischen Regression durchgeführt mit stats::glm(), wo die model1$aic ist 4013.232, aber mit dem Stack Overflow Antwort, die ich gefunden habe, ich 30613.03 erhalten.

Meine Frage ist - weiß jemand, wie Log-Likelihood von einer logistischen Regression von Hand mit der Ausgabe, die ich oben in R aufgeführt habe, zu berechnen?

Answer 1

Keine statistischen Erkenntnisse hier, nur die Lösung sehe ich von glm.fit. Dies funktioniert nur, wenn Sie nicht Gewichte angegeben haben, während die Modelle passend

get_logLik <- function(s_model, family = binomial(logit)) { 
    n <- length(s_model$y) 
    wt <- rep(1, n) # or s_model$prior_weights if field exists 
    deviance <- sum(family$dev.resids(s_model$y, s_model$fitted.values, wt)) 
    mod_rank <- sum(!is.na(s_model$coefficients)) # or s_model$rank if field exists 

    aic <- family$aic(s_model$y, rep(1, n), s_model$fitted.values, wt, deviance) + 2 * mod_rank 
    log_lik <- mod_rank - aic/2 
    return(log_lik) 
} 

Zum Beispiel ...

model <- glm(vs ~ mpg, mtcars, family = binomial(logit)) 
logLik(model) 
# 'log Lik.' -12.76667 (df=2) 

sparse_model <- model[c("theta", "coefficients", "linear.predictors", "fitted.values", "y", "P", "residuals", "cov", "converged")] 
get_logLik(sparse_model) 
#[1] -12.76667