如何正确地将由`lm`拟合的线性模型通过`dput`函数保存到ASCII文件中，并在后续重新创建该模型？

Question

我想将lm对象保存到文件中，并在另一个程序中重新加载它。 我知道可以使用saveRDS/readRDS通过写入/读取二进制文件来实现，但我希望有一个ASCII格式的文件。更一般地说，我想知道为什么我通常用于读取dput输出的习惯用法没有按照我的预期那样工作。

以下是制作简单拟合的示例以及成功和失败重新创建模型的示例：

dat_train <- data.frame(x=1:4, z=c(1, 2.1, 2.9, 4))
fit <- lm(z ~ x, dat_train)
rm(dat_train) # Just to make sure fit is not dependent upon `dat_train existence`

dat_score <- data.frame(x=c(1.5, 3.5))

## This works (of course)
predict(fit, dat_score)
#    1    2 
# 1.52 3.48

二进制文件保存可行：

## https://dev59.com/QW435IYBdhLWcg3w9FBY
saveRDS(fit, "model.RDS")
fit2 <- readRDS("model.RDS")
predict(fit2, dat_score)
#    1    2 
# 1.52 3.48

这样做（将dput应用于R会话而不是文件）：

fit2 <- eval(dput(fit))
predict(fit2, dat_score)
#    1    2 
# 1.52 3.48

但如果我将文件持久化到磁盘，我就不知道如何恢复正常状态：

dput(fit, file = "model.R")
fit3 <- source("model.R")$value

# Error in is.data.frame(data): object 'dat_train' not found

predict(fit3, dat_score)
# Error in predict(fit3, dat_score): object 'fit3' not found

试图显式使用 eval 也不起作用：

## https://dev59.com/SGox5IYBdhLWcg3wpV3F
dput(fit, file="model.R")
fit4 <- eval(parse(text=paste(readLines("model.R"), collapse=" ")))

# Error in is.data.frame(data): object 'dat_train' not found

predict(fit4, dat_score)
# Error in predict(fit4, dat_score): object 'fit4' not found

以上两种情况，我期望fit3和fit4都能正常工作，但它们无法重新编译为可用于predict()的lm对象。

请问有谁可以建议我如何将模型保存到一个具有类似ASCII结构的structure(...)文件中，并将其作为lm对象重新读取以便在predict()中使用？另外，请问我的当前方法为什么不起作用？

- mpettis

---

我尝试使用dget重新读取文件，但是我得到了错误信息“Error in eval(expr, envir, enclos) : object 'z' not found.”。dput/dget组合应该可以工作，但在这种情况下似乎并不起作用。 - Dave2e

Answer 1

步骤1：

您需要控制去解析选项：

dput(fit, control = c("quoteExpressions", "showAttributes"), file = "model.R") 

您可以在?deparseOpts中阅读有关所有可能选项的更多信息。

---

"quoteExpressions"将所有调用/表达式/语言都用quote包装起来，以便在您以后重新解析时不对它们进行评估。注意：

• source正在执行解析；
• 您拟合的“lm”对象中的call字段是一个调用：

fit$call
# lm(formula = z ~ x, data = dat_train)

如果没有使用“quoteExpressions”选项，R将尝试在解析期间评估lm调用。如果进行评估，则会拟合线性模型，并且R将尝试找到dat_train，而在新的R会话中该对象将不存在。

---

“showAttributes”是另一个必需的选项，因为“lm”对象具有类属性。您肯定不希望放弃所有类属性并仅导出一个普通的“list”对象，对吗？此外，“lm”对象中的许多元素，如model（模型框架）、qr（紧凑QR矩阵）和terms（术语信息）等都具有属性。您想要保留它们所有。

---

如果您不设置control，则默认设置为：

control = c("keepNA", "keepInteger", "showAttributes")

将被使用。如您所见，没有“quoteExpressions”，因此您会遇到麻烦。

您还可以指定“keepInteger”和“keepNA”，但我不认为需要“lm”对象。

------

第二步：

上述步骤将使source正常工作。您可以恢复您的模型：

fit1 <- source("model.R")$value

然而，它尚未准备好为像summary和predict这样的通用函数工作。 为什么？

关键问题在于fit1中的terms对象实际上不是一个真正的“terms”对象，而只是一个公式（它甚至不是一个公式，而只是一个没有“formula”类的“语言”对象！）。 只需比较fit$terms和fit1$terms，就能看出差异。 不要惊讶；我们之前设置了“quoteExpressions”。 虽然这绝对有助于防止call的评估，但它对terms产生了副作用。 因此，我们需要尽可能重建terms。

幸运的是，做以下操作已经足够：

fit1$terms <- terms.formula(fit1$terms)

尽管这仍然无法恢复fit$terms中的所有信息（例如变量类别丢失），但它很容易成为一个有效的“terms”对象。

为什么“terms”对象很关键？因为所有通用函数都依赖于它。您可能不需要更多了解，因为这确实很技术性，所以我会在这里停下。

完成此操作后，我们就可以成功使用predict（和summary）：

predict(fit1)  ## no `newdata` given, using model frame `fit1$model`
#   1    2    3    4 
#1.03 2.01 2.99 3.97 

predict(fit1, dat_score)  ## with `newdata`
#   1    2 
#1.52 3.48 

-------

结论：

虽然我向您展示了如何使事情正常运作，但我并不推荐一般情况下这样做。例如，在将模型拟合到大型数据集时，“lm”对象将非常大，例如，残差、拟合值是长向量，而qr和model是巨大的矩阵/数据框。所以，请考虑一下。

- Zheyuan Li

Answer 2

这是一个重要的更新！

如前面所提到的，最具挑战性的部分是尽可能恢复$terms。建议使用terms.formula的方法适用于OP的示例，但不适用于以下使用bs()和poly()的示例：

dat <- data.frame(x1 = runif(20), x2 = runif(20), x3 = runif(20), y = rnorm(20))
library(splines)
fit <- lm(y ~ bs(x1, df = 3) + poly(x2, degree = 3) + x3, data = dat)
rm(dat)

如果我们遵循前面的答案：

dput(fit, control = c("quoteExpressions", "showAttributes"), file = "model.R") 
fit1 <- source("model.R")$value
fit1$terms <- terms.formula(fit1$terms)

我们将看到summary.lm和anova.lm的工作是正确的，但predict.lm不正确：

predict(fit1, newdata = data.frame(x1 = 0.5, x2 = 0.5, x3 = 0.5))

bs(x1, df = 3) 函数未找到。

这是因为 $terms 的 ".Environment" 属性缺失。我们需要

environment(fit1$terms) <- .GlobalEnv

现在再次运行上面的predict，我们会看到一个不同的错误：

在poly(x2, degree = 3)中出现错误：

'degree'必须小于唯一点的数量

这是因为我们缺少对bs()和poly()进行安全/正确预测所需的"predvars"属性。

解决方案是，我们需要额外地这些特殊属性：

dput(attr(fit$terms, "predvars"), control = "quoteExpressions", file = "predvars.R")

然后阅读并添加它。

attr(fit1$terms, "predvars") <- source("predvars.R")$value

现在运行predict正常工作。

请注意，$terms的"dataClass"属性也缺失了，但这似乎对任何通用函数都没有造成问题。

- Zheyuan Li