【发布时间】：2018-10-03 13:39:45
【问题描述】：
假设我有一个依赖 y(t) 和参数 p1、p2 和 p3 的数据
这可能会影响值y(t)。
我创建了 3 个线性方程，它们取决于以下组合
参数p1 和p2 - p3 对y(t) 没有影响，这意味着它遵循随机分配。
您可以在问题的末尾找到一个可重现的示例。
这三个方程是
p1 p2 Equation
1 1 5 + 3t
2 1 1 - t
2 2 3 + t
包含随机数据的 3 个方程的图如下所示：
现在，如果我根据我的随机数据调用lm()（有关公式，请参阅here），我会得到以下结果。
lm(formula = y ~ .^2, data = mydata)
Residuals:
Min 1Q Median 3Q Max
-1.14707 -0.22785 0.00157 0.23099 1.10528
Coefficients: (6 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.83711 0.17548 27.565 <2e-16 ***
t 2.97316 0.02909 102.220 <2e-16 ***
p12 -3.86697 0.21487 -17.997 <2e-16 ***
p22 2.30617 0.20508 11.245 <2e-16 ***
p23 NA NA NA NA
p32 0.16518 0.21213 0.779 0.4375
p33 0.23450 0.22594 1.038 0.3012
t:p12 -4.00574 0.03119 -128.435 <2e-16 ***
t:p22 2.01230 0.03147 63.947 <2e-16 ***
t:p23 NA NA NA NA
t:p32 0.01155 0.03020 0.383 0.7027
t:p33 0.02469 0.03265 0.756 0.4508
p12:p22 NA NA NA NA
p12:p23 NA NA NA NA
p12:p32 -0.10368 0.21629 -0.479 0.6325
p12:p33 -0.11728 0.21386 -0.548 0.5843
p22:p32 -0.20871 0.19633 -1.063 0.2896
p23:p32 NA NA NA NA
p22:p33 -0.44250 0.22322 -1.982 0.0495 *
p23:p33 NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4112 on 136 degrees of freedom
Multiple R-squared: 0.9988, Adjusted R-squared: 0.9987
F-statistic: 8589 on 13 and 136 DF, p-value: < 2.2e-16
如果我只想考虑具有重要意义的参数，我会主张忽略接近零的参数。如果我理解正确，零参数不会导致“新行”。然后我得到以下简化模型（为了便于阅读，对值进行了四舍五入）：
Estimate
(Intercept) 5 ***
t 3 ***
p12 -4 ***
p22 2 ***
t:p12 -4 ***
t:p22 2 ***
然后我将根据估计重建理论模型如下 上面（只有非常重要的参数！）：
p1 p2 Equation Result
1 1 5+3t 5+3t
1 2 5+3t+p22+t:p22*t 7+5t
2 1 5+3t+p12+t:p12*t 1-t
2 2 5+3t+p22+t:p22*t+p12+t:p12*t 3+t
现在，7 + 5t 显然是错误的，但我不确定原因。
我猜，lm 依次添加参数，从而对应的模型
y ~ t:p2不包含在上面的模型中？
这个question 和其中的引用可能是相关的，但我没有查看lm 的结果 - 所以没有什么。
可重现的例子：
r <- generate_3lines(sigma = 0.5, slopes = c(3, 1, -1), offsets = c(5, 3, 1))
t_m <- r$t_m; y_m <- r$y_m; y_t <- r$y_t; rm(r)
mydata <- generate_randomdata(t_m, y_m, y_t)
# What the raw data looks like:
plot(t_m[[1]], y_t[[1]], type = "l", lty = 3, col = "black", main = "Raw data",
xlim = c(0, 10), ylim = c(min(mydata$y), max(mydata$y)), xlab = "t", ylab = "y")
lines(t_m[[2]], y_t[[2]], col = "black", lty = 3)
lines(t_m[[3]], y_t[[3]], col = "black", lty = 3)
points(x = mydata$t, y = mydata$y)
fit <- lm(y ~ .^2, data = mydata) # Not all levels / variables are linearly
print(summary(fit))
和功能
generate_3lines <- function(sigma = 0.5, slopes = c(3, 1, -1), offsets = c(5, 3, 1)) {
t <- seq(0,10, length.out = 1000) # large sample of x values
t_m <- list()
y_m <- list()
y_t <- list()
for (i in 1:3) {
set.seed(33*i)
t_m[[i]] <- sort(sample(t, 50, replace = F))
set.seed(33*i)
noise <- rnorm(10, 0, sigma)
y_m[[i]] <- slopes[i]*t_m[[i]] + offsets[i] + noise
y_t[[i]] <- slopes[i]*t_m[[i]] + offsets[i]
}
return(list(t_m = t_m, y_m = y_m, y_t = y_t))
}
generate_randomdata <- function(t_m, y_m, y_t) {
# Final data set
df1 <- data.frame(t = t_m[[1]], y = y_m[[1]], p1 = rep(1), p2 = rep(1),
p3 = sample(c(1, 2, 3), length(t_m[[1]]), replace = T))
df2 <- data.frame(t = t_m[[2]], y = y_m[[2]], p1 = rep(2), p2 = rep(2),
p3 = sample(c(1, 2, 3), length(t_m[[1]]), replace = T))
df3 <- data.frame(t = t_m[[3]], y = y_m[[3]], p1 = rep(2), p2 = rep(3),
p3 = sample(c(1, 2, 3), length(t_m[[1]]), replace = T))
mydata <- rbind(df1, df2, df3)
mydata$p1 <- factor(mydata$p1)
mydata$p2 <- factor(mydata$p2)
mydata$p3 <- factor(mydata$p3)
mydata <- mydata[sample(nrow(mydata)), ]
return(mydata)
}
@MrFlick 输入后编辑：问题现在也在Cross Validated
评论：看来，ggplot 中的拟合并不是真正自动化的，请参阅 here
【问题讨论】：
-
这真的是关于编程的问题吗？似乎您只需要帮助解释统计模型。 Cross Validated 更适合统计问题。
-
@MrFlick 我无法想象没有 r 知识的人知道如何解释 r 的输出。
-
尽管有相反的谣言，但大多数统计学家并不使用纸和笔来完成所有工作。许多人使用计算机，甚至许多人使用 R。大多数人应该能够解释
lm()调用的结果。但实际上，您不能只从模型中挑选重要系数。 -
@MrFlick 我会争辩说我没有选择参数，而是我将参数（接近 0）设置为 0，因此它们不会导致不同的 y 值。
-
@MrFlick：感谢您将我指向 Cross Validated。我发布了问题there。