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被parametrised给予更直接地从值估计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")