5

我正在寻找有关如何在所有出现的独立/解释变量都是 NA（即x3以下）时处理线性回归中的 NA 的建议。

我知道显而易见的解决方案是从模型中排除有问题的自变量/解释变量，但我正在遍历多个区域，并且不希望每个区域都有不同的功能形式。

下面是一些示例数据：

set.seed(23409)
n <- 100

time <- seq(1,n, 1)
x1 <- cumsum(runif(n))           
y  <- .8*x1 + rnorm(n, mean=0, sd=2)
x2 <- seq(1,n, 1)       
x3 <- rep(NA, n)            
df <- data.frame(y=y, time=time, x1=x1, x2=x2, x3=x3)

# Quick plot of data
library(ggplot2)
library(reshape2)
df.melt <-melt(df, id=c("time"))

p <- ggplot(df.melt, aes(x=time, y=value)) + 
  geom_line() + facet_grid(variable ~ .)
p

我已阅读文档lm并尝试了各种na.action设置但均未成功：

lm(y~x1+x2+x3, data=df, singular.ok=TRUE)

lm(y~x1+x2+x3, data=df, na.action=na.omit)
lm(y~x1+x2+x3, data=df, na.action=na.exclude)

lm(y~x1+x2+x3, data=df, singular.ok=TRUE, na.exclude=na.omit)
lm(y~x1+x2+x3, data=df, singular.ok=TRUE, na.exclude=na.exclude)

有没有办法让 lm 无错误地运行，并简单地返回一个系数，以反映相关变量缺乏解释能力（即零或 NA）的解释性反映？

4

2 回答 2

3

您将无法包含包含所有NA值的列。它会做一些奇怪的事情model.matrix

 x1 <- 1:5
 x2 <- rep(NA,5)

 model.matrix(~x1+x2) 
     (Intercept) x1 x2TRUE
attr(,"assign")
[1] 0 1 2
attr(,"contrasts")
attr(,"contrasts")$x2
[1] "contr.treatment"

因此，您的替代方法是根据数据以编程方式创建模型公式。

就像是...

make_formula <- function(variables, data, response = 'y'){
   if(missing(data)){stop('data not specified')}
   using <-  Filter(variables,f= function(i) !all(is.na(data[[i]])))

   deparse(reformulate(using, response))
 }

 variables <- c('x1','x2','x3')

make_formula(variables, data =df)

[1] "y ~ x1 + x2"

我曾经deparse返回一个字符串，这样environment在函数中创建公式就不会出现问题。lm可以愉快地接受一个有效公式的字符串。

于 2013-03-14T04:00:43.300 回答
3

这是一个想法：

set.seed(23409)
n <- 100

time <- seq(1,n, 1)
x1 <- cumsum(runif(n))           
y  <- .8*x1 + rnorm(n, mean=0, sd=2)
x2 <- seq(1,n, 1)       
x3 <- rep(NA, n)            
df <- data.frame(y=y, time=time, x1=x1, x2=x2, x3=x3)

replaceNA<-function(x){
  if(all(is.na(x))){
    rep(0,length(x)) 
  } else x

} 

lm(y~x1+x2+x3, data= data.frame(lapply(df,replaceNA)))
Call:
lm(formula = y ~ x1 + x2 + x3, data = data.frame(lapply(df, replaceNA)))

Coefficients:
(Intercept)           x1           x2           x3  
    0.05467      1.01133     -0.10613           NA  

lm(y~x1+x2, data=df)
Call:
lm(formula = y ~ x1 + x2, data = df)

Coefficients:
(Intercept)           x1           x2  
    0.05467      1.01133     -0.10613

因此，您将仅包含NA's 的变量替换为仅包含 0 的变量。你得到系数值 NA，但模型拟合的所有相关部分都是相同的（期望 qr 分解，但如果需要相关信息，可以轻松修改）。请注意，组件summary(fit)$alias（请参阅 参考资料?alias）可能有用。

这似乎与您的另一个问题有关：Replace lm coefficients in [r]

于 2013-03-14T19:53:23.113 回答