Model liniowy, w którym dane są niepewne, z wykorzystaniem R.

Powiedzmy, że mam dane niepewne. Na przykład:

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

Naturą niepewności może być na przykład powtarzanie pomiarów lub eksperymentów lub niepewność przyrządu pomiarowego.

Chciałbym dopasować do niej krzywą za pomocą R, co normalnie bym zrobił lm. Nie bierze to jednak pod uwagę niepewności danych, gdy daje mi to niepewność co do współczynników dopasowania, aw konsekwencji przedziałów prognozowania. Patrząc na dokumentację, lmstrona ma to:

... wagi mogą być użyte do wskazania, że ​​różne obserwacje mają różne wariancje ...

To sprawia, że ​​myślę, że może to ma coś wspólnego z tym. Znam teorię robienia tego ręcznie, ale zastanawiałem się, czy można to zrobić za pomocą tej lmfunkcji. Jeśli nie, czy jest jakaś inna funkcja (lub pakiet), która jest w stanie to zrobić?

EDYTOWAĆ

Widząc niektóre komentarze, oto wyjaśnienie. Weź ten przykład:

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)

Daje mi:

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

Zasadniczo moje współczynniki wynoszą a = 39,8 ± 22,3, b = 92,0 ± 9,3, c = -4,3 ± 0,8. Powiedzmy teraz, że dla każdego punktu danych błąd wynosi 20. Użyję weights = rep(20,10)w lmwywołaniu, a zamiast tego otrzymuję:

Residual standard error: 84.87 on 7 degrees of freedom

ale błędy standardowe współczynników się nie zmieniają.

Ręcznie wiem, jak to zrobić, obliczając macierz kowariancji za pomocą algebry macierzy i umieszczając w niej wagi / błędy oraz obliczając przedziały ufności za pomocą tego. Czy istnieje sposób, aby to zrobić w samej funkcji lm lub w jakiejkolwiek innej funkcji?

Ten typ modelu jest w rzeczywistości znacznie bardziej powszechny w niektórych gałęziach nauki (np. Fizyka) i inżynierii niż „normalna” regresja liniowa. Zatem w narzędziach fizyki, takich jak ROOTdopasowanie tego typu, jest trywialne, podczas gdy regresja liniowa nie jest natywnie wdrażana! Fizycy nazywają to po prostu „dopasowaniem” lub chi-kwadratem minimalizującym dopasowanie.

Normalne model regresji liniowej zakłada się, że istnieje całkowita wariancja dołączony do każdego pomiaru. Następnie maksymalizuje prawdopodobieństwo lub równoważnie jego logarytm Stąd nazwa najmniejszych kwadratów - maksymalne prawdopodobieństwo to to samo, co minimalizowanie sumy kwadratów, a jest nieistotną stałą, o ile jest stała. Przy pomiarach, które mają różne znane niepewności, będziesz chciał zmaksymalizować $\sigma$

$$L \propto \prod_i e^{-\frac{1}{2} \left( \frac{y_i-(ax_i+b)}{\sigma} \right)^2}$$ $$\log(L) = \mathrm{constant} - \frac{1}{2\sigma^2} \sum_i (y_i-(ax_i+b))^2$$ $\sigma$ $$L \propto \prod e^{-\frac{1}{2} \left( \frac{y-(ax+b)}{\sigma_i} \right)^2}$$ lub równoważnie jego logarytm Tak , w rzeczywistości chcesz zważyć pomiary odwrotną wariancją , a nie wariancją. Ma to sens - dokładniejszy pomiar ma mniejszą niepewność i powinien mieć większą wagę. Zauważ, że jeśli ta waga jest stała, nadal odejmuje się od sumy. Nie wpływa więc na wartości szacunkowe, ale powinno wpływać na standardowe błędy, zaczerpnięte z drugiej pochodnej . $$\log(L) = \mathrm{constant} - \frac{1}{2} \sum \left( \frac{y_i-(ax_i+b)}{\sigma_i} \right)^2$$ $1/\sigma_i^2$$\log(L)$

Tutaj jednak dochodzimy do kolejnej różnicy między fizyką / nauką a całością statystyki. Zazwyczaj w statystykach można się spodziewać korelacji między dwiema zmiennymi, ale rzadko będzie to dokładne. Z drugiej strony w fizyce i innych naukach często oczekuje się, że korelacja lub związek będzie dokładny, choćby nie w przypadku nieznośnych błędów pomiaru (np. , a nie ). Twój problem wydaje się bardziej pasować do przypadku fizyki / inżynierii. W konsekwencji interpretacja niepewności związanej z twoimi pomiarami i wag nie jest dokładnie taka sama, jak tego chcesz. Przyjmie ciężary, ale nadal uważa, że ​​istnieje ogólny$F=ma$$F=ma+\epsilon$lm$\sigma^2$w celu uwzględnienia błędu regresji, który nie jest tym, czego chcesz - chcesz, aby błędy pomiaru były jedynym rodzajem błędu. (Końcowym wynikiem lminterpretacji jest to, że liczą się tylko względne wartości wag, dlatego stałe masy dodane podczas testu nie miały żadnego wpływu). Tutaj pytanie i odpowiedź mają więcej szczegółów:

Wagi i błąd standardowy

Istnieje kilka możliwych rozwiązań podanych w tych odpowiedziach. W szczególności sugeruje tam anonimową odpowiedź

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

Zasadniczo lmskaluje macierz kowariancji w oparciu o jej oszacowany i chcesz to cofnąć. Następnie możesz uzyskać potrzebne informacje z poprawionej macierzy kowariancji. Spróbuj tego, ale spróbuj to dwukrotnie sprawdzić, jeśli możesz, korzystając z ręcznej algebry liniowej. I pamiętajcie, że wagi powinny być odwrotnymi wariancjami.$\sigma$

EDYTOWAĆ

Jeśli często robisz tego rodzaju rzeczy, możesz rozważyć użycie ROOT(co wydaje się robić to natywnie, lma glmnie robić). Oto krótki przykład tego, jak to zrobić ROOT. Po pierwsze, ROOTmoże być używany przez C ++ lub Python, a jego ogromne pobieranie i instalacja. Możesz wypróbować go w przeglądarce za pomocą notatnika Jupiter, klikając link tutaj , wybierając „Binder” po prawej stronie i „Python” po lewej stronie.

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()

Wprowadziłem pierwiastki kwadratowe jako niepewności dotyczące wartości . Moc wyjściowa dopasowania to$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    

i powstaje ładna fabuła:

Instalator ROOT może również radzić sobie z niepewnościami wartości , co prawdopodobnie wymagałoby jeszcze większego włamania . Jeśli ktoś zna natywny sposób robienia tego w R, byłbym zainteresowany, aby się tego nauczyć.$x$lm

DRUGA EDYCJA

Druga odpowiedź z tego samego poprzedniego pytania autorstwa @Wolfgang daje jeszcze lepsze rozwiązanie: rmanarzędzie z metaforpakietu (pierwotnie zinterpretowałem tekst w tej odpowiedzi, aby nie obliczyć przechwytywania, ale tak nie jest). Przyjmując wariancje w pomiarach y po prostu 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

To zdecydowanie najlepsze czyste narzędzie R dla tego typu regresji, jakie znalazłem.