r中带有公式的箱线图
boxplots with formulas in r
我不明白公式在 R 中是如何工作的，什么是公式。
因此，如果我有一个包含每月时间序列的向量，我如何创建它的箱线图，其中数据按季节性逻辑划分？我想要12盒，一个月一盒。
在存储级别，公式是parse tree。解析树编码对 `~`() 函数的调用，使用一个或两个参数。对于单边公式，它采用一个参数来表示 RHS，而对于双面公式，它采用两个参数来表示公式的 LHS 和 RHS。
应该提到的是，对嵌入在公式的解析树存储表示中的 `~`() 的调用实际上没有任何意义。通常，`~`() 函数实际上不做任何事情，除了允许创建公式对象，无论是显式（例如 `~`(a+b,c/d)）还是使用 R 语言提供的语法糖功能（例如a+b~c/d)。在公式的解析树存储表示的顶层编码中使用 `~`() 函数是一个相当随意且无关紧要的实现细节。我稍后会对此进行扩展。
R语言使得将解析树分解成递归列表结构成为可能，这可以帮助我们检查和理解这些解析树的结构。
我写了一个简短的递归函数可以做到这一点：
ptunwrap <- function(x) if (typeof(x)=='language') lapply(as.list(x),ptunwrap) else x;
那么让我们看一个例子：
f1 <- a+b~c/d;
f1;
## a + b ~ c/d
ptunwrap(f1);
## [[1]]
## `~`
##
## [[2]]
## [[2]][[1]]
## `+`
##
## [[2]][[2]]
## a
##
## [[2]][[3]]
## b
##
##
## [[3]]
## [[3]][[1]]
## `/`
##
## [[3]][[2]]
## c
##
## [[3]][[3]]
## d
##
当分析树包含函数调用时，函数调用表示为递归列表结构中的单个列表节点。函数的符号作为该列表的第一个列表组件嵌入，其参数作为列表的后续列表组件嵌入。
上面的解析树中有3个函数调用。如前所述，顶级列表表示对 `~`() 函数的调用。第二个顶级列表组件进一步分支到另一个列表，该列表包含对 `+`() 函数的调用，该函数本身有两个参数，符号 a 和 b。第三个组件类似，表示对 `/`() 函数的调用并再次采用两个参数，符号 c 和 d.
重要的是要理解，虽然解析树始终表示语法上有效的 R 代码，并且具有能力 被评估以产生单个结果值，但它不是 需要 来评估解析树。完全有可能创建一个解析树而不对其求值。
那么，创建一个您从不求值的解析树的目的是什么？在 R 中，通常这样做是为了使用方便的语法将某些信息传递给函数。
作为一个随机示例，data.table 包允许使用 := 运算符使用以下语法糖将列添加到现有 data.table：
library(data.table);
dt <- data.table(a=1:3);
dt[,b:=a*2L];
dt;
## a b
## 1: 1 2
## 2: 2 4
## 3: 3 6
在内部，data.table 使用非标准参数评估来检索 `[.data.table`() 函数的第二个参数（技术上在函数定义中的第三个；语法糖形式中的第二个）的解析树，通常被称为“j 参数”，因为参数名称是 j。如果需要，您实际上可以直接在 R 中检查源代码本身。这是最相关的代码片段：
data.table:::`[.data.table`;
## function (x, i, j, by, keyby, with = TRUE, nomatch = getOption("datatable.nomatch"),
## mult = "all", roll = FALSE, rollends = if (roll == "nearest") c(TRUE,
## TRUE) else if (roll >= 0) c(FALSE, TRUE) else c(TRUE,
## FALSE), which = FALSE, .SDcols, verbose = getOption("datatable.verbose"),
## allow.cartesian = getOption("datatable.allow.cartesian"),
## drop = NULL, on = NULL)
## {
##
## ... snip ...
##
## if (!missing(j)) {
## jsub = substitute(j)
## if (is.call(jsub))
## jsub = construct(deconstruct_and_eval(replace_dot(jsub),
## parent.frame(), parent.frame()))
## }
##
## ... snip ...
##
我们可以看到他们正在使用 substitute(j) 来检索 j 参数的解析树。对于上面的演示，这是他们将得到的：
ptunwrap(substitute(b:=a*2L));
## [[1]]
## `:=`
##
## [[2]]
## b
##
## [[3]]
## [[3]][[1]]
## `*`
##
## [[3]][[2]]
## a
##
## [[3]][[3]]
## [1] 2
##
稍后在代码中，他们测试顶级函数符号是否为 :=，这是支持向 data.table。如果是这样，他们将测试 LHS 是否由单个裸词组成，该裸词被视为要添加（或修改或删除）的列的名称。请注意，在这种情况下，实际上 不可能 他们实际评估解析树的 LHS，因为它包含 [=] 中尚不存在的符号141=]。但是，RHS 最终会被评估以生成要添加到新名称下的 data.table 的列向量。
所以应该清楚的是，公式可以在 R 中的各种上下文中使用，并且它们并不总是被评估。有时会简单地检查解析树以检索从调用者传递给被调用者的信息。即使在 被 评估的上下文中，有时也只有 LHS 或 RHS（或两者）会被独立评估，忽略当时嵌入在解析树中的顶级函数符号它已创建。
转到boxplot()函数，让我们看一下formula参数的文档：
a formula, such as y ~ grp, where y is a numeric vector of data values to be split into groups according to the grouping variable grp (usually a factor).
在这种情况下，公式的两边最终被独立评估，LHS 提供数据向量，RHS 提供分组定义。
证明这一点的一个好方法如下：
boxplot(1:9~1:9%%3L);
注意公式的两边是如何由文字表达式组成的：
1:9;
## [1] 1 2 3 4 5 6 7 8 9
1:9%%3L;
## [1] 1 2 0 1 2 0 1 2 0
在内部，boxplot() 必须独立评估公式的每一边以生成数据和分组向量，几乎就像您将两个表达式作为单独的参数传递一样。
那么，让我们创建一个月度时间序列箱线图的简单演示：
N <- 36L; df <- data.frame(date=seq(as.Date('2016-01-01'),by='month',len=N),y=rnorm(N));
df;
## date y
## 1 2016-01-01 -1.56004488
## 2 2016-02-01 0.65699747
## 3 2016-03-01 0.05729631
## 4 2016-04-01 -0.02092276
## 5 2016-05-01 0.46673530
## 6 2016-06-01 -0.18652580
## 7 2016-07-01 0.06228650
## 8 2016-08-01 1.54452267
## 9 2016-09-01 1.06643594
## 10 2016-10-01 -1.51178160
## 11 2016-11-01 0.82904673
## 12 2016-12-01 0.37667201
## 13 2017-01-01 -0.10135801
## 14 2017-02-01 0.94692462
## 15 2017-03-01 -1.60781946
## 16 2017-04-01 0.47189753
## 17 2017-05-01 -1.32869317
## 18 2017-06-01 -0.49821455
## 19 2017-07-01 0.54474606
## 20 2017-08-01 0.47565264
## 21 2017-09-01 -0.97494730
## 22 2017-10-01 -1.22781588
## 23 2017-11-01 -0.34919086
## 24 2017-12-01 -0.78153843
## 25 2018-01-01 -0.59355220
## 26 2018-02-01 -2.58287605
## 27 2018-03-01 1.42148186
## 28 2018-04-01 -1.01278176
## 29 2018-05-01 -0.80961662
## 30 2018-06-01 0.19793126
## 31 2018-07-01 -1.03072915
## 32 2018-08-01 -0.87896416
## 33 2018-09-01 -2.36216655
## 34 2018-10-01 1.82708221
## 35 2018-11-01 0.05579195
## 36 2018-12-01 1.28612246
boxplot(y~months(date),df);
需要的话可以研究下源码，里面需要追溯S3的查找过程：
boxplot;
## function (x, ...)
## UseMethod("boxplot")
## <bytecode: 0x600b50760>
## <environment: namespace:graphics>
graphics:::boxplot.formula;
## function (formula, data = NULL, ..., subset, na.action = NULL)
## {
## if (missing(formula) || (length(formula) != 3L))
## stop("'formula' missing or incorrect")
## m <- match.call(expand.dots = FALSE)
## if (is.matrix(eval(m$data, parent.frame())))
## m$data <- as.data.frame(data)
## m$... <- NULL
## m$na.action <- na.action
## m[[1L]] <- quote(stats::model.frame)
## mf <- eval(m, parent.frame())
## response <- attr(attr(mf, "terms"), "response")
## boxplot(split(mf[[response]], mf[-response]), ...)
## }
## <bytecode: 0x6035c67f8>
## <environment: namespace:graphics>
它几乎是令人恼火的迂回和复杂，但 graphics:::boxplot.formula() 有效地检索（通过 match.call()) the parse tree that caused its own invocation, massages it a little, most notably replacing its own function symbol boxplot.formula with stats::model.frame, and then evaluates that new parse tree, thereby calling stats::model.frame()。该函数本身非常复杂，并涉及进一步的 S3 查找，但这是最相关的代码：
model.frame;
## function (formula, ...)
## UseMethod("model.frame")
## <bytecode: 0x601464b18>
## <environment: namespace:stats>
model.frame.default;
## function (formula, data = NULL, subset = NULL, na.action = na.fail,
## drop.unused.levels = FALSE, xlev = NULL, ...)
## {
##
## ... snip ...
##
## if (!inherits(formula, "terms"))
## formula <- terms(formula, data = data)
## env <- environment(formula)
## rownames <- .row_names_info(data, 0L)
## vars <- attr(formula, "variables")
## predvars <- attr(formula, "predvars")
## if (is.null(predvars))
## predvars <- vars
## varnames <- sapply(vars, function(x) paste(deparse(x, width.cutoff = 500),
## collapse = " "))[-1L]
## variables <- eval(predvars, data, env)
##
## ... snip ...
##
因此，最终，它从公式对象中检索各个表达式，并使用 eval() 以给定的 data.frame 和公式的闭包环境作为上下文对它们进行评估，从而给出结果向量:
attr(terms(y~months(date),data=df),'variables');
## list(y, months(date))
eval(attr(terms(y~months(date),data=df),'variables'),df);
## [[1]]
## [1] -1.56004488 0.65699747 0.05729631 -0.02092276 0.46673530 -0.18652580 0.06228650 1.54452267 1.06643594 -1.51178160 0.82904673 0.37667201 -0.10135801 0.94692462 -1.60781946 0.47189753 -1.32869317 -0.49821455 0.54474606
## [20] 0.47565264 -0.97494730 -1.22781588 -0.34919086 -0.78153843 -0.59355220 -2.58287605 1.42148186 -1.01278176 -0.80961662 0.19793126 -1.03072915 -0.87896416 -2.36216655 1.82708221 0.05579195 1.28612246
##
## [[2]]
## [1] "January" "February" "March" "April" "May" "June" "July" "August" "September" "October" "November" "December" "January" "February" "March" "April" "May" "June" "July"
## [20] "August" "September" "October" "November" "December" "January" "February" "March" "April" "May" "June" "July" "August" "September" "October" "November" "December"
##
为了重申之前提出的观点，请注意 `~`() 函数在评估过程中是如何无处可寻的。这是 R 公式的任意实现细节，它们将 `~`() 函数编码为公式对象的解析树存储表示中的顶级函数符号。公式边的实际求值（如果发生）不涉及该函数的求值。
最后，让我们考虑一下如果您实际评估包含公式存储表示的整个解析树会发生什么。回想一下，`~`() 函数只是根据其参数创建一个公式。因此，评估公式具有吐出与刚刚评估的完全相同的公式的有趣效果：
f1;
## a + b ~ c/d
eval(f1);
## a + b ~ c/d
eval(eval(f1));
## a + b ~ c/d
我已经写了几个关于解析树和公式的其他答案，您可能会感兴趣：
-- 重点介绍解析树涉及的数据类型 -- 关注公式和函数的闭包环境
我不明白公式在 R 中是如何工作的，什么是公式。 因此，如果我有一个包含每月时间序列的向量，我如何创建它的箱线图，其中数据按季节性逻辑划分？我想要12盒，一个月一盒。
在存储级别，公式是parse tree。解析树编码对 `~`() 函数的调用，使用一个或两个参数。对于单边公式，它采用一个参数来表示 RHS，而对于双面公式，它采用两个参数来表示公式的 LHS 和 RHS。
应该提到的是，对嵌入在公式的解析树存储表示中的 `~`() 的调用实际上没有任何意义。通常，`~`() 函数实际上不做任何事情，除了允许创建公式对象，无论是显式（例如 `~`(a+b,c/d)）还是使用 R 语言提供的语法糖功能（例如a+b~c/d)。在公式的解析树存储表示的顶层编码中使用 `~`() 函数是一个相当随意且无关紧要的实现细节。我稍后会对此进行扩展。
R语言使得将解析树分解成递归列表结构成为可能，这可以帮助我们检查和理解这些解析树的结构。
我写了一个简短的递归函数可以做到这一点：
ptunwrap <- function(x) if (typeof(x)=='language') lapply(as.list(x),ptunwrap) else x;
那么让我们看一个例子：
f1 <- a+b~c/d;
f1;
## a + b ~ c/d
ptunwrap(f1);
## [[1]]
## `~`
##
## [[2]]
## [[2]][[1]]
## `+`
##
## [[2]][[2]]
## a
##
## [[2]][[3]]
## b
##
##
## [[3]]
## [[3]][[1]]
## `/`
##
## [[3]][[2]]
## c
##
## [[3]][[3]]
## d
##
当分析树包含函数调用时，函数调用表示为递归列表结构中的单个列表节点。函数的符号作为该列表的第一个列表组件嵌入，其参数作为列表的后续列表组件嵌入。
上面的解析树中有3个函数调用。如前所述，顶级列表表示对 `~`() 函数的调用。第二个顶级列表组件进一步分支到另一个列表，该列表包含对 `+`() 函数的调用，该函数本身有两个参数，符号 a 和 b。第三个组件类似，表示对 `/`() 函数的调用并再次采用两个参数，符号 c 和 d.
重要的是要理解，虽然解析树始终表示语法上有效的 R 代码，并且具有能力 被评估以产生单个结果值，但它不是 需要 来评估解析树。完全有可能创建一个解析树而不对其求值。
那么，创建一个您从不求值的解析树的目的是什么？在 R 中，通常这样做是为了使用方便的语法将某些信息传递给函数。
作为一个随机示例，data.table 包允许使用 := 运算符使用以下语法糖将列添加到现有 data.table：
library(data.table);
dt <- data.table(a=1:3);
dt[,b:=a*2L];
dt;
## a b
## 1: 1 2
## 2: 2 4
## 3: 3 6
在内部，data.table 使用非标准参数评估来检索 `[.data.table`() 函数的第二个参数（技术上在函数定义中的第三个；语法糖形式中的第二个）的解析树，通常被称为“j 参数”，因为参数名称是 j。如果需要，您实际上可以直接在 R 中检查源代码本身。这是最相关的代码片段：
data.table:::`[.data.table`;
## function (x, i, j, by, keyby, with = TRUE, nomatch = getOption("datatable.nomatch"),
## mult = "all", roll = FALSE, rollends = if (roll == "nearest") c(TRUE,
## TRUE) else if (roll >= 0) c(FALSE, TRUE) else c(TRUE,
## FALSE), which = FALSE, .SDcols, verbose = getOption("datatable.verbose"),
## allow.cartesian = getOption("datatable.allow.cartesian"),
## drop = NULL, on = NULL)
## {
##
## ... snip ...
##
## if (!missing(j)) {
## jsub = substitute(j)
## if (is.call(jsub))
## jsub = construct(deconstruct_and_eval(replace_dot(jsub),
## parent.frame(), parent.frame()))
## }
##
## ... snip ...
##
我们可以看到他们正在使用 substitute(j) 来检索 j 参数的解析树。对于上面的演示，这是他们将得到的：
ptunwrap(substitute(b:=a*2L));
## [[1]]
## `:=`
##
## [[2]]
## b
##
## [[3]]
## [[3]][[1]]
## `*`
##
## [[3]][[2]]
## a
##
## [[3]][[3]]
## [1] 2
##
稍后在代码中，他们测试顶级函数符号是否为 :=，这是支持向 data.table。如果是这样，他们将测试 LHS 是否由单个裸词组成，该裸词被视为要添加（或修改或删除）的列的名称。请注意，在这种情况下，实际上 不可能 他们实际评估解析树的 LHS，因为它包含 [=] 中尚不存在的符号141=]。但是，RHS 最终会被评估以生成要添加到新名称下的 data.table 的列向量。
所以应该清楚的是，公式可以在 R 中的各种上下文中使用，并且它们并不总是被评估。有时会简单地检查解析树以检索从调用者传递给被调用者的信息。即使在 被 评估的上下文中，有时也只有 LHS 或 RHS（或两者）会被独立评估，忽略当时嵌入在解析树中的顶级函数符号它已创建。
转到boxplot()函数，让我们看一下formula参数的文档：
a formula, such as y ~ grp, where y is a numeric vector of data values to be split into groups according to the grouping variable grp (usually a factor).
在这种情况下，公式的两边最终被独立评估，LHS 提供数据向量，RHS 提供分组定义。
证明这一点的一个好方法如下：
boxplot(1:9~1:9%%3L);
注意公式的两边是如何由文字表达式组成的：
1:9;
## [1] 1 2 3 4 5 6 7 8 9
1:9%%3L;
## [1] 1 2 0 1 2 0 1 2 0
在内部，boxplot() 必须独立评估公式的每一边以生成数据和分组向量，几乎就像您将两个表达式作为单独的参数传递一样。
那么，让我们创建一个月度时间序列箱线图的简单演示：
N <- 36L; df <- data.frame(date=seq(as.Date('2016-01-01'),by='month',len=N),y=rnorm(N));
df;
## date y
## 1 2016-01-01 -1.56004488
## 2 2016-02-01 0.65699747
## 3 2016-03-01 0.05729631
## 4 2016-04-01 -0.02092276
## 5 2016-05-01 0.46673530
## 6 2016-06-01 -0.18652580
## 7 2016-07-01 0.06228650
## 8 2016-08-01 1.54452267
## 9 2016-09-01 1.06643594
## 10 2016-10-01 -1.51178160
## 11 2016-11-01 0.82904673
## 12 2016-12-01 0.37667201
## 13 2017-01-01 -0.10135801
## 14 2017-02-01 0.94692462
## 15 2017-03-01 -1.60781946
## 16 2017-04-01 0.47189753
## 17 2017-05-01 -1.32869317
## 18 2017-06-01 -0.49821455
## 19 2017-07-01 0.54474606
## 20 2017-08-01 0.47565264
## 21 2017-09-01 -0.97494730
## 22 2017-10-01 -1.22781588
## 23 2017-11-01 -0.34919086
## 24 2017-12-01 -0.78153843
## 25 2018-01-01 -0.59355220
## 26 2018-02-01 -2.58287605
## 27 2018-03-01 1.42148186
## 28 2018-04-01 -1.01278176
## 29 2018-05-01 -0.80961662
## 30 2018-06-01 0.19793126
## 31 2018-07-01 -1.03072915
## 32 2018-08-01 -0.87896416
## 33 2018-09-01 -2.36216655
## 34 2018-10-01 1.82708221
## 35 2018-11-01 0.05579195
## 36 2018-12-01 1.28612246
boxplot(y~months(date),df);
需要的话可以研究下源码，里面需要追溯S3的查找过程：
boxplot;
## function (x, ...)
## UseMethod("boxplot")
## <bytecode: 0x600b50760>
## <environment: namespace:graphics>
graphics:::boxplot.formula;
## function (formula, data = NULL, ..., subset, na.action = NULL)
## {
## if (missing(formula) || (length(formula) != 3L))
## stop("'formula' missing or incorrect")
## m <- match.call(expand.dots = FALSE)
## if (is.matrix(eval(m$data, parent.frame())))
## m$data <- as.data.frame(data)
## m$... <- NULL
## m$na.action <- na.action
## m[[1L]] <- quote(stats::model.frame)
## mf <- eval(m, parent.frame())
## response <- attr(attr(mf, "terms"), "response")
## boxplot(split(mf[[response]], mf[-response]), ...)
## }
## <bytecode: 0x6035c67f8>
## <environment: namespace:graphics>
它几乎是令人恼火的迂回和复杂，但 graphics:::boxplot.formula() 有效地检索（通过 match.call()) the parse tree that caused its own invocation, massages it a little, most notably replacing its own function symbol boxplot.formula with stats::model.frame, and then evaluates that new parse tree, thereby calling stats::model.frame()。该函数本身非常复杂，并涉及进一步的 S3 查找，但这是最相关的代码：
model.frame;
## function (formula, ...)
## UseMethod("model.frame")
## <bytecode: 0x601464b18>
## <environment: namespace:stats>
model.frame.default;
## function (formula, data = NULL, subset = NULL, na.action = na.fail,
## drop.unused.levels = FALSE, xlev = NULL, ...)
## {
##
## ... snip ...
##
## if (!inherits(formula, "terms"))
## formula <- terms(formula, data = data)
## env <- environment(formula)
## rownames <- .row_names_info(data, 0L)
## vars <- attr(formula, "variables")
## predvars <- attr(formula, "predvars")
## if (is.null(predvars))
## predvars <- vars
## varnames <- sapply(vars, function(x) paste(deparse(x, width.cutoff = 500),
## collapse = " "))[-1L]
## variables <- eval(predvars, data, env)
##
## ... snip ...
##
因此，最终，它从公式对象中检索各个表达式，并使用 eval() 以给定的 data.frame 和公式的闭包环境作为上下文对它们进行评估，从而给出结果向量:
attr(terms(y~months(date),data=df),'variables');
## list(y, months(date))
eval(attr(terms(y~months(date),data=df),'variables'),df);
## [[1]]
## [1] -1.56004488 0.65699747 0.05729631 -0.02092276 0.46673530 -0.18652580 0.06228650 1.54452267 1.06643594 -1.51178160 0.82904673 0.37667201 -0.10135801 0.94692462 -1.60781946 0.47189753 -1.32869317 -0.49821455 0.54474606
## [20] 0.47565264 -0.97494730 -1.22781588 -0.34919086 -0.78153843 -0.59355220 -2.58287605 1.42148186 -1.01278176 -0.80961662 0.19793126 -1.03072915 -0.87896416 -2.36216655 1.82708221 0.05579195 1.28612246
##
## [[2]]
## [1] "January" "February" "March" "April" "May" "June" "July" "August" "September" "October" "November" "December" "January" "February" "March" "April" "May" "June" "July"
## [20] "August" "September" "October" "November" "December" "January" "February" "March" "April" "May" "June" "July" "August" "September" "October" "November" "December"
##
为了重申之前提出的观点，请注意 `~`() 函数在评估过程中是如何无处可寻的。这是 R 公式的任意实现细节，它们将 `~`() 函数编码为公式对象的解析树存储表示中的顶级函数符号。公式边的实际求值（如果发生）不涉及该函数的求值。
最后，让我们考虑一下如果您实际评估包含公式存储表示的整个解析树会发生什么。回想一下，`~`() 函数只是根据其参数创建一个公式。因此，评估公式具有吐出与刚刚评估的完全相同的公式的有趣效果：
f1;
## a + b ~ c/d
eval(f1);
## a + b ~ c/d
eval(eval(f1));
## a + b ~ c/d
我已经写了几个关于解析树和公式的其他答案，您可能会感兴趣：