predict.lm () mit einer unbekannten Faktorstufe in Testdaten

Ich passe ein Modell an, um Daten zu faktorisieren und zu prognostizieren. Wenn die newdata in ” predict.lm() eine einzige Faktorenebene enthalten, die dem Modell unbekannt ist, schlägt die function ” predict.lm() fehl und gibt einen Fehler zurück.

Gibt es eine gute Möglichkeit, predict.lm() eine Vorhersage für jene Faktorstufen zu erhalten, die das Modell kennt, und NA für unbekannte Faktorstufen, anstatt nur einen Fehler?

Beispielcode:

 foo <- data.frame(response=rnorm(3),predictor=as.factor(c("A","B","C"))) model <- lm(response~predictor,foo) foo.new <- data.frame(predictor=as.factor(c("A","B","C","D"))) predict(model,newdata=foo.new) 

Ich möchte, dass der allerletzte Befehl drei “echte” Vorhersagen zurückgibt, die den Faktorstufen “A”, “B” und “C” entsprechen, und eine NA , die dem unbekannten Pegel “D” entspricht.

Die function wurde von MorgenBall aufgeräumt und erweitert. Es ist jetzt auch in sperrorest implementiert.

Zusatzfunktionen

  • löscht nicht verwendete Faktorstufen, anstatt nur die fehlenden Werte auf NA .
  • gibt eine Nachricht an den Benutzer aus, dass Faktorstufen gelöscht wurden
  • prüft auf das Vorhandensein von Faktorvariablen in test_data und gibt originale Daten zurück.frame falls nicht vorhanden
  • funktioniert nicht nur für lm , glm und auch für glmmPQL

Hinweis: Die hier gezeigte function kann sich im Laufe der Zeit ändern (verbessern).

 #' @title remove_missing_levels #' @description Accounts for missing factor levels present only in test data #' but not in train data by setting values to NA #' #' @import magrittr #' @importFrom gdata unmatrix #' @importFrom stringr str_split #' #' @param fit fitted model on training data #' #' @param test_data data to make predictions for #' #' @return data.frame with matching factor levels to fitted model #' #' @keywords internal #' #' @export remove_missing_levels < - function(fit, test_data) { # https://stackoverflow.com/a/39495480/4185785 # drop empty factor levels in test data test_data %>% droplevels() %>% as.data.frame() -> test_data # 'fit' object structure of 'lm' and 'glmmPQL' is different so we need to # account for it if (any(class(fit) == "glmmPQL")) { # Obtain factor predictors in the model and their levels factors < - (gsub("[-^0-9]|as.factor|\\(|\\)", "", names(unlist(fit$contrasts)))) # do nothing if no factors are present if (length(factors) == 0) { return(test_data) } map(fit$contrasts, function(x) names(unmatrix(x))) %>% unlist() -> factor_levels factor_levels %>% str_split(":", simplify = TRUE) %>% extract(, 1) -> factor_levels model_factors < - as.data.frame(cbind(factors, factor_levels)) } else { # Obtain factor predictors in the model and their levels factors <- (gsub("[-^0-9]|as.factor|\\(|\\)", "", names(unlist(fit$xlevels)))) # do nothing if no factors are present if (length(factors) == 0) { return(test_data) } factor_levels <- unname(unlist(fit$xlevels)) model_factors <- as.data.frame(cbind(factors, factor_levels)) } # Select column names in test data that are factor predictors in # trained model predictors <- names(test_data[names(test_data) %in% factors]) # For each factor predictor in your data, if the level is not in the model, # set the value to NA for (i in 1:length(predictors)) { found <- test_data[, predictors[i]] %in% model_factors[ model_factors$factors == predictors[i], ]$factor_levels if (any(!found)) { # track which variable var <- predictors[i] # set to NA test_data[!found, predictors[i]] <- NA # drop empty factor levels in test data test_data %>% droplevels() -> test_data # issue warning to console message(sprintf(paste0("Setting missing levels in '%s', only present", " in test data but missing in train data,", " to 'NA'."), var)) } } return(test_data) } 

Wir können diese function auf das Beispiel in der Frage wie folgt anwenden:

 predict(model,newdata=remove_missing_levels (fit=model, test_data=foo.new)) 

Während ich versuchte, diese function zu verbessern, stieß ich auf die Tatsache, dass SL Lernmethoden wie lm , glm usw. die gleichen Ebenen in Train & Test svm , während ML Lernmethoden ( svm , randomForest ) fehlschlagen, wenn die Level entfernt werden. Diese Methoden benötigen alle Ebenen in Training und Test.

Eine allgemeine Lösung ist ziemlich schwer zu erreichen, da jedes angepasste Modell eine andere Möglichkeit hat, ihre Faktor-Level-Komponente zu speichern ( fit$xlevels für lm und fit$contrasts für glmmPQL ). Zumindest scheint es bei lm verwandten Modellen konsistent zu sein.

Sie müssen die zusätzlichen Ebenen vor jeder Berechnung entfernen, wie:

 > id < - which(!(foo.new$predictor %in% levels(foo$predictor))) > foo.new$predictor[id] < - NA > predict(model,newdata=foo.new) 1 2 3 4 -0.1676941 -0.6454521 0.4524391 NA 

Dies ist eine allgemeinere Art, es zu tun, es wird alle Ebenen, die nicht in den ursprünglichen Daten auftreten, auf NA setzen. Wie Hadley in den Kommentaren erwähnt, hätten sie sich dafür entscheiden können, dies in die function predict() einzufügen, was sie jedoch nicht taten

Warum Sie das tun müssen, wird offensichtlich, wenn Sie sich die Berechnung selbst ansehen. Intern werden die Vorhersagen wie folgt berechnet:

 model.matrix(~predictor,data=foo) %*% coef(model) [,1] 1 -0.1676941 2 -0.6454521 3 0.4524391 

Am unteren Rand haben Sie beide Modellmatrizen. Sie sehen, dass die für foo.new eine extra Spalte hat, so dass Sie die Matrixberechnung nicht mehr verwenden können. Wenn Sie das neue Dataset zum Modellieren verwenden würden, würden Sie auch ein anderes Modell erhalten, das eine zusätzliche Dummy-Variable für das zusätzliche Level enthält.

 > model.matrix(~predictor,data=foo) (Intercept) predictorB predictorC 1 1 0 0 2 1 1 0 3 1 0 1 attr(,"assign") [1] 0 1 1 attr(,"contrasts") attr(,"contrasts")$predictor [1] "contr.treatment" > model.matrix(~predictor,data=foo.new) (Intercept) predictorB predictorC predictorD 1 1 0 0 0 2 1 1 0 0 3 1 0 1 0 4 1 0 0 1 attr(,"assign") [1] 0 1 1 1 attr(,"contrasts") attr(,"contrasts")$predictor [1] "contr.treatment" 

Sie können nicht nur die letzte Spalte aus der Modellmatrix löschen, denn auch wenn Sie dies tun, werden beide anderen Ebenen noch beeinflusst. Der Code für Level A wäre (0,0). Für B dies (1,0), für C dies (0,1) … und für D es wieder (0,0)! Ihr Modell würde also annehmen, dass A und D das gleiche Niveau haben, wenn es die letzte Dummy-Variable naiv fallen lässt.

Zu einem eher theoretischen Teil: Es ist möglich, ein Modell ohne alle Ebenen zu bauen. Wie ich vorher erklären wollte, ist dieses Modell nur für die Ebenen gültig, die Sie beim Erstellen des Modells verwendet haben. Wenn Sie auf neue Ebenen stoßen, müssen Sie ein neues Modell erstellen, um die zusätzlichen Informationen aufzunehmen. Wenn Sie das nicht tun, können Sie nur die zusätzlichen Ebenen aus dem Dataset löschen. Aber dann verlierst du im Grunde genommen alle Informationen, die darin enthalten sind, also wird es im Allgemeinen nicht als gute Übung betrachtet.

Wenn Sie nach dem Erstellen Ihres LM – Modells mit den fehlenden Ebenen in Ihren Daten fertig werden wollen, aber bevor Sie Predicted aufrufen (da wir nicht genau wissen, welche Level im Voraus fehlen könnten), habe ich hier eine function erstellt, die alle Ebenen nicht in der Modell zu NA – die Vorhersage liefert dann auch NA und Sie können dann eine alternative Methode zur Vorhersage dieser Werte verwenden.

Objekt wird Ihre lm Ausgabe von lm sein (…, data = trainData)

Daten sind der Datenrahmen, für den Sie Vorhersagen erstellen möchten

 missingLevelsToNA< -function(object,data){ #Obtain factor predictors in the model and their levels ------------------ factors<-(gsub("[-^0-9]|as.factor|\\(|\\)", "",names(unlist(object$xlevels)))) factorLevels<-unname(unlist(object$xlevels)) modelFactors<-as.data.frame(cbind(factors,factorLevels)) #Select column names in your data that are factor predictors in your model ----- predictors<-names(data[names(data) %in% factors]) #For each factor predictor in your data if the level is not in the model set the value to NA -------------- for (i in 1:length(predictors)){ found<-data[,predictors[i]] %in% modelFactors[modelFactors$factors==predictors[i],]$factorLevels if (any(!found)) data[!found,predictors[i]]<-NA } data } 

Klingt nach zufälligen Effekten. Schauen Sie in etwas wie glmer (lme4-Paket). Mit einem Bayes-Modell erhalten Sie Effekte, die sich 0 annähern, wenn nur wenige Informationen zur Schätzung verwendet werden. Warnung, dass Sie selbst eine Vorhersage machen müssen, anstatt vorherzusagen ().

Alternativ können Sie einfach Dummy-Variablen für die Ebenen erstellen, die Sie in das Modell einbeziehen möchten, z. B. eine Variable 0/1 für Montag, eine für Dienstag, eine für Mittwoch usw. Sonntag wird automatisch aus dem Modell entfernt, wenn es alle enthält 0’s. Aber eine 1 in der Sonntagsspalte in den anderen Daten zu haben, wird den Vorhersageschritt nicht durchfallen lassen. Es wird nur angenommen, dass der Sonntag einen Effekt hat, der an den anderen Tagen durchschnittlich ist (was wahr sein kann oder auch nicht).

Eine der Annahmen von linearen / logistischen Regressionen ist zu wenig oder keine Multi-Kollinearität; Wenn also die Prädiktorvariablen idealerweise unabhängig voneinander sind, muss das Modell nicht alle möglichen Faktorstufen sehen. Ein neuer Faktor-Level (D) ist ein neuer Prädiktor und kann auf NA gesetzt werden, ohne die Vorhersagefähigkeit der verbleibenden Faktoren A, B, C zu beeinflussen. Deshalb sollte das Modell immer noch in der Lage sein, Vorhersagen zu treffen. Aber das Hinzufügen der neuen Ebene D verwirft das erwartete Schema. Das ist das ganze Problem. Einstellung NA behebt das.

Das lme4 Paket wird neue Ebenen behandeln, wenn Sie das Flag allow.new.levels=TRUE beim Aufruf von predict .

Beispiel: Wenn Ihr Wochentag-Faktor in einem Variablen- dow und einem kategorischen Ergebnis b_fail , könnten Sie laufen

M0 < - lmer(b_fail ~ x + (1 | dow), data=df.your.data, family=binomial(link='logit')) M0.preds <- predict(M0, df.new.data, allow.new.levels=TRUE)

Dies ist ein Beispiel mit einer zufälligen logistischen Regressionsanalyse. Natürlich können Sie eine reguläre Regression durchführen ... oder die meisten GLM-Modelle. Wenn Sie den Bayes-Pfad weiter hinunter gehen wollen, schauen Sie sich Gelman & Hill's exzellentes Buch und die Stan- Infrastruktur an.

Eine schnelle Lösung für Split-Tests besteht darin, seltene Werte als “andere” zu rekodieren. Hier ist eine Implementierung:

 rare_to_other < - function(x, fault_factor = 1e6) { # dirty dealing with rare levels: # recode small cells as "other" before splitting to train/test, # assuring that lopsided split occurs with prob < 1/fault_factor # (Nb not fully kosher, but useful for quick and dirty exploratory). if (is.factor(x) | is.character(x)) { min.cell.size = log(fault_factor, 2) + 1 xfreq <- sort(table(x), dec = T) rare_levels <- names(which(xfreq < min.cell.size)) if (length(rare_levels) == length(unique(x))) { warning("all levels are rare and recorded as other. make sure this is desirable") } if (length(rare_levels) > 0) { message("recoding rare levels") if (is.factor(x)) { altx < - as.character(x) altx[altx %in% rare_levels] <- "other" x <- as.factor(altx) return(x) } else { # is.character(x) x[x %in% rare_levels] <- "other" return(x) } } else { message("no rare levels encountered") return(x) } } else { message("x is neither a factor nor a character, doing nothing") return(x) } } 

Bei data.table würde der Aufruf beispielsweise lauten:

 dt[, (xcols) := mclapply(.SD, rare_to_other), .SDcol = xcols] # recode rare levels as other 

wo xcols ist eine beliebige Teilmenge von colnames(dt) .