有没有办法通过从 R 中的模型中获取参数来构建真正的回归方程?

Is there any way to construct real regression equation by taking parameters from models in R?

数据是:

d <- data.frame(x = rnorm(100, 0, 1),
            y = rnorm(100, 0, 1),
            z = rnorm(100, 0, 1))

适合 5 个模型的函数

library(splines)
func <-function(d){
  fit1 <- lm( y~ x + z, data = d)
  fit2 <- lm( y~x + I(z^2), data = d)
  fit3 <- lm( y~poly(x,3) + z, data = d)
  fit4 <- lm( y~ns(x, 3) + z, data = d)
  l <- list(fit1, fit2, fit3, fit4)
  names(l) <- paste0("fit", 1:4) 
  return(l)
}

mods <- func(d) 
mods[[1]]

stargazer(mods, type="text)

我想通过从拟合模型中获取参数并在可能的情况下在 R 中自动查找变量,以每个模型的实数格式构建实数回归方程。例如:对于 fit1 模型,截距 = -0.20612,x = 0.17443,x = 0.03203。然后等式将是这样的:y = -0.206 + 0.174x + 0.032z 等,并想在 table 中列出所有模型的这些等式以及非常常见的有用统计数据,如 R2、P 值、adj.R2、观察等 stargazer 没有向我展示我想要的输出。所以我想确定是否有任何方法可以在 R 中执行此操作而无需在 excel?

中手动执行

提前致谢!

假设 mods 显示在末尾的注释中,并且我们想要的是公式的文本表示的字符向量,其中系数被替换,我们有以下内容。

fit2text 函数采用拟合对象并输出带有公式文本表示的字符串。 round 参数给出了结果中系数四舍五入的位数。 rmI 参数,如果为 TRUE,则删除任何 I(...) 并只保留 ... 内部假设,为了便于实现,内部表达式不包含任何括号。如果为 FALSE,则 I 不会被删除。

其他统计数据可以从summary(mods[[1]])broom::glance(mods[[1]])

中提取
fit2text <- function(fit, round = 2, rmI = TRUE) {
  fo <- formula(fit)  
  resp <- all.vars(fo)[1]
  co <- round(coef(fit), round)
  labs <- c(if (terms(fit, "intercept") == 1) "", labels(fit))
  p <- gsub("\+ *-", "- ", paste(resp, "~ ", paste(paste(co, labs), collapse = " + ")))
  p2 <- if (rmI) gsub("I\(([^)]+)\)", "\1", p) else p
  gsub(" +", " ", p2)
}
sapply(mods, fit2text)

给予:

                                                           fit1 
                                  "y ~ -0.11 - 0.05 x + 0.03 z" 
                                                           fit2 
                                "y ~ -0.07 - 0.05 x - 0.04 z^2" 
                                                           fit3 
"y ~ -0.11 - 0.43 poly(x, 3) - 1.05 z + 0.27 + 0.04 poly(x, 3)" 
                                                           fit4 
    "y ~ -0.55 + 0.23 ns(x, 3) + 0.79 z - 0.25 + 0.04 ns(x, 3)" 

备注

问题中的代码不可重现,因为缺少库调用,它使用没有 set.seed 的随机数,并且代码中还有一些其他错误。为清楚起见,我们提供了以下可重现的代码,用于为上述答案提供输入。

library(splines)
set.seed(123)

d <- data.frame(x = rnorm(100, 0, 1),
            y = rnorm(100, 0, 1),
            z = rnorm(100, 0, 1))

# function to fit 5 models
func <-function(d){
  fit1 <- lm( y~ x + z, data = d)
  fit2 <- lm( y~x + I(z^2), data = d)

  fit3 <- lm( y~poly(x,3) + z, data = d)
  fit4 <- lm( y~ns(x, 3) + z, data = d)
  l <- list(fit1, fit2, fit3, fit4)
  names(l) <- paste0("fit", 1:4) 
  return(l)
}

mods <- func(d) 

问题是您的每个模型都不是表格数据的理想选择,例如拟合 3 returns 4 个估计而拟合 1 returns 仅 3

如果您对列表感到满意,我建议它们是存储此类信息的好方法

library(broom)
library(tidyverse)
library(splines)

d <- data.frame(x = rnorm(100, 0, 1),
                y = rnorm(100, 0, 1),
                z = rnorm(100, 0, 1))

func <-function(d){
  fit1 <- lm( y~ x + z, data = d)
  fit2 <- lm( y~x + I(z^2), data = d)
  fit3 <- lm( y~poly(x,3) + z, data = d)
  fit4 <- lm( y~ns(x, 3) + z, data = d)
  l <- list(fit1, fit2, fit3, fit4)
  names(l) <- paste0("fit", 1:4) 
  return(l)
}

mods <- func(d) 

list_representation<- map(mods,tidy)  

我们可以map通过mods使用@J.R.的函数broom::glance到模型R2, P-value,和 adj.R2.

library(purrr)
library(broom)
map_dfr(mods, 
        function(x) data.frame('Eq'=regEq(lmObj = x, dig = 3), broom::glance(x), stringsAsFactors = FALSE), 
        .id='Model') 

  Model                                                                              Eq    r.squared adj.r.squared    sigma  statistic   p.value df    logLik      AIC      BIC
1  fit1                                                   y = 0.091 - 0.022*x - 0.027*z 0.0012601436   -0.01933243 1.028408 0.06119408 0.9406769  3 -143.1721 294.3441 304.7648
2  fit2                                              y = 0.093 - 0.022*x - 0.003*I(z^2) 0.0006154188   -0.01999045 1.028740 0.02986619 0.9705843  3 -143.2043 294.4087 304.8294
3  fit3 y = 0.093 - 0.248*poly(x, 3)1 - 0.186*poly(x, 3)2 - 0.581*poly(x, 3)3 - 0.031*z 0.0048717358   -0.03702840 1.037296 0.11627016 0.9764662  5 -142.9909 297.9819 313.6129
4  fit4        y = 0.201 + 0.08*ns(x, 3)1 - 0.385*ns(x, 3)2 - 0.281*ns(x, 3)3 - 0.031*z 0.0032813558   -0.03868575 1.038125 0.07818877 0.9887911  5 -143.0708 298.1416 313.7726
  deviance df.residual
1 102.5894          97
2 102.6556          97
3 102.2184          95
4 102.3818          95