Wie geht R mit fehlenden Werten in lm um?

32

Ich möchte einen Vektor B für jede der Spalten in einer Matrix A regressieren. Dies ist trivial, wenn keine Daten fehlen. Wenn die Matrix A jedoch fehlende Werte enthält, darf meine Regression für A nur Zeilen enthalten, in denen alle enthalten sind Werte sind vorhanden (das Standardverhalten von na.omit ). Dies führt zu falschen Ergebnissen für Spalten ohne fehlende Daten. Ich kann die Spaltenmatrix B gegen einzelne Spalten der Matrix A regressieren, aber ich habe Tausende von Regressionen zu tun, und dies ist unerschwinglich langsam und unelegant. Die Funktion na.exclude scheint für diesen Fall entwickelt worden zu sein, aber ich kann sie nicht zum Laufen bringen. Was mache ich hier falsch? Verwenden Sie R 2.13 unter OSX, wenn es darauf ankommt.

A = matrix(1:20, nrow=10, ncol=2)
B = matrix(1:10, nrow=10, ncol=1)
dim(lm(A~B)$residuals)
# [1] 10 2 (the expected 10 residual values)

# Missing value in first column; now we have 9 residuals
A[1,1] = NA  
dim(lm(A~B)$residuals)
#[1]  9 2 (the expected 9 residuals, given na.omit() is the default)

# Call lm with na.exclude; still have 9 residuals
dim(lm(A~B, na.action=na.exclude)$residuals)
#[1]  9 2 (was hoping to get a 10x2 matrix with a missing value here)

A.ex = na.exclude(A)
dim(lm(A.ex~B)$residuals)
# Throws an error because dim(A.ex)==9,2
#Error in model.frame.default(formula = A.ex ~ B, drop.unused.levels = TRUE) : 
#  variable lengths differ (found for 'B')
David Quigley
quelle
1
Was meinen Sie mit "Ich kann jede Zeile einzeln berechnen"?
chl
Tut mir leid, ich wollte sagen "Ich kann die Spaltenmatrix B gegen die Spalten in A einzeln regressieren", was bedeutet, dass nacheinander Aufrufe an lm erfolgen. Bearbeitet, um dies widerzuspiegeln.
David Quigley
1
Einzelne Aufrufe von lm / regression sind kein guter Weg, um eine Regression durchzuführen (gemäß der Definition der Regression, bei der die partielle Auswirkung jedes Prädiktors auf eine Antwort / ein Ergebnis in Anbetracht des Zustands des anderen ermittelt wird) Variablen)
KarthikS

Antworten:

23

Edit: Ich habe deine Frage falsch verstanden. Es gibt zwei Aspekte:

a) na.omitund na.excludebeide löschen nacheinander in Bezug auf Prädiktoren und Kriterien. Sie unterscheiden sich nur darin, dass Extraktorfunktionen ihre Ausgabe für die ausgelassenen Fälle mit s auffüllen residuals()oder fitted()auffüllen und somit eine Ausgabe mit der gleichen Länge wie die Eingabevariablen haben.NAna.exclude

> N    <- 20                               # generate some data
> y1   <- rnorm(N, 175, 7)                 # criterion 1
> y2   <- rnorm(N,  30, 8)                 # criterion 2
> x    <- 0.5*y1 - 0.3*y2 + rnorm(N, 0, 3) # predictor
> y1[c(1, 3,  5)] <- NA                    # some NA values
> y2[c(7, 9, 11)] <- NA                    # some other NA values
> Y    <- cbind(y1, y2)                    # matrix for multivariate regression
> fitO <- lm(Y ~ x, na.action=na.omit)     # fit with na.omit
> dim(residuals(fitO))                     # use extractor function
[1] 14  2

> fitE <- lm(Y ~ x, na.action=na.exclude)  # fit with na.exclude
> dim(residuals(fitE))                     # use extractor function -> = N
[1] 20  2

> dim(fitE$residuals)                      # access residuals directly
[1] 14  2

b) Das eigentliche Problem liegt nicht in diesem Unterschied zwischen na.omitund na.exclude. Sie möchten anscheinend keine zufällige Löschung, bei der Kriterienvariablen berücksichtigt werden, was beide tun.

> X <- model.matrix(fitE)                  # design matrix
> dim(X)                                   # casewise deletion -> only 14 complete cases
[1] 14  2

X+=(XX)-1XXβ^=X+Y.H=XX+Y.^=HY.XY.Es führt also kein Weg daran vorbei, für jedes Kriterium separate Regressionen festzulegen. Sie können versuchen, den Overhead von zu vermeiden, lm()indem Sie folgendermaßen vorgehen:

> Xf <- model.matrix(~ x)                    # full design matrix (all cases)
# function: manually calculate coefficients and fitted values for single criterion y
> getFit <- function(y) {
+     idx   <- !is.na(y)                     # throw away NAs
+     Xsvd  <- svd(Xf[idx , ])               # SVD decomposition of X
+     # get X+ but note: there might be better ways
+     Xplus <- tcrossprod(Xsvd$v %*% diag(Xsvd$d^(-2)) %*% t(Xsvd$v), Xf[idx, ])
+     list(coefs=(Xplus %*% y[idx]), yhat=(Xf[idx, ] %*% Xplus %*% y[idx]))
+ }

> res <- apply(Y, 2, getFit)    # get fits for each column of Y
> res$y1$coefs
                   [,1]
(Intercept) 113.9398761
x             0.7601234

> res$y2$coefs
                 [,1]
(Intercept) 91.580505
x           -0.805897

> coefficients(lm(y1 ~ x))      # compare with separate results from lm()
(Intercept)           x 
113.9398761   0.7601234 

> coefficients(lm(y2 ~ x))
(Intercept)           x 
  91.580505   -0.805897

X+HQ.RY.lm()

caracal
quelle
Das macht Sinn, wenn ich verstehe, wie na.exclude funktionieren sollte. Wenn Sie jedoch> X.both = cbind (X1, X2) und dann> dim (lm (X.both ~ Y, na.action = na.exclude) $ residuals) aufrufen, erhalten Sie immer noch 94 Residuen anstelle von 97 und 97.
David Quigley
Das ist eine Verbesserung, aber wenn Sie sich die Residuen ansehen (lm (X.both ~ Y, na.action = na.exclude)), sehen Sie, dass jede Spalte sechs fehlende Werte hat, obwohl die fehlenden Werte in Spalte 1 von X sind. beide stammen aus anderen Proben als in Spalte 2. na.exclude behält also die Form der Residuenmatrix bei, aber unter der Haube bildet sich R anscheinend nur mit Werten zurück, die in allen Zeilen von X.both vorhanden sind. Es mag einen guten statistischen Grund dafür geben, aber für meine Anwendung ist es ein Problem.
David Quigley
@ David Ich hatte deine Frage falsch verstanden. Ich denke, ich verstehe jetzt Ihren Standpunkt und habe meine Antwort bearbeitet, um sie zu adressieren.
caracal
5

Ich kann mir zwei Möglichkeiten vorstellen. Man kombiniert die Daten mit na.excludeund trennt dann die Daten wieder:

A = matrix(1:20, nrow=10, ncol=2)
colnames(A) <- paste("A",1:ncol(A),sep="")

B = matrix(1:10, nrow=10, ncol=1)
colnames(B) <- paste("B",1:ncol(B),sep="")

C <- cbind(A,B)

C[1,1] <- NA
C.ex <- na.exclude(C)

A.ex <- C[,colnames(A)]
B.ex <- C[,colnames(B)]

lm(A.ex~B.ex)

Eine andere Möglichkeit besteht darin, das dataArgument zu verwenden und eine Formel zu erstellen.

Cd <- data.frame(C)
fr <- formula(paste("cbind(",paste(colnames(A),collapse=","),")~",paste(colnames(B),collapse="+"),sep=""))

lm(fr,data=Cd)

Cd[1,1] <-NA

lm(fr,data=Cd,na.action=na.exclude)

Wenn Sie viel Regression betreiben, sollte der erste Weg schneller sein, da weniger Hintergrundmagie ausgeführt wird. Wenn Sie jedoch nur Koeffizienten und Residuen benötigen, empfehle ich die Verwendung von lsfit, die viel schneller ist als lm. Der zweite Weg ist ein bisschen netter, aber der Versuch, auf meinem Laptop eine Zusammenfassung der resultierenden Regression zu erstellen, löst einen Fehler aus. Ich werde versuchen zu sehen, ob dies ein Fehler ist.

mpiktas
quelle
Danke, aber lm (A.ex ~ B.ex) in Ihrem Code passt 9 Punkte gegen A1 (richtig) und 9 Punkte gegen A2 (unerwünscht). Es gibt 10 gemessene Punkte für B1 und A2; Ich werfe einen Punkt in der Regression von B1 gegen A2 aus, weil der entsprechende Punkt in A1 fehlt. Wenn das nur so ist, kann ich das akzeptieren, aber ich versuche nicht, R dazu zu bringen, das zu tun.
David Quigley
@ David, oh, es sieht so aus, als hätte ich dein Problem falsch verstanden. Ich werde das Update später posten.
mpiktas
1

Das folgende Beispiel zeigt, wie Vorhersagen und Residuen erstellt werden, die dem ursprünglichen Datenrahmen entsprechen (unter Verwendung der Option "na.action = na.exclude" in lm (), um anzugeben, dass NAs in den Residuen- und Vorhersagevektoren platziert werden sollen, in denen sich der ursprüngliche Datenrahmen befindet Es wird auch gezeigt, wie angegeben werden kann, ob Vorhersagen nur Beobachtungen enthalten sollten, bei denen sowohl erklärende als auch abhängige Variablen vollständig waren (dh ausschließlich Vorhersagen in der Stichprobe), oder Beobachtungen, bei denen die erklärenden Variablen vollständig waren und daher eine Xb-Vorhersage möglich ist ( dh einschließlich einer Vorhersage außerhalb der Stichprobe für Beobachtungen mit vollständigen erklärenden Variablen, denen jedoch die abhängige Variable fehlte).

Ich benutze cbind, um die vorhergesagten und verbleibenden Variablen zum ursprünglichen Datensatz hinzuzufügen.

## Set up data with a linear model
N <- 10
NXmissing <- 2 
X <- runif(N, 0, 10)
Y <- 6 + 2*X + rnorm(N, 0, 1)
## Put in missing values (missing X, missing Y, missing both)
X[ sample(1:N , NXmissing) ] <- NA
Y[ sample(which(is.na(X)), 1)]  <- NA
Y[ sample(which(!is.na(X)), 1)]  <- NA
(my.df <- data.frame(X,Y))

## Run the regression with na.action specified to na.exclude
## This puts NA's in the residual and prediction vectors
my.lm  <- lm( Y ~ X, na.action=na.exclude, data=my.df)

## Predict outcome for observations with complete both explanatory and
## outcome variables, i.e. observations included in the regression
my.predict.insample  <- predict(my.lm)

## Predict outcome for observations with complete explanatory
## variables.  The newdata= option specifies the dataset on which
## to apply the coefficients
my.predict.inandout  <- predict(my.lm,newdata=my.df)

## Predict residuals 
my.residuals  <- residuals(my.lm)

## Make sure that it binds correctly
(my.new.df  <- cbind(my.df,my.predict.insample,my.predict.inandout,my.residuals))

## or in one fell swoop

(my.new.df  <- cbind(my.df,yhat=predict(my.lm),yhato=predict(my.lm,newdata=my.df),uhat=residuals(my.lm)))
Michael Ash
quelle