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

9

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?

Lista gimnastyczna
źródło
Jeśli znasz rozkład danych, możesz przeładować go za pomocą bootpakietu w R. Następnie możesz pozwolić regresji liniowej na zestaw danych ładowania początkowego.
Ferdi,
lmużyje znormalizowanych wariancji jako wag, a następnie przyjmie, że model jest statystycznie poprawny do oszacowania niepewności parametrów. Jeśli uważasz, że tak nie jest (paski błędów zbyt małe lub zbyt duże), nie powinieneś ufać żadnym szacunkom niepewności.
Pascal,
Zobacz także to pytanie tutaj: stats.stackexchange.com/questions/113987/...
jwimberley,

Odpowiedzi:

14

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ć σ

Lie12(yi(axi+b)σ)2
log(L)=constant12σ2i(yi(axi+b))2
σ
Le12(y(ax+b)σi)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)=constant12(yi(axi+b)σi)2
1/σi2log(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ólnyF=maF=ma+ϵlmσ2w 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.σ

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 toy

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:

quadfit

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ć.xlm

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.

jwimberley
źródło
Myślę, że cofnięcie skalowania jest zasadniczo błędne lm. Jeśli to zrobisz, statystyki sprawdzania poprawności, takie jak chi-kwadrat, zostaną wyłączone. Jeśli dyspersja twoich reszt nie pasuje do twoich słupków błędów, coś jest nie tak w modelu statystycznym (albo wybór modelu, albo słupki błędów albo normalna hipoteza ...). W obu przypadkach niepewności parametrów będą niewiarygodne !!!
Pascal
@PascalPERNOT Nie myślałem o tym; Pomyślę o twoich komentarzach. Szczerze mówiąc, zgadzam się w ogólnym sensie, że uważam, że najlepszym rozwiązaniem jest użycie oprogramowania fizyki lub oprogramowania inżynierskiego z gwarancją poprawnego rozwiązania tego problemu, a nie włamanie się w lmcelu uzyskania prawidłowego wyniku. (Jeśli ktoś jest ciekawy, pokażę, jak to zrobić ROOT).
jwimberley,
1
Jedną potencjalną zaletą podejścia statystycznego do problemu jest to, że pozwala on na łączenie oszacowań wariancji między obserwacjami na różnych poziomach. Jeśli leżąca u podstaw wariancja jest stała lub ma określony związek z pomiarami, jak w procesach Poissona, wówczas analiza zostanie zazwyczaj poprawiona w porównaniu z tym, co otrzymujesz z (zazwyczaj nierealistycznego) założenia, że ​​zmierzona wariancja dla każdego punktu danych jest poprawna, a zatem niesprawiedliwie ważona niektóre punkty danych. W danych PO sądzę, że założenie o stałej wariancji może być lepsze.
EdM,
1
@ jwimberley Zakładam, że zapewnia przeskalowanie błędu standardowego ważonych reszt do 1, przed obliczeniem macierzy kowariancji parametrów. Możesz to sprawdzić, mnożąc swoje wagi przez i zobaczyć, jak wpływa to na wyjście „Błąd standardowy resztek”. W twoim przykładzie zmienia się on z 1.088 na 1. Jeśli twoje ustawienie jest poprawne statystycznie, skalowanie ma jedynie niewielki wpływ na niepewności parametrów ...σσ2
Pascal,
1
Dobra dyskusja na te tematy znajduje się w rozdziale 8 Andreona, S. i Weavera, B. (2015) Bayesowskich metod nauk fizycznych. Skoczek. springer.com/us/book/9783319152868
Tony Ladson