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

【问题标题】：Using experimental uncertainty when performing linear regression in R在 R 中执行线性回归时使用实验不确定性
【发布时间】：2016-03-28 21:07:22
【问题描述】：

我有一个想要拟合多项式的实验数据集。数据包括自变量、因变量和后者测量中的不确定性，例如

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 数据集的第三列是不确定的，因此应该在回归分析中使用它？

【问题讨论】：

请注意，这里的数据是“模拟的”，基于但不是使用真实的东西。
您希望如何使用不确定性？作为另一个自变量？作为别的东西？请帮自己一个忙，不要养成附加数据的习惯。它会使您的环境变得混乱，并可能导致分组数据和排序等问题。
@Heroka 在图形分析程序（Origin，Igor，...）中使用一个列作为不确定性是非常标准的：我不是统计学家，所以我不知道它是怎么回事除此之外使用。关于“附加”，我猜你的意思是 attach(data)：我从 The R Book（第二版，eg 第 467 页）中得到了这个，所以假设（d ) 这是标准的。
不确定性超出了我的知识范围，抱歉。关于附加：它不是标准的，但经常教（尤其是在初学者课程中）。在我看来，这不是一个好习惯，因为随着时间的推移，你的分析会变得更加复杂。（分组操作、排序、导出数据和多个数据集都更容易，所有内容都包含在一个对象中。它可以防止你犯难以调试的错误。here 是关于这个主题的更多讨论。

标签： r regression lm

【解决方案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))

【讨论】：