使用lapply和公式估计具有不同权重的lm。

Question

与如何使用lapply和公式以及在函数中的lapply调用update，为什么不起作用？稍微有些不同：

我试图使用复制权重估计模型。为了得到正确的标准误差，我需要对每个复制权重版本都使用相同的回归模型进行估计。因为我需要估计许多不同的模型，并且不想总是编写单独的循环，所以我尝试编写一个函数，在其中指定回归的数据、回归公式和带有复制权重的数据。虽然当我在函数中显式指定公式并将其放在lapply()命令内而不是作为函数输入时，该函数可以正常工作（下面的函数tryout），但是一旦我将回归公式指定为函数输入（下面的函数tryout2），它就会出现问题。

这里提供一个可重现的例子：

library(tidyverse)
set.seed(123)
lm.dat <- data.frame(id=1:500,
                     x1=sample(1:100, replace=T, size=500),
                     x2=runif(n=500, min=0, max=20)) %>%
  mutate(y=0.2*x1+1.5*x2+rnorm(n=500, mean=0, sd=5))
repweights <- data.frame(id=1:500)
set.seed(123)
for (i in 1:200) {
  repweights[,i+1] <- runif(n=500, min=0, max=10)
  names(repweights)[i+1] <- paste0("hrwgt", i)
}

这两个函数定义如下：

trythis <- function(data, weightsdata, weightsN){
  rep <- as.list(1:weightsN)
  res <- lapply(rep, function(x) lm(data=data, formula=y~x1+x2, weights=weightsdata[,x]))
  return(res)
}
results1 <- trythis(data=lm.dat, weightsdata=repweights[-1], weightsN=200)

trythis2 <- function(LMformula, data, weightsdata, weightsN){
  rep <- as.list(1:weightsN)
  res <- lapply(rep, function(x) lm(data=data, formula=LMformula, weights=weightsdata[,x]))
  return(res)
}

第一个函数可行，但应用第二个函数会导致错误：

trythis2(LMformula = y~x1+x2, data=lm.dat, weightsN=200, weightsdata = repweights[-1])
 Error in eval(extras, data, env) : object 'weightsdata' not found

- andreasf

Answer 1

公式有一个关联的环境，其中可以找到引用变量。在您的情况下，您传递的公式具有调用框架的环境。要访问函数内的变量，您需要将公式分配给本地框架，以便它可以找到正确的变量：

trythis3 <- function(LMformula, data, weightsdata, weightsN){
  rep <- as.list(1:weightsN)
  res <- lapply(rep, function(x) {
    environment(LMformula) <- sys.frames()[[length(sys.frames())]]
    lm(data = data, formula = LMformula, weights = weightsdata[,x])
  })
  return(res)
}

trythis3(LMformula = y~x1+x2, data = lm.dat, weightsN = 200, 
         weightsdata = repweights[-1])

这导致

#> [[1]]
#> 
#> Call:
#> lm(formula = LMformula, data = data, weights = weightsdata[, 
#>     x])
#> 
#> Coefficients:
#> (Intercept)           x1           x2  
#>      1.2932       0.1874       1.4308  
#> 
#> 
#> [[2]]
#> 
#> Call:
#> lm(formula = LMformula, data = data, weights = weightsdata[, 
#>     x])
#> 
#> Coefficients:
#> (Intercept)           x1           x2  
#>      1.2932       0.1874       1.4308  
#> 
#> 
#> [[3]]
#> 
#> Call:
#> lm(formula = LMformula, data = data, weights = weightsdata[, 
#>     x])
#> 
#> Coefficients:
#> (Intercept)           x1           x2  
#>      1.2932       0.1874       1.4308  

...etc

- Allan Cameron