r - R中的线性回归(正态和对数数据)

Question

我想在 R 中对正态和双对数图中的数据进行线性回归。

对于 正常数据数据集可能如下:

lin <- data.frame(x = c(0:6), y = c(0.3, 0.1, 0.9, 3.1, 5, 4.9, 6.2))
plot (lin$x, lin$y)

在那里我想计算只为数据点 2、3 和 4 的线性回归画一条线。

对于 双对数数据数据集可能如下:

data = data.frame(
    x=c(1:15),
    y=c(
        1.000, 0.742, 0.623, 0.550, 0.500, 0.462, 0.433,
        0.051, 0.043, 0.037, 0.032, 0.028, 0.025, 0.022, 0.020
      )
    )
plot (data$x, data$y, log="xy")

在这里，我想为数据集 1:7 和 8:15 绘制回归线。

我可以计算斜率和 y 偏移量以及拟合参数(R^2，p 值)吗？

对于正态数据和对数数据，它是如何完成的？

谢谢你的帮助，

斯文

Answer 1

在 R 中，线性最小二乘模型通过 lm() 拟合功能。使用公式界面我们可以使用subset选择用于拟合实际模型的数据点的参数，例如:

lin <- data.frame(x = c(0:6), y = c(0.3, 0.1, 0.9, 3.1, 5, 4.9, 6.2))
linm <- lm(y ~ x, data = lin, subset = 2:4)

给予:

R> linm

Call:
lm(formula = y ~ x, data = lin, subset = 2:4)

Coefficients:
(Intercept)            x  
     -1.633        1.500  

R> fitted(linm)
         2          3          4 
-0.1333333  1.3666667  2.8666667

至于双日志，我猜你有两种选择； i) 像我们上面所做的那样估计两个单独的模型，或者 ii) 通过 ANCOVA 估计。使用 log() 在公式中完成对数转换.

通过两个独立的模型:

logm1 <- lm(log(y) ~ log(x), data = dat, subset = 1:7)
logm2 <- lm(log(y) ~ log(x), data = dat, subset = 8:15)

或者通过 ANCOVA，我们需要一个指标变量

dat <- transform(dat, ind = factor(1:15 <= 7))
logm3 <- lm(log(y) ~ log(x) * ind, data = dat)

您可能会问这两种方法是否等效？好吧，我们可以通过模型系数来证明这一点。

R> coef(logm1)
  (Intercept)        log(x) 
-0.0001487042 -0.4305802355 
R> coef(logm2)
(Intercept)      log(x) 
  0.1428293  -1.4966954

所以对于单独的模型，两个斜率是 -0.4306 和 -1.4967。 ANCOVA 模型的系数为:

R> coef(logm3)
   (Intercept)         log(x)        indTRUE log(x):indTRUE 
     0.1428293     -1.4966954     -0.1429780      1.0661152

我们如何调和两者？好吧，我设置的方式 ind , logm3参数化以提供从 logm2 估计的更直接的值; logm2 的截距和 logm3与 log(x) 的系数相同.获得等价于系数的值
的 logm1 ，我们需要做一个操作，首先进行拦截:

R> coefs[1] + coefs[3]
  (Intercept) 
-0.0001487042

其中 indTRUE 的系数是第 1 组的平均值与第 2 组的平均值的差异。对于斜率:

R> coefs[2] + coefs[4]
    log(x) 
-0.4305802

这与我们得到的 logm1 相同并且基于组 2 的斜率 ( coefs[2] ) 修改为组 1 的斜率差异 ( coefs[4] )。

至于绘图，一个简单的方法是通过 abline()对于简单模型。例如。对于正常数据示例:

plot(y ~ x, data = lin)
abline(linm)

对于日志数据，我们可能需要更有创意，这里的一般解决方案是预测数据范围并绘制预测:

pdat <- with(dat, data.frame(x = seq(from = head(x, 1), to = tail(x,1), 
                                     by = 0.1))
pdat <- transform(pdat, yhat = c(predict(logm1, pdat[1:70,, drop = FALSE]), 
                                 predict(logm2, pdat[71:141,, drop = FALSE])))

可以通过对 yhat 取幂来绘制原始比例

plot(y ~ x, data = dat)
lines(exp(yhat) ~ x, dat = pdat, subset = 1:70, col = "red")
lines(exp(yhat) ~ x, dat = pdat, subset = 71:141, col = "blue")

或在对数刻度上:

plot(log(y) ~ log(x), data = dat)
lines(yhat ~ log(x), dat = pdat, subset = 1:70, col = "red")
lines(yhat ~ log(x), dat = pdat, subset = 71:141, col = "blue")

例如...

这种通用解决方案也适用于更复杂的 ANCOVA 模型。在这里我像以前一样创建一个新的 pdat 并添加一个指标

pdat <- with(dat, data.frame(x = seq(from = head(x, 1), to = tail(x,1), 
                                     by = 0.1)[1:140],
                             ind = factor(rep(c(TRUE, FALSE), each = 70))))
pdat <- transform(pdat, yhat = predict(logm3, pdat))

请注意我们如何从对 predict() 的单个调用中获得我们想要的所有预测。因为使用了ANCOVA来贴合 logm3 .我们现在可以像以前一样绘制:

plot(y ~ x, data = dat)
lines(exp(yhat) ~ x, dat = pdat, subset = 1:70, col = "red")
lines(exp(yhat) ~ x, dat = pdat, subset = 71:141, col = "blue")