Regressions-Subset-Algorithmus in R

Question

Ich möchte einen Regressions-Subset-Algorithmus in R für ein "Beta-Regressions" -Modell erstellen. Es gibt ein Paket betareg in R, das für Beta-Regressionen passt, und was mich interessiert, ist das Modell, das die "Log-Wahrscheinlichkeit" maximiert.

Grundsätzlich funktioniert dies, indem Sie das beste k-Faktor-Regressionsmodell für k = 1,2, ..., p auswählen, wobei p die Anzahl der Variablen ist, die Sie haben.

Zum Beispiel, wenn ich x_1, x_2, x_3 als meine Variablen und y als meine Antwort habe. Ich möchte etwas haben, das ist:

Schritt 1: Finden Sie die besten 1-Faktor-Modell

mod1 <- betareg(y~x_1, data = test) 
mod1.sum <- summary(mod1) 

mod2 <- betareg(y~x_2, data = test) 
mod2.sum <- summary(mod2) 

mod3 <- betareg(y~x_3, data = test) 
mod3.sum <- summary(mod3) 

jetzt, dass ich alle Modelle passen müssen, möchte ich die Log-Likelihood von jedem vergleichen:

likelihoods <- c(mod1.sum$loglik, mod2.sum$loglik, mod3.sum$loglik) 
which.max(likelihoods) 

Schritt 2: Finden Sie den besten Faktor zum Hinzufügen zum besten 1-Faktor-Modell. Angenommen, x_1 war im vorherigen Schritt der beste. Dann vergleichen wir in diesem Schritt das Modell mit x_1 und x_2 mit dem Modell mit x_1 und x_3 mit demjenigen mit der größten Log-Wahrscheinlichkeit.

Schritt 3: Nehmen Sie die besten beiden Variablen als gegeben, finden Sie die dritte Variable, die den größten Anstieg zur Protokollwahrscheinlichkeit beiträgt.

Schritt 4: Geben Sie das beste 1-Faktor-Modell, das beste 2-Faktor-Modell, ..., das beste p-Faktor-Modell, die eingeschlossenen Faktoren und ihre entsprechenden Log-Wahrscheinlichkeiten zurück.

Ich kämpfe dies effizient zu tun, wenn p groß ist, sagen rund 40

Answer 1

Neben meiner anderen Antwort, die zeigt, wie mit kofnGA für Beta-Regression Best-Subset Auswahl zu tun, ich bin auch ein Beispiel dafür, wie vorwärts Auswahl von Hand zu tun.

Wir wieder durch das Laden des Pakets und Daten beginnen.

library("betareg") 
data("FoodExpenditure", package = "betareg") 

Ich bin auch Listen mit der Antwort sowie alle Regressoren Einrichtung (einschließlich 40 zufällig diejenigen (Beachten Sie, dass im Gegensatz zu den anderen bin ich Halten Sie den Abschnitt in x, die hier bequemer ist.

)

fe_data <- list(
    y = with(FoodExpenditure, food/income), 
    x = model.matrix(~ income + persons, data = FoodExpenditure) 
) 
set.seed(123) 
fe_data$x <- cbind(fe_data$x, 
    matrix(rnorm(40 * nrow(fe_data$x)), ncol = 40)) 
colnames(fe_data$x)[4:43] <- paste0("x", 1:40) 

Dann setzen wir wieder eine Funktion für die negative Log-Likelihood-up (aber nicht das Intercept manuell schließen müssen, weil sie noch in x) ist.

nll <- function(v, data) -betareg.fit(x = data$x[, v, drop = FALSE], 
    y = data$y)$loglik 

Dann speichern wir den Index aller möglichen Regressoren vall und unsere Suche mit dem Intercept (v <- 1) und die entsprechenden negativen Log-Likelihood (n) initialisieren.

vall <- 1:ncol(fe_data$x) 
v <- 1 
n <- nll(v, data = fe_data) 

Und dann iterieren wir unsere Vorwärtssuche für 15 Iterationen (numerische Instabilitäten auf diesem kleinen Datensatz für eine höhere Anzahl von Variablen zu vermeiden). Wir wählen immer diese zusätzliche Variable, die den negativen log-likelihood meisten abnimmt:

for(i in 1:15) { 
    vi <- vall[-v] 
    ni <- sapply(vi, function(vii) nll(v = c(v, vii), data = fe_data)) 
    v <- c(v, vi[which.min(ni)]) 
    n <- c(n, ni[which.min(ni)]) 
} 

Die Reihenfolge, in der die Variablen ausgewählt werden, ist die folgende. Beachten Sie, dass die echten Regressoren zuerst ausgewählt werden, gefolgt von den zufälligen Rauschregressoren. (Aber versuchen zu set.seed(1), über dem Zufall Regressoren vor den realen umfassen wird ...)

colnames(fe_data$x)[v] 
## [1] "(Intercept)" "income"  "persons"  "x28"   "x18"   
## [6] "x29"   "x22"   "x11"   "x5"   "x8"   
## [11] "x38"   "x24"   "x13"   "x23"   "x36"   
## [16] "x16"   

Die entsprechende Abnahme des negativen Log-Likelihood und zugehörige BIC kann sichtbar gemacht werden:

m <- seq_along(v) 
plot(m, n, type = "b", 
    xlab = "Number of regressors", ylab = "Log-likelihood") 
plot(m, n + log(nrow(fe_data$x)) * (m + 1), type = "b", 
    xlab = "Number of regressors", ylab = "BIC") 

Das würde also das Modell mit den drei echten Regressoren als das beste BIC-Modell wählen.

Answer 2

Hier ist eine alternative Lösung ohne betareg zu verwenden. Die Ausgabe ist ähnlich und für die Berücksichtigung Ihrer Probleme.

Hier ist der Datensatz ich benutzt hatte:

set.seed(12345) 
dat <- data.frame(y=runif(50), x_1=runif(50), x_2=runif(50), x_3=runif(50)) 

Sprünge Bibliothek Mit einer Liste aller möglichen Kombinationen zu erstellen:

library(leaps) 
subs<-regsubsets(y~., data=dat, nbest=10, nvmax=100, really.big=T) 
subs<-summary(subs)$which[,-1] 
all.mods<-lapply(1:nrow(subs), function(x)paste("y", paste(names(which(subs[x,])), collapse="+"), sep="~")) 

all.mods 

[[1]] 
[1] "y~x_2" 

[[2]] 
[1] "y~x_1" 

[[3]] 
[1] "y~x_3" 

[[4]] 
[1] "y~x_2+x_3" 

[[5]] 
[1] "y~x_1+x_2" 

[[6]] 
[1] "y~x_1+x_3" 

[[7]] 
[1] "y~x_1+x_2+x_3" 

lineare Regression für alle Modelle Laufen:

all.lm<-lapply(all.mods, function(x)lm(as.formula(x), data=dat)) 

Überprüfen Sie logLikihood für jedes Modell:

lapply(all.lm, logLik) 

[[1]] 
'log Lik.' -7.051835 (df=3) 

[[2]] 
'log Lik.' -9.288776 (df=3) 

[[3]] 
'log Lik.' -9.334048 (df=3) 

[[4]] 
'log Lik.' -6.904604 (df=4) 

[[5]] 
'log Lik.' -7.051584 (df=4) 

[[6]] 
'log Lik.' -9.215915 (df=4) 

[[7]] 
'log Lik.' -6.888849 (df=5) 

Answer 3

Nach meinem besten Wissen gibt es keine dedizierte effiziente Implementierung der besten Teilmenge Auswahl für Beta-Regression (in R oder anders). Es gibt jedoch einige generische Implementierungen, die annähernde Lösungen dafür bereitstellen, z. B. basierend auf genetischen Algorithmen wie dem kofnGA Paket (auf CRAN und veröffentlicht in JSS). Siehe das Beispiel unten. (Um vorwärts Suche anstelle von Best-Subset Auswahl zu benutzen, siehe meine andere Antwort.)

Alternativ können Sie ein (verallgemeinerte) lineares Modell verwenden, das in etwa, was betareg tut und verwenden Sie die Subgruppenauswahl für die . Zum Beispiel könnten Sie die Antwort logit-transformieren (d. H., qlogis(y)) und dann die Auswahl der besten Untergruppe unter Verwendung einer linearen Regression über leaps (CRAN) oder lmSubsets (R-Forge). Oder Sie könnten einen GLM mit family = quasibinomial verwenden und glmulti (CRAN, JSS) verwenden. Und dann könnten Sie das Ergebnis der besten Teilmenge aus diesem ungefähren Modell verwenden und es in einer Beta-Regression verwenden. Natürlich gibt dies nicht das beste Beta-Regressionsergebnis, aber es könnte ein nützlicher Ausgangspunkt für weitere Analysen sein.

Daher zurück zum direkten genetischen Algorithmus für die Beta-Regression. Um darzustellen, wie dies mit kofnGA wir erste Lastpakete und Beispieldaten erreicht werden:

library("betareg") 
library("kofnGA") 
data("FoodExpenditure", package = "betareg") 

Dann bauen wir eine Liste mit den Antwortvariablen y und einer Regressormatrix x. Beachten Sie, dass wir den Abschnitt hier weglassen, um ihn später in das Modell zu zwingen (d. H. Der Abschnitt sollte nicht der Auswahl unterliegen).

fe_data <- list(
    y = with(FoodExpenditure, food/income), 
    x = model.matrix(~ income + persons, data = FoodExpenditure)[, -1] 
) 

Zusätzlich zu den beiden Regressoren oben Wir richten nun 40 Zufallsrauschen Variablen auf den Regressormatrix

fe_data$x <- cbind(fe_data$x, 
    matrix(rnorm(40 * nrow(fe_data$x)), ncol = 40)) 
colnames(fe_data$x)[3:42] <- paste0("x", 1:40) 

Nun fügen wir kofnGA verwenden, um das beste Modell mit 2 Regressoren aus der Auswahl mögliche 42 Regressoren (plus der immer eingeschlossene Abschnitt). Als kofnGA minimiert eine Obejctive wir verwenden die negative Log-Likelihood von betareg zur Verfügung gestellt. Die Arbeitspferd Funktion betareg.fit statt betareg verwendet wird unnötige Formel zu vermeiden Parsen usw.

nll <- function(v, data) -betareg.fit(x = cbind(1, data$x[, v]), 
    y = data$y)$loglik 

Schließlich führen wir den genetischen Algorithmus nur für 100 Generationen eine gewisse Rechenzeit in diesem kurzen Beispiel zu speichern:

set.seed(1) 
ga <- kofnGA(n = 42, k = 2, OF = nll, data = fe_data, ngen = 100) 

die resultierende Ausgabe ist

summary(ga) 
## Genetic algorithm search, 100 generations 
## Number of unique solutions in the final population: 1 
## 
## Objective function values: 
##      average minimum 
## Initial population -36.56597 -41.74583 
## Final population -45.33351 -45.33351 
## 
## Best solution (found at generation 1): 
## 1 2 

So wird in dieser sehr einfachen künstlichen Einrichtung in der Tat den genetischen Algorithmus pi cks die ersten 2 Regressoren (von den echten Daten) und nicht irgendwelche der irrelevanten zufälligen 40 Regressoren, die wir hinzufügten. So konnten wir gehen jetzt voran und setzen ein richtiges Beta Regressionsmodell mit den Regressoren

colnames(fe_data$x)[ga$bestsol] 
## [1] "income" "persons" 

usw. Beachten Sie, dass die Beta-Regression hier einfach verwendet wird, verwendet einen festen Präzisionsparameter (mit Log-Link). Wenn Sie eine variable Dispersion wünschen, müssen Sie nll entsprechend ändern.