在R中执行线性回归时使用实验不确定性

Question

我有一个实验数据集，希望将其拟合为多项式。该数据包含一个自变量、一个因变量以及后者测量中的不确定性，例如：

2000  0.2084272   0.002067834
2500  0.207078125 0.001037248
3000  0.2054202   0.001959138
3500  0.203488075 0.000328942
4000  0.2013152   0.000646088
4500  0.198933825 0.001375657
5000  0.196375    0.000908696
5500  0.193668575 0.00014721
6000  0.1908432   0.000526976
6500  0.187926325 0.001217318
7000  0.1849442   0.000556495
7500  0.181921875 0.000401883
8000  0.1788832   0.001446992
8500  0.175850825 0.001235017
9000  0.1728462   0.001676249
9500  0.169889575 0.001011735
10000 0.167       0.000326678

(列 x、y、±y )。

使用上述内容，我可以进行多项式拟合，例如：

mydata = read.table("example.txt")
model <- lm(V2~V1+I(V1^2)+I(V1^3)+I(V1^4), data = mydata)

但这并没有利用不确定性的值。我该如何告诉R，数据集的第三列是一个不确定因素，因此应在回归分析中使用？

- Joseph Wright

---

请注意，这里的数据是“模拟”的，基于但不使用真实数据。 - Joseph Wright

你希望如何使用不确定性？作为另一个独立变量？还是其他？请自己开心一点，不要养成附加数据的习惯。这会使您的环境杂乱无章，并可能导致组合数据和排序等问题。 - Heroka

@Heroka 在图形分析程序（Origin，Igor，...）中，将一列用作不确定性在分析中是相当标准的：我不是统计学家，所以不知道它在此之外如何使用。关于“附加”，我猜你指的是attach(data)：我从《R语言实战》（第二版，例如第467页）中得到了这个，所以认为它是标准的。 - Joseph Wright

对不起，我对不确定性的事情了解有限。关于附加操作：这不是标准做法，但在初学者课程中经常教授。以我个人的观点来看，这不是一个好习惯，因为随着时间的推移，你的分析会变得更加复杂。（分组操作、排序、导出数据和多个数据集都更容易通过将所有内容包含在一个对象中来实现。这样可以避免出现难以调试的错误。这里有更多关于这个主题的讨论。） - Heroka

Answer 1

当因变量的测量误差与自变量不相关时，估计的系数是无偏的，但标准误差太小。这是我使用的参考文献（第1页和第2页）：http://people.stfx.ca/tleo/econ370term2lec4.pdf 我认为你只需要调整lm（）计算出来的标准误差。这就是我在下面的代码中尝试做的事情。我不是统计学专业人士，所以你可能想发帖到交叉验证并寻求更好的直觉。

对于下面的示例，我假设“不确定性”列是标准差（或标准误差）。为简单起见，我将模型更改为：y〜x。

# initialize dataset
df <- data.frame(
        x = c(2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000,7500,8000,8500,9000,9500,10000),
        y = c(0.2084272,0.207078125,0.2054202,0.203488075,0.2013152,0.198933825, 0.196375,0.193668575, 0.1908432, 0.187926325,   0.1849442,   0.181921875, 0.1788832, 0.175850825, 0.1728462,0.169889575,  0.167),
        y.err = c(0.002067834, 0.001037248,  0.001959138, 0.000328942, 0.000646088, 0.001375657, 0.000908696, 0.00014721, 0.000526976, 0.001217318, 0.000556495, 0.000401883, 0.001446992, 0.001235017, 0.001676249, 0.001011735, 0.000326678)
)

df

#  model regression
model <- lm(y~x, data = df)
summary(model)

#  get the variance of the measurement error
#  thanks to: http://schools-wikipedia.org/wp/v/Variance.htm
#  law of total variance
pooled.var.y.err <- mean((df$y.err)^2) + var((df$y.err)^2)

# get variance of beta from model
#   thanks to: http://stats.stackexchange.com/questions/44838/how-are-the-standard-errors-of-coefficients-calculated-in-a-regression
X <- cbind(1, df$x)
#     (if you add more variables to the model you need to modify the following line)
var.betaHat <- anova(model)[[3]][2] * solve(t(X) %*% X) 

# add betaHat variance to measurement error variance
var.betaHat.adj <- var.betaHat + pooled.var.y.err

# calculate adjusted standard errors 
sqrt(diag(var.betaHat.adj))

# compare to un-adjusted standard errors
sqrt(diag(var.betaHat))

- Andrew Leidner