Jak R radzi sobie z brakującymi wartościami w lm?

32

Chciałbym regresować wektor B względem każdej kolumny w macierzy A. Jest to trywialne, jeśli nie ma brakujących danych, ale jeśli macierz A zawiera brakujące wartości, to moja regresja w stosunku do A jest ograniczona i obejmuje tylko wiersze, w których wszystkie wartości są obecne (domyślne zachowanie na.omit ). To powoduje nieprawidłowe wyniki dla kolumn bez brakujących danych. Mogę regresować macierz kolumn B względem pojedynczych kolumn macierzy A, ale mam tysiące regresji do wykonania, a to jest zbyt wolne i nieeleganckie. Na.exclude funkcja wydaje się być zaprojektowany dla tej sprawy, ale nie mogę tego dokonać. Co robię tutaj źle? Używanie R 2.13 na OSX, jeśli ma to znaczenie.

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
źródło
1
Co rozumiesz przez „Mogę obliczyć każdy wiersz osobno”?
chl
Przepraszam, miałem na myśli: „Mogę regresować macierz kolumny B indywidualnie względem kolumn w kolumnie A”, co oznacza pojedyncze wywołania do lm. Edytowane, aby to odzwierciedlić.
David Quigley,
1
Jednorazowe połączenia z lm / regresją nie są świetnym sposobem na zrobienie regresji (przechodząc przez definicję regresji, która polega na znalezieniu częściowego wpływu każdego predyktora na odpowiedź / wynik, biorąc pod uwagę stan innego zmienne)
KarthikS

Odpowiedzi:

23

Edycja: źle zrozumiałem twoje pytanie. Istnieją dwa aspekty:

a) na.omiti na.excludeoba dokonują przypadkowego usunięcia w odniesieniu zarówno do predyktorów, jak i kryteriów. Różnią się tylko tym, że funkcje ekstraktora, takie jak residuals()lub fitted()wypełniają swoje wyjście za pomocą NAs dla pominiętych przypadków na.exclude, dzięki czemu mają wyjście o tej samej długości co zmienne wejściowe.

> 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) Prawdziwy problem nie tkwi w tej różnicy między, na.omiti na.excludewydaje się, że nie chcesz usuwania przypadków z uwzględnieniem zmiennych kryteriów, co robią oba.

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

X+=(XX)1XXβ^=X+YH=XX+Y^=HYXY, więc nie ma mowy o dopasowaniu osobnych regresji dla każdego kryterium. Możesz spróbować uniknąć narzutów lm(), wykonując czynności w następujący sposób:

> 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+HQRYlm()

karakal
źródło
Ma to sens, biorąc pod uwagę moje zrozumienie, jak powinien działać. Jeśli jednak wywołasz> X.both = cbind (X1, X2), a następnie> dim (lm (X.both ~ Y, na.action = na.wyklucz) $ resztki) nadal otrzymujesz 94 reszt, zamiast 97 i 97.
David Quigley,
Jest to poprawa, ale jeśli spojrzysz na wartości resztkowe (lm (X.both ~ Y, na.action = na.exclude)), zobaczysz, że w każdej kolumnie brakuje sześciu brakujących wartości, nawet jeśli brakuje wartości w kolumnie 1 X. oba pochodzą z innych próbek niż te w kolumnie 2. Wyklucza więc zachowanie kształtu matrycy reszt, ale pod maską R najwyraźniej cofa się tylko z wartościami obecnymi we wszystkich rzędach X. obu. Może to mieć dobry powód statystyczny, ale dla mojej aplikacji jest to problem.
David Quigley
@ David Nie zrozumiałem twojego pytania. Myślę, że teraz rozumiem twój punkt widzenia i zredagowałem moją odpowiedź, aby rozwiązać ten problem.
caracal
5

Mogę wymyślić dwa sposoby. Jednym z nich jest połączenie danych, na.excludea następnie ponowne rozdzielenie danych:

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)

Innym sposobem jest użycie dataargumentu i utworzenie formuły.

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)

Jeśli wykonujesz dużo regresji, pierwszy sposób powinien być szybszy, ponieważ wykonuje się mniej magii w tle. Chociaż jeśli potrzebujesz tylko współczynników i reszt, sugeruję użycie lsfit, co jest znacznie szybsze niż lm. Drugi sposób jest nieco ładniejszy, ale na moim laptopie próba podsumowania wynikowej regresji powoduje błąd. Spróbuję sprawdzić, czy to błąd.

mpiktas
źródło
Dzięki, ale lm (A.ex ~ B.ex) w twoim kodzie pasuje 9 punktów w stosunku do A1 (poprawnie) i 9 punktów w stosunku do A2 (niepożądany). Istnieje 10 punktów pomiarowych zarówno dla B1, jak i A2; Wyrzucam jeden punkt w regresji B1 przeciwko A2, ponieważ w A1 brakuje odpowiedniego punktu. Jeśli to po prostu sposób, w jaki to działa, mogę to zaakceptować, ale nie o to staram się zmusić R.
David Quigley,
@ David, och, wygląda na to, że źle zrozumiałem twój problem. Opublikuję poprawkę później.
mpiktas
1

Poniższy przykład pokazuje, jak tworzyć prognozy i reszty zgodne z oryginalną ramką danych (przy użyciu opcji „na.action = na.exclude” w lm () w celu określenia, że ​​NA należy umieścić w wektorach reszt i predykcji, w których oryginalna ramka danych brakowało wartości. Pokazuje także, jak określić, czy przewidywania powinny obejmować tylko obserwacje, w których zarówno zmienne objaśniające, jak i zależne były kompletne (tj. przewidywania ściśle w próbie) lub obserwacje, w których zmienne objaśniające były kompletne, a zatem możliwe jest przewidywanie Xb ( tj. łącznie z prognozowaniem poza próbą dla obserwacji, które miały pełne zmienne objaśniające, ale brakowało zmiennej zależnej).

Korzystam z cbind, aby dodać przewidywane i resztkowe zmienne do oryginalnego zestawu danych.

## 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
źródło