Model linear dimana data memiliki ketidakpastian, menggunakan R

Katakanlah saya memiliki data yang memiliki ketidakpastian. Sebagai contoh:

X  Y
1  10±4
2  50±3
3  80±7
4  105±1
5  120±9

Sifat ketidakpastian dapat berupa pengukuran ulang atau eksperimen, atau misalnya ketidakpastian instrumen pengukuran.

Saya ingin menyesuaikan kurva menggunakan R, sesuatu yang biasanya saya lakukan lm. Namun, ini tidak memperhitungkan ketidakpastian dalam data saat itu memberi saya ketidakpastian dalam koefisien fit, dan akibatnya interval prediksi. Melihat dokumentasi, lmhalaman memiliki ini:

... bobot dapat digunakan untuk menunjukkan bahwa pengamatan yang berbeda memiliki varian yang berbeda ...

Jadi itu membuat saya berpikir bahwa mungkin ini ada hubungannya dengan itu. Saya tahu teori melakukannya secara manual, tetapi saya bertanya-tanya apakah mungkin melakukan itu dengan lmfungsinya. Jika tidak, apakah ada fungsi lain (atau paket) yang mampu melakukan ini?

EDIT

Melihat beberapa komentar, berikut adalah beberapa klarifikasi. Ambil contoh ini:

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)

Memberi saya:

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

Jadi pada dasarnya, koefisien saya adalah a = 39,8 ± 22,3, b = 92,0 ± 9,3, c = -4,3 ± 0,8. Sekarang katakanlah untuk setiap titik data, kesalahannya adalah 20. Saya akan menggunakan weights = rep(20,10)dalam lmpanggilan dan saya mendapatkan ini sebagai gantinya:

Residual standard error: 84.87 on 7 degrees of freedom

tetapi kesalahan std pada koefisien tidak berubah.

Secara manual, saya tahu bagaimana cara melakukannya dengan menghitung matriks kovarians menggunakan aljabar matriks dan meletakkan bobot / kesalahan di sana, dan menurunkan interval kepercayaan menggunakan itu. Jadi apakah ada cara untuk melakukannya dalam fungsi lm itu sendiri, atau fungsi lainnya?

Jenis model ini sebenarnya jauh lebih umum di cabang ilmu tertentu (misalnya fisika) dan teknik daripada regresi linier "normal". Jadi, dalam alat fisika seperti ROOT, melakukan jenis fit ini sepele, sedangkan regresi linier tidak diterapkan secara asli! Fisikawan cenderung menyebut ini hanya "fit" atau fit meminimalkan chi-square.

$\sigma$

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

lm bobot dan kesalahan standar

Ada beberapa solusi yang mungkin diberikan dalam jawaban di sana. Secara khusus, jawaban anonim di sana menyarankan menggunakan

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

lm$\sigma$

EDIT

Jika Anda sering melakukan hal semacam ini, Anda mungkin mempertimbangkan untuk menggunakan ROOT(yang tampaknya melakukan ini secara sementara lmdan glmtidak). Berikut adalah contoh singkat tentang bagaimana melakukan ini ROOT. Pertama, ROOTdapat digunakan melalui C ++ atau Python, dan ini merupakan unduhan dan instalasi yang sangat besar. Anda dapat mencobanya di browser menggunakan notebook Jupiter, mengikuti tautan di sini , memilih "Binder" di sebelah kanan, dan "Python" di sebelah kiri.

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    

dan plot yang bagus dihasilkan:

$x$lm

EDIT KEDUA

Jawaban lain dari pertanyaan sebelumnya yang sama oleh @ Wolfgang memberikan solusi yang lebih baik: rmaalat dari metaforpaket (saya awalnya menafsirkan teks dalam jawaban itu berarti tidak menghitung intersep, tapi bukan itu masalahnya). Mengambil varians dalam pengukuran y menjadi sekadar 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

Ini jelas merupakan alat R murni terbaik untuk jenis regresi yang saya temukan.