數據具有不確定性的線性模型，使用 R
假設我有一些不確定的數據。例如：
X Y 1 10±4 2 50±3 3 80±7 4 105±1 5 120±9例如，不確定性的性質可以是重複測量或實驗，或測量儀器的不確定性。
我想用 R 擬合一條曲線，通常我會用
lm. 但是，當它給我擬合係數的不確定性以及預測區間的不確定性時，這並沒有考慮到數據中的不確定性。查看文檔，該lm頁面具有以下內容：…權重可用於表示不同的觀測值具有不同的方差…
所以這讓我覺得這可能與它有關。我知道手動執行此操作的理論，但我想知道是否可以使用該
lm功能執行此操作。如果沒有，是否有任何其他功能（或包）能夠做到這一點？編輯
看到一些評論，這裡是一些澄清。舉個例子：
x <- 1:10 y <- c(131.4,227.1,245,331.2,386.9,464.9,476.3,512.2,510.8,532.9) mod <- lm(y ~ x + I(x^2)) summary(mod)給我：
Residuals: Min 1Q Median 3Q Max -32.536 -8.022 0.087 7.666 26.358 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 39.8050 22.3210 1.783 0.11773 x 92.0311 9.3222 9.872 2.33e-05 *** I(x^2) -4.2625 0.8259 -5.161 0.00131 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 18.98 on 7 degrees of freedom Multiple R-squared: 0.986, Adjusted R-squared: 0.982 F-statistic: 246.7 on 2 and 7 DF, p-value: 3.237e-07所以基本上，我的係數是 a=39.8±22.3，b=92.0±9.3，c=-4.3±0.8。現在讓我們說對於每個數據點，錯誤是 20。我將
weights = rep(20,10)在lm調用中使用，我得到了這個：Residual standard error: 84.87 on 7 degrees of freedom但係數上的標準誤差不會改變。
手動，我知道如何使用矩陣代數計算協方差矩陣並將權重/誤差放在那裡，並使用它推導置信區間。那麼有沒有辦法在 lm 函數本身或任何其他函數中做到這一點？
這種類型的模型實際上在某些科學分支（例如物理學）和工程學中比“正常”線性回歸更常見。因此，在諸如 之類的物理工具中
ROOT，進行這種類型的擬合是微不足道的，而線性回歸併不是本機實現的！物理學家傾向於將其稱為“擬合”或卡方最小化擬合。正態線性回歸模型假設存在總體方差附加到每個測量。然後最大化可能性
或等價於它的對數
因此，最小二乘的名稱——最大化似然性與最小化平方和相同，並且是一個不重要的常數，只要它是常數。對於具有不同已知不確定性的測量，您需要最大化
或等價於它的對數
因此，您實際上想通過逆方差對測量值進行加權，而不是方差。這是有道理的——更準確的測量具有更小的不確定性，應該給予更多的權重。請注意，如果此權重是恆定的，它仍然會從總和中扣除。因此，它不會影響估計值，但應該會影響標準誤差，取自. 但是，在這裡我們遇到了物理學/科學與統計學之間的另一個區別。通常在統計中，您期望兩個變量之間可能存在相關性，但很少會是準確的。另一方面，在物理學和其他科學中，您通常期望相關性或關係是準確的，如果不是因為討厭的測量錯誤（例如， 不是）。您的問題似乎更多地屬於物理/工程案例。因此，
lm對與您的測量結果和權重相關的不確定性的解釋與您想要的並不完全相同。它會承受重量，但它仍然認為有一個整體考慮回歸誤差，這不是你想要的——你希望你的測量誤差是唯一的誤差。（ 的lm解釋的最終結果是，只有權重的相對值很重要，這就是為什麼您作為測試添加的恆定權重沒有效果的原因）。這裡的問答有更多細節：那裡的答案中有幾個可能的解決方案。特別是，那裡的匿名答案建議使用
vcov(mod)/summary(mod)$sigma^2基本上，
lm根據估計的協方差矩陣縮放，並且您想要撤消此操作。然後，您可以從校正後的協方差矩陣中獲得所需的信息。試試這個，但如果可以使用手動線性代數，請嘗試仔細檢查。請記住，權重應該是逆方差。編輯
如果您經常做這種事情，您可能會考慮使用
ROOT（這似乎在本地執行此操作lm，glm但不這樣做）。下面是一個簡短的示例，說明如何在ROOT. 首先，ROOT可以通過 C++ 或 Python 使用，而且下載和安裝量很大。您可以使用 Jupiter notebook 在瀏覽器中嘗試，點擊此處的鏈接，在右側選擇“Binder”，在左側選擇“Python”。import ROOT from array import array import math x = range(1,11) xerrs = [0]*10 y = [131.4,227.1,245,331.2,386.9,464.9,476.3,512.2,510.8,532.9] yerrs = [math.sqrt(i) for i in y] graph = ROOT.TGraphErrors(len(x),array('d',x),array('d',y),array('d',xerrs),array('d',yerrs)) graph.Fit("pol2","S") c = ROOT.TCanvas("test","test",800,600) graph.Draw("AP") c.Draw()我已經把平方根作為不確定性價值觀。擬合的輸出是
Welcome to JupyROOT 6.07/03 **************************************** Minimizer is Linear Chi2 = 8.2817 NDf = 7 p0 = 46.6629 +/- 16.0838 p1 = 88.194 +/- 8.09565 p2 = -3.91398 +/- 0.78028並產生了一個不錯的情節：
ROOT fitter 也可以處理不確定性值，這可能需要對
lm. 如果有人知道在 R 中執行此操作的本地方法，我將有興趣學習它。第二次編輯
@Wolfgang 上一個問題的另一個答案給出了一個更好的解決方案：包中的
rma工具metafor（我最初將該答案中的文本解釋為它沒有計算截距，但事實並非如此）。將測量值 y 中的方差簡化為 y：> rma(y~x+I(x^2),y,method="FE") Fixed-Effects with Moderators Model (k = 10) Test for Residual Heterogeneity: QE(df = 7) = 8.2817, p-val = 0.3084 Test of Moderators (coefficient(s) 2,3): QM(df = 2) = 659.4641, p-val < .0001 Model Results: estimate se zval pval ci.lb ci.ub intrcpt 46.6629 16.0838 2.9012 0.0037 15.1393 78.1866 ** x 88.1940 8.0956 10.8940 <.0001 72.3268 104.0612 *** I(x^2) -3.9140 0.7803 -5.0161 <.0001 -5.4433 -2.3847 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1對於我發現的這種類型的回歸，這絕對是最好的純 R 工具。
