¿Cómo maneja R los valores perdidos en lm?

Question

Me gustaría hacer una regresión de un vector B contra cada una de las columnas de una matriz A. Esto es trivial si no hay datos perdidos, pero si la matriz A contiene valores perdidos, entonces mi regresión contra A está restringida a incluir sólo las filas donde todos los valores están presentes (el valor por defecto na.omitir comportamiento). Esto produce resultados incorrectos para las columnas sin datos perdidos. Puedo hacer una regresión de la matriz de columnas B contra columnas individuales de la matriz A, pero tengo que hacer miles de regresiones, y esto es prohibitivamente lento y poco elegante. El na.excluir parece estar diseñada para este caso, pero no consigo que funcione. ¿Qué estoy haciendo mal? Usando R 2.13 en OSX, si importa.

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

Answer 1

Edición: He entendido mal su pregunta. Hay dos aspectos:

a) na.omit y na.exclude ambos hacen la eliminación de casos con respecto a los predictores y a los criterios. Sólo se diferencian en que las funciones del extractor como residuals() o fitted() rellenarán su salida con NA s para los casos omitidos con na.exclude , teniendo así una salida de la misma longitud que las variables de entrada.

> 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) El verdadero problema no es esta diferencia entre na.omit y na.exclude Parece que no quieres una eliminación por casos que tenga en cuenta las variables de criterio, que es lo que hacen ambos.

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

Los resultados de la regresión dependen de las matrices $X^{+} = (X' X)^{-1} X'$ (pseudoinverso de la matriz de diseño $X$ , coeficientes $\hat{\beta} = X^{+} Y$ ) y la matriz del sombrero $H = X X^{+}$ , valores ajustados $\hat{Y} = H Y$ ). Si no desea la eliminación de casos, necesita una matriz de diseño diferente $X$ para cada columna de $Y$ por lo que no hay forma de evitar ajustar regresiones separadas para cada criterio. Se puede intentar evitar la sobrecarga de lm() haciendo algo parecido a lo siguiente:

> 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

Tenga en cuenta que puede haber formas numéricamente mejores de calcular $X^{+}$ y $H$ , se podría comprobar un $QR$ -en lugar de la descomposición. El enfoque SVD se explica aquí en SE . No he cronometrado el enfoque anterior con matrices grandes $Y$ contra el uso real de lm() .

Answer 2

Se me ocurren dos formas. Una es combinar los datos utilizando el na.exclude y luego volver a separar los datos:

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)

Otra forma es utilizar el data y crear una fórmula.

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)

Si está haciendo mucha regresión, la primera forma debería ser más rápida, ya que se realiza menos magia de fondo. Aunque si sólo necesitas los coeficientes y los residuos te sugiero que utilices lsfit que es mucho más rápido que lm . La segunda forma es un poco más bonita, pero en mi portátil al intentar hacer un resumen de la regresión resultante da un error. Intentaré ver si se trata de un error.

Answer 3

El siguiente ejemplo muestra cómo hacer predicciones y residuos que se ajustan al marco de datos original (utilizando la opción "na.action=na.exclude" en lm() para especificar que los NA deben colocarse en los vectores de residuos y predicciones donde el marco de datos original tenía valores perdidos. También muestra cómo especificar si las predicciones deben incluir sólo las observaciones en las que tanto las variables explicativas como las dependientes estaban completas (es decir, predicciones estrictamente dentro de la muestra) o las observaciones en las que las variables explicativas estaban completas y, por lo tanto, la predicción Xb es posible, (es decir, incluyendo la predicción fuera de la muestra para las observaciones que tenían variables explicativas completas pero faltaba la variable dependiente).

Utilizo cbind para añadir las variables predichas y residuales al conjunto de datos original.

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