r - R 的 lm 函数如何处理因子级别(在 C_Cdqrls 中？)？

Question

或者换句话说:在这种情况下使用哪种算法？我猜他们使用判别分析作为描述，例如在第 4.4 章。在詹姆斯等。阿尔。 “使用 R 中的应用程序进行统计学习简介”？

从评论中输入后，我还可以将问题重述如下:

魔法的第一部分出现在ans <- .External2(C_modelmatrix, t, data) (在 model.matrix.default 中)其中模型根据因子水平变化 => 我想我理解这部分。

第二部分还涉及z <- .Call(C_Cdqrls, x, y, tol, FALSE)我没想到，线性回归和判别分析在数学层面上是一样的。
我是否错过了一些明显的东西 ?再次，我的 stats包是二进制文件，我无权访问源代码...

我在此 article 中找到了非常有用的解释，但在某些时候它只说明

... This [factor] deconstruction can be a complex task, so we will not go into details lest it take us too far afield...

我在文档中找不到任何内容，也无法理解使用 debug(lm) 的情况。
我使用可重现的示例理解的内容:

n <- 10
p <- 6
set.seed(1)
x <- seq(0, 20, length.out = n) + rnorm(n, 0, 1)
y <- c(1:3)
y <- sample(y, n, replace = TRUE)
z <- 10*y*x + 10*y + 10 + rnorm(n, 0, 1)
debug(lm)
fit <- lm(z ~ x*y)

后 mt <- attr(mf, "terms")看起来像

mt
# ...
# attr(,"dataClasses")
#         z         x         y 
# "numeric" "numeric" "numeric" 

而之后

n <- 10
p <- 6
set.seed(1)
x <- seq(0, 20, length.out = n) + rnorm(n, 0, 1)
y <- c(1:3)
y <- sample(y, n, replace = TRUE)
z <- 10*y*x + 10*y + 10 + rnorm(n, 0, 1)
y <- as.factor(y)
debug(lm)
fit <- lm(z ~ x*y)

和 mt <- attr(mf, "terms")好像

mt
# ...
# attr(,"dataClasses")
#         z         x         y 
# "numeric" "numeric"  "factor"

但后来似乎，他们总是打电话 lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...)还有 z <- .Call(C_Cdqrls, x, y, tol, FALSE)我认为只有在没有因素的情况下才有效。
上面的链接很好地解释了模型矩阵和 qr 分解的所有内容，我认为这在因素的情况下不起作用。

编辑 : x <- model.matrix(mt, mf, contrasts)之后的模型矩阵已经不同了。如果是数字

x
   (Intercept)          x y       x:y
1            1 -0.6264538 3 -1.879361
2            1  2.4058655 1  2.405866
3            1  3.6088158 2  7.217632
4            1  8.2619475 1  8.261947
5            1  9.2183967 1  9.218397
6            1 10.2906427 2 20.581285
7            1 13.8207624 1 13.820762
8            1 16.2938803 2 32.587761
9            1 18.3535591 3 55.060677
10           1 19.6946116 2 39.389223
attr(,"assign")
[1] 0 1 2 3

在因素的情况下

x
   (Intercept)          x y2 y3      x:y2       x:y3
1            1 -0.6264538  0  1  0.000000 -0.6264538
2            1  2.4058655  0  0  0.000000  0.0000000
3            1  3.6088158  1  0  3.608816  0.0000000
4            1  8.2619475  0  0  0.000000  0.0000000
5            1  9.2183967  0  0  0.000000  0.0000000
6            1 10.2906427  1  0 10.290643  0.0000000
7            1 13.8207624  0  0  0.000000  0.0000000
8            1 16.2938803  1  0 16.293880  0.0000000
9            1 18.3535591  0  1  0.000000 18.3535591
10           1 19.6946116  1  0 19.694612  0.0000000
attr(,"assign")
[1] 0 1 2 2 3 3
attr(,"contrasts")
attr(,"contrasts")$`y`
[1] "contr.treatment"

编辑 2 : 部分问题也可以找到 here

Answer 1

借助 this question 的答案，我意识到，答案很简单:

如果因子属于变量(预测变量)，则model.matrix只会变大。因此很明显，C_Cdqrls可以处理模型矩阵。

仅当因变量包含因子、线性回归或 lm不能正常工作，判别分析是一种可能性。 (乍一看，stats::glm 似乎使用了 logit 模型。

来自 Wikipedia :

Discriminant function analysis is very similar to logistic regression, and both can be used to answer the same research questions. Logistic regression does not have as many assumptions and restrictions as discriminant analysis. However, when discriminant analysis’ assumptions are met, it is more powerful than logistic regression. Unlike logistic regression, discriminant analysis can be used with small sample sizes. It has been shown that when sample sizes are equal, and homogeneity of variance/covariance holds, discriminant analysis is more accurate. With all this being considered, logistic regression has become the common choice, since the assumptions of discriminant analysis are rarely met.

例子:

x <- seq(0, 10, length.out = 21)
y <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
y <- as.factor(y)
df <- data.frame(x = x, y = y)

# see ??numeric and the ‘Warning’ section in factor:
plot(x, as.numeric(levels(y))[y], ylim = c(0, 1.2))

fit <- lm(y ~ x, data = df)
print(summary(fit))

fit_glm <- stats::glm(y ~ x, family = binomial(link = "logit"), data = df, control = list(maxit = 50))
print(summary(fit_glm))

df$glm.probs <- stats::predict(fit_glm, newdata = df, type = "response")
df$glm.pred = ifelse(glm.probs > 0.5, 1, 0)
points(x, df$glm.pred + 0.05, col = "red")