【问题标题】:How to change generic lm() coefficient names in a list of lm objects for use in predict()?如何更改用于预测()的lm对象列表中的通用lm()系数名称?
【发布时间】:2020-06-06 12:11:21
【问题描述】:

我正在尝试编写一个函数,该函数从数据框中创建一个 lm 对象列表,我的数据集中每个唯一日期都有一个不同的 lm。然后,我希望能够将这些 lms 传递给 predict() 并使用新的数据框来生成预测值和置信区间。

数据如下所示:

Date        ppm area
10/18/2019  0   0
10/18/2019  0   0
10/18/2019  0.1 438.9804
10/18/2019  0.1 447.1784
10/18/2019  0.1 443.7794
10/18/2019  1   3232.2088
10/18/2019  1   3206.6672
10/18/2019  1   3206.232
10/24/2019  0   0
10/24/2019  0   15.98
10/24/2019  0   0
10/24/2019  0   0
10/24/2019  0.1 379.387
10/24/2019  0.1 325.5268
10/24/2019  0.1 325.8126
10/24/2019  0.1 310.5972
10/24/2019  1   3259.366
10/24/2019  1   3218.0836
10/24/2019  1   3192.7076

第一部分似乎很简单——编写一个为每个日期创建不同 lm 的函数:

standard.lm= function(standards,
                          date_field = "date",
                          peak_field,
                          std_field,
                          peak_field2 = NA){
  library(tidyverse)
  library(broom)


  y = standards %>% nest(-date_field) %>%
    mutate(fit = map(data, ~lm(.[[std_field]] ~ .[[peak_field]], data = .)))

    return(y)  }

然后我可以运行命令:

test = standard.lm(standard_data, std_field = "std.ppm", peak_field = "area")

这可以很好地为每个日期生成lms,但问题是系数被命名为.[[peak_field]]而不是"area"

这给我带来了一个问题,因为我想将这些 lm 对象传递给 predict() 以从区域测量中预测 ppm 值。我在下一个数据表中的列将命名为area,我无法将其重命名为.[[peak_field]]。我尝试这样的事情,我得到一个错误:

a = c(1300.1, 1400.3, 1500.9)
df = data.frame(area = a)
df$std.ppm = predict(test$fit[[1]], newdata = df)

$<-.data.frame(*tmp*, std.ppm, value = c(1 = -0.00299110569401364, : 替换有 8 行,数据有 3 另外:警告信息: 'newdata' 有 3 行,但找到的变量有 8 行

发生这种情况是因为 predict() 正在寻找名为 .[[peak_field]] 的列而不是识别区域,并且正在预测原始输入 lm 数据的值而不是我希望它预测的数据。

所以基本上我正在寻找解决这个问题的解决方案。最好的解决方案是让我在第一个函数中最初创建 lm 对象时指定系数名称,但我可以接受允许我指定在 predict() 中使用哪一列的东西

【问题讨论】:

标签: r tidyverse lm predict


【解决方案1】:

您可以尝试使用定义的 y 和 x 变量在函数中创建公式:

standard.lm= function(standards,date_field = "Date",
                      peak_field,std_field,peak_field2 = NA){
  lm_form = as.formula(paste(std_field,"~",peak_field))
  #another away
  #lm_form = substitute(y~x,list(y=as.name(std_field),x=as.name(peak_field)))
  y = standards %>% nest(data=-one_of(date_field)) %>%
    mutate(fit = map(data, ~lm(lm_form, data = .)))

    return(y)  }

我们对其进行测试:

standard_data = structure(list(Date = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("10/18/2019", 
"10/24/2019"), class = "factor"), ppm = c(0, 0, 0.1, 0.1, 0.1, 
1, 1, 1, 0, 0, 0, 0, 0.1, 0.1, 0.1, 0.1, 1, 1, 1), area = c(0, 
0, 438.9804, 447.1784, 443.7794, 3232.2088, 3206.6672, 3206.232, 
0, 15.98, 0, 0, 379.387, 325.5268, 325.8126, 310.5972, 3259.366, 
3218.0836, 3192.7076)), class = "data.frame", row.names = c(NA, 
-19L))

mdl = standard.lm(standard_data, std_field = "ppm", peak_field = "area")

predict(mdl$fit[[1]], data.frame(area=c(1300.1,1400.3)))
        1         2 
0.3897161 0.4215205 

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2021-11-06
    • 2014-05-06
    • 2012-11-06
    • 1970-01-01
    • 2016-09-20
    • 1970-01-01
    • 1970-01-01
    • 2017-11-26
    相关资源
    最近更新 更多