2013-07-11 7 views
1

Как функция stepAIC() от пакета MASS имеет проблемы при использовании внутри функции, я использую ее с do.call() (описано here). Моя проблема звучит очень просто, но я не нашел для нее решения: когда я использую do.call() для модели lm() с несколькими растровыми слоями, все слои сохраняются в пределах модели. Если я хочу напечатать summary() модели, он записывает все слои в выходном файле и становится очень запутанным. Как получить «нормальный» вывод summary, как я мог бы получить без использования do.call?R: Преобразование do.call() - summary in summary

Вот краткий пример:

Создать список растровых слоев:

xz.list <- lapply(1:5,function(x){ 
    r1 <- raster(ncol=3, nrow=3) 
    values(r1) <- 1:ncell(r1) 
    r1 
}) 

Преобразовать их в data.frame:

xz<-getValues(stack(xz.list)) 

xz <- as.data.frame(xz) 

Использование do.call для lm модели:

fit1<-do.call("lm", list(xz[,1] ~ . , data = xz)) 

summary() результат выглядит следующим образом:

summary(fit1) 

Call: 
lm(formula = xz[, 1] ~ ., data = structure(list(layer.1 = 1:9, 
    layer.2 = 1:9, layer.3 = 1:9, layer.4 = 1:9, layer.5 = 1:9), .Names = c("layer.1", 
"layer.2", "layer.3", "layer.4", "layer.5"), row.names = c(NA, 
-9L), class = "data.frame")) 

Residuals: 
     Min   1Q  Median   3Q  Max 
-9.006e-16 -2.472e-16 -2.031e-16 -1.370e-16 1.724e-15 

Coefficients: (4 not defined because of singularities) 
      Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.184e-15 5.784e-16 2.047e+00 0.0798 . 
layer.1  1.000e+00 1.028e-16 9.729e+15 <2e-16 *** 
layer.2   NA   NA  NA  NA  
layer.3   NA   NA  NA  NA  
layer.4   NA   NA  NA  NA  
layer.5   NA   NA  NA  NA  
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 7.962e-16 on 7 degrees of freedom 
Multiple R-squared:  1, Adjusted R-squared:  1 
F-statistic: 9.465e+31 on 1 and 7 DF, p-value: < 2.2e-16 

Это не выглядит плохо в этом маленьком примере, но это становится беспорядок, когда вы используете 10 или более raster слоев около 32к значений каждого. Так что я хотел бы сделать вывод выглядеть, как я бы просто использовать summary(lm) функцию без do.call:

fit<-lm(xz[,1] ~ . , data=xz) 
summary(fit) 

Call: 
lm(formula = xz[, 1] ~ ., data = xz) 

Residuals: 
     Min   1Q  Median   3Q  Max 
-9.006e-16 -2.472e-16 -2.031e-16 -1.370e-16 1.724e-15 

Coefficients: (4 not defined because of singularities) 
      Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.184e-15 5.784e-16 2.047e+00 0.0798 . 
layer.1  1.000e+00 1.028e-16 9.729e+15 <2e-16 *** 
layer.2   NA   NA  NA  NA  
layer.3   NA   NA  NA  NA  
layer.4   NA   NA  NA  NA  
layer.5   NA   NA  NA  NA  
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 7.962e-16 on 7 degrees of freedom 
Multiple R-squared:  1, Adjusted R-squared:  1 
F-statistic: 9.465e+31 on 1 and 7 DF, p-value: < 2.2e-16 

ответ

1

Вы можете переопределить вашу lm функцию следующим образом:

lm <- function(form, ...) { fm <- stats::lm(form,...); 
          fm$call <- form; fm } 

тестирования:

fit2<-do.call("lm", list(xz[,1] ~ . , data = xz)) 

summary(fit2) 

Call: 
xz[, 1] ~ . 

Residuals: 
     Min   1Q  Median   3Q  Max 
-9.006e-16 -2.472e-16 -2.031e-16 -1.370e-16 1.724e-15 

Coefficients: (4 not defined because of singularities) 
      Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.184e-15 5.784e-16 2.047e+00 0.0798 . 
layer.1  1.000e+00 1.028e-16 9.729e+15 <2e-16 *** 
layer.2   NA   NA  NA  NA  
layer.3   NA   NA  NA  NA  
layer.4   NA   NA  NA  NA  
layer.5   NA   NA  NA  NA  
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 7.962e-16 on 7 degrees of freedom 
Multiple R-squared:  1, Adjusted R-squared:  1 
F-statistic: 9.465e+31 on 1 and 7 DF, p-value: < 2.2e-16