Bagaimana R menangani nilai yang hilang dalam lm?

32

Saya ingin mundur vektor B terhadap masing-masing kolom dalam matriks A. Ini sepele jika tidak ada data yang hilang, tetapi jika matriks A berisi nilai yang hilang, maka regresi saya terhadap A dibatasi untuk menyertakan hanya baris di mana semua nilai ada ( perilaku default na.omit ). Ini menghasilkan hasil yang salah untuk kolom tanpa data yang hilang. Saya bisa mundur matriks B kolom terhadap kolom individu dari matriks A, tapi saya punya ribuan regresi yang harus dilakukan, dan ini sangat lambat dan tidak elegan. Fungsi na.exclude tampaknya dirancang untuk kasus ini, tapi saya tidak bisa membuatnya berfungsi. Apa yang saya lakukan salah di sini? Menggunakan R 2.13 pada OSX, jika itu penting.

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
sumber
1
Apa yang Anda maksud dengan "Saya bisa menghitung setiap baris secara individual"?
chl
Maaf, dimaksudkan untuk mengatakan "Saya bisa mundur matriks kolom B terhadap kolom di A secara individual", yang berarti panggilan satu per satu ke lm. Diedit untuk mencerminkan ini.
David Quigley
1
Panggilan satu per satu ke lm / regresi bukanlah cara yang bagus untuk melakukan regresi (mengikuti definisi regresi, yaitu untuk menemukan efek parsial dari masing-masing prediktor pada respons / hasil mengingat keadaan lain variabel)
KarthikS

Jawaban:

23

Sunting: Saya salah mengerti pertanyaan Anda. Ada dua aspek:

a) na.omitdan na.excludekeduanya melakukan penghapusan dengan santai sehubungan dengan prediktor dan kriteria. Mereka hanya berbeda dalam fungsi ekstraktor yang suka residuals()atau fitted()akan NAmengisi output mereka dengan s untuk kasus yang dihilangkan na.exclude, sehingga memiliki output dengan panjang yang sama dengan variabel input.

> 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) Masalah sebenarnya bukan dengan perbedaan antara na.omitdan na.exclude, Anda tampaknya tidak ingin penghapusan dengan santai yang memperhitungkan variabel kriteria, yang keduanya lakukan.

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

Hasil regresi tergantung pada matriks (pseudoinverse dari matriks desain , koefisien ) dan topi matriks , nilai yang dipasang ). Jika Anda tidak ingin penghapusan dengan santai, Anda memerlukan matriks desain berbeda untuk setiap kolom , jadi tidak ada jalan lain untuk menyesuaikan regresi terpisah untuk setiap kriteria. Anda dapat mencoba menghindari overhead dengan melakukan sesuatu di sepanjang baris berikut ini: X β = X + Y H = X X + Y = H Y X YX+=(XX)1XXβ^=X+YH=XX+Y^=HYXYlm()

> 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

Perhatikan bahwa mungkin ada cara yang lebih baik secara numerik untuk menghitung dan , Anda dapat memeriksa dekomposisi- sebagai gantinya. Pendekatan SVD dijelaskan di sini di SE . Saya belum waktunya pendekatan di atas dengan matriks besar terhadap benar-benar menggunakan . H Q R YX+HQRYlm()

caracal
sumber
Itu masuk akal mengingat pemahaman saya tentang bagaimana na.exclude seharusnya bekerja. Namun, jika Anda memanggil> X.both = cbind (X1, X2) dan kemudian> redup (lm (X.both ~ Y, na.action = na.exclude) $ residual) Anda masih mendapatkan 94 residu, bukan 97 dan 97.
David Quigley
Itu peningkatan, tetapi jika Anda melihat residual (lm (X.both ~ Y, na.action = na.exclude)), Anda melihat bahwa setiap kolom memiliki enam nilai yang hilang, meskipun nilai yang hilang di kolom 1 dari X. keduanya berasal dari sampel yang berbeda dari yang ada di kolom 2. Jadi na.exclude mempertahankan bentuk matriks residu, tetapi di bawah tenda R tampaknya hanya mundur dengan nilai yang ada di semua baris X. keduanya. Mungkin ada alasan statistik yang bagus untuk ini, tetapi untuk aplikasi saya ini adalah masalah.
David Quigley
@ David Saya telah salah mengerti pertanyaan Anda. Saya pikir saya sekarang mengerti maksud Anda, dan telah mengedit jawaban saya untuk mengatasinya.
caracal
5

Saya bisa memikirkan dua cara. Salah satunya adalah menggabungkan data menggunakan na.excludedan kemudian memisahkan data lagi:

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)

Cara lain adalah dengan menggunakan dataargumen dan membuat formula.

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)

Jika Anda melakukan banyak regresi, cara pertama harus lebih cepat, karena lebih sedikit sihir latar belakang yang dilakukan. Meskipun jika Anda hanya perlu koefisien dan residu saya sarankan menggunakan lsfit, yang jauh lebih cepat daripada lm. Cara kedua sedikit lebih baik, tetapi pada laptop saya mencoba untuk melakukan ringkasan tentang regresi yang menghasilkan kesalahan. Saya akan mencoba melihat apakah ini bug.

mpiktas
sumber
Terima kasih, tetapi lm (A.ex ~ B.ex) dalam kode Anda cocok dengan 9 poin terhadap A1 (benar) dan 9 poin terhadap A2 (tidak diinginkan). Ada 10 titik yang diukur untuk B1 dan A2; Saya membuang satu titik dalam regresi B1 terhadap A2 karena titik yang sesuai hilang dalam A1. Jika itu hanya cara kerjanya, saya dapat menerimanya, tetapi bukan itu yang saya coba lakukan untuk mendapatkan R.
David Quigley
@ David, oh, sepertinya saya salah mengerti masalah Anda. Saya akan memposting perbaikannya nanti.
mpiktas
1

Contoh berikut menunjukkan cara membuat prediksi dan residu yang sesuai dengan kerangka data asli (menggunakan opsi "na.action = na.exclude" dalam lm () untuk menentukan bahwa NA harus ditempatkan dalam vektor residual dan prediksi di mana bingkai data asli memiliki nilai yang hilang. Hal ini juga menunjukkan bagaimana menentukan apakah prediksi harus mencakup hanya pengamatan di mana variabel penjelas dan variabel dependen lengkap (yaitu, prediksi sampel ketat) atau observasi di mana variabel penjelas lengkap, dan karenanya prediksi Xb dimungkinkan, ( yaitu, termasuk prediksi out-of-sample untuk pengamatan yang memiliki variabel penjelas lengkap tetapi tidak ada variabel dependen).

Saya menggunakan cbind untuk menambahkan prediksi dan variabel sisa ke dataset asli.

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