Линейная модель, где данные имеют неопределенность, используя 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

В общем, мои коэффициенты а = 39,8 ± 22,3, b = 92,0 ± 9,3, с = -4,3 ± 0,8. Теперь допустим, что для каждой точки данных ошибка равна 20. Я буду использовать weights = rep(20,10)в lmвызове, и вместо этого получу следующее:

Residual standard error: 84.87 on 7 degrees of freedom

но стандартные ошибки на коэффициенты не меняются.

Вручную, я знаю, как сделать это с помощью вычисления ковариационной матрицы с использованием матричной алгебры и размещения в ней весов / ошибок, а также получения доверительных интервалов, используя это. Так есть ли способ сделать это в самой функции lm или любой другой функции?

Этот тип модели на самом деле гораздо чаще встречается в определенных отраслях науки (например, в физике) и технике, чем "нормальная" линейная регрессия. Таким образом, в таких физических инструментах, как подобный ROOTподход, тривиально, а линейная регрессия изначально не реализована! Физики имеют тенденцию называть это просто «соответствием» или уменьшением соответствия хи-квадрат.

$\sigma$

$F=ma$$F=ma+\epsilon$lm$\sigma^2$lm

лм весов и стандартная ошибка

В ответах есть несколько возможных решений. В частности, анонимный ответ предлагает использовать

vcov(mod)/summary(mod)$sigma^2

lm$\sigma$

РЕДАКТИРОВАТЬ

Если вы много делаете такого рода вещи, вы можете подумать об использовании ROOT(который, кажется, делает это изначально, lmа glmне делает). Вот краткий пример того, как это сделать ROOT. Во-первых, ROOTего можно использовать через C ++ или Python, и его можно скачать и установить. Вы можете попробовать его в браузере, используя блокнот Jupiter, перейдя по ссылке здесь , выбрав «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()

$y$

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    

и получается хороший сюжет:

$x$lm

ВТОРОЕ РЕДАКТИРОВАНИЕ

Другой ответ от того же предыдущего вопроса от @Wolfgang дает еще лучшее решение: rmaинструмент из metaforпакета (я первоначально интерпретировал текст в этом ответе как означающий, что он не вычислял перехват, но это не так). Принимая дисперсии в измерениях у просто:

> 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 для этого типа регрессии, который я нашел.