在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)
或通过协方差分析,我们需要变
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")
我该如何获得一个名为num的单个数字? 'str(coef(daten_fit))'给出: '命名的数字0.8 - attr(*,“names”)= chr“x”',但它是否可能问coef以名称“x ”。我尝试了不同的想法,比如'coeff(daten_fit)$ x'或'coeff(daten_fit)[1]'或'attr(coeff(daten_fit),“x”)',但是没有任何工作。是不可能通过名字得到价值? – 2011-06-08 17:59:08
@Sven它是一个命名的数字向量。请注意,它是'coef()'与一个“f”而不是'coeff()'与两个。 '我的答案中的'logm2'对象为''coef(logm2)[“(Intercept)”]和'coef(logm2)[“log(x)”]你不能在原子向量(一种基本类型的向量)上使用'$',它只能用于列表(以及当然是列表的数据帧)。 – 2011-06-08 18:10:25