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

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')

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^{+} = (X' X)^{-1} X'$$X$$\hat{\beta} = X^{+} Y$$H = X X^{+}$$\hat{Y} = H Y$$X$$Y$, 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^{+}$$H$$QR$$Y$lm()

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.

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)))