Comment R gère-t-il les valeurs manquantes dans lm?

32

Je voudrais régresser un vecteur B par rapport à chacune des colonnes d'une matrice A. C'est trivial s'il n'y a pas de données manquantes, mais si la matrice A contient des valeurs manquantes, ma régression par rapport à A est contrainte d'inclure uniquement les lignes où tout des valeurs sont présentes (le comportement na.omit par défaut ). Cela produit des résultats incorrects pour les colonnes sans données manquantes. Je peux régresser la matrice de colonnes B contre des colonnes individuelles de la matrice A, mais j'ai des milliers de régressions à faire, ce qui est excessivement lent et inélégant. La fonction na.exclude semble être conçue pour ce cas, mais je ne peux pas la faire fonctionner. Qu'est-ce que je fais mal ici? Utiliser R 2.13 sur OSX, si cela est important.

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
la source
1
Qu'entendez-vous par «je peux calculer chaque ligne individuellement»?
chl
Désolé, je voulais dire "Je peux régresser la matrice de colonnes B par rapport aux colonnes de A individuellement", ce qui signifie des appels uniques à lm. Modifié pour refléter cela.
David Quigley
1
Les appels ponctuels à lm / régression ne sont pas un excellent moyen de procéder à une régression (en suivant la définition de la régression, qui consiste à trouver l'effet partiel de chaque prédicteur sur une réponse / un résultat compte tenu de l'état des autres variables)
KarthikS

Réponses:

23

Edit: j'ai mal compris votre question. Il y a deux aspects:

a) na.omitet les na.excludedeux suppriment également les prédicteurs et les critères. Ils diffèrent seulement par le fait que l'extracteur fonctionne comme residuals()ou fitted()remplira leur sortie avec NAs pour les cas omis avec na.exclude, ayant ainsi une sortie de la même longueur que les variables d'entrée.

> 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) Le vrai problème n'est pas avec cette différence entre na.omitet na.exclude, vous ne semblez pas vouloir une suppression également qui prend en compte les variables de critère, ce que les deux font.

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

Les résultats de la régression dépendent des matrices (pseudo-inverse de la matrice de conception X , coefficientsX+=(XX)-1XX) et la matrice de chapeauH=XX+, valeurs ajustées Y =HY). Si vous ne voulez pas la suppression de casas, vous avez besoin d'une matrice de conception différenteXpour chaque colonne deYβ^=X+YH=XX+Y^=HYXY, il n'y a donc aucun moyen d'ajuster des régressions distinctes pour chaque critère. Vous pouvez essayer d'éviter les frais généraux de lm()en faisant quelque chose comme suit:

> 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

Notez qu'il pourrait y avoir de meilleures façons numériquement de calculer et H , vous pouvez plutôt vérifier une décomposition Q R. L'approche SVD est expliquée ici sur SE . Je n'ai pas chronométré l'approche ci-dessus avec de grandes matrices Y par rapport à l'utilisation réelle .X+HQRYlm()

caracal
la source
Cela a du sens étant donné ma compréhension du fonctionnement de na.exclude. Cependant, si vous appelez> X.both = cbind (X1, X2) puis> dim (lm (X.both ~ Y, na.action = na.exclude) $ residuels), vous obtenez toujours 94 résidus, au lieu de 97 et 97.
David Quigley
C'est une amélioration, mais si vous regardez les résidus (lm (X.both ~ Y, na.action = na.exclude)), vous voyez que chaque colonne a six valeurs manquantes, même si les valeurs manquantes dans la colonne 1 de X. les deux proviennent d'échantillons différents de ceux de la colonne 2. Ainsi, na.exclude préserve la forme de la matrice des résidus, mais sous le capot, R ne régresse apparemment qu'avec les valeurs présentes dans toutes les lignes de X.both. Il peut y avoir une bonne raison statistique à cela, mais pour mon application, c'est un problème.
David Quigley
@ David J'avais mal compris votre question. Je pense que je vois maintenant votre point, et j'ai édité ma réponse pour y répondre.
caracal
5

Je peux penser à deux façons. La première consiste à combiner les données à l'aide de la na.exclude, puis à nouveau séparer les données:

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)

Une autre façon consiste à utiliser l' dataargument et à créer une formule.

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 vous effectuez beaucoup de régression, la première méthode devrait être plus rapide, car moins de magie d'arrière-plan est effectuée. Bien que si vous n'avez besoin que de coefficients et de résidus, je suggère d'utiliser lsfit, ce qui est beaucoup plus rapide que lm. La deuxième façon est un peu plus agréable, mais sur mon ordinateur portable, essayer de faire un résumé sur la régression résultante génère une erreur. Je vais essayer de voir s'il s'agit d'un bug.

mpiktas
la source
Merci, mais lm (A.ex ~ B.ex) dans votre code correspond à 9 points contre A1 (correct) et 9 points contre A2 (indésirable). Il y a 10 points mesurés pour B1 et A2; Je rejette un point dans la régression de B1 contre A2 car le point correspondant manque dans A1. Si c'est juste la façon dont cela fonctionne, je peux l'accepter, mais ce n'est pas ce que j'essaie de faire faire à R.
David Quigley
@David, oh, on dirait que j'ai mal compris votre problème. Je publierai le correctif plus tard.
mpiktas
1

L'exemple suivant montre comment faire des prédictions et des résidus qui sont conformes à la trame de données d'origine (en utilisant l'option "na.action = na.exclude" dans lm () pour spécifier que les NA doivent être placés dans les vecteurs résiduels et de prédiction où la trame de données d'origine Il indique également comment spécifier si les prévisions ne doivent inclure que des observations où les variables explicatives et dépendantes sont complètes (c'est-à-dire des prédictions strictement dans l'échantillon) ou des observations où les variables explicatives sont complètes et, par conséquent, la prédiction Xb est possible ( c'est-à-dire, y compris la prédiction hors échantillon pour les observations qui avaient des variables explicatives complètes mais manquaient la variable dépendante).

J'utilise cbind pour ajouter les variables prédites et résiduelles à l'ensemble de données d'origine.

## 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
la source