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ę, lm
strona 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 lm
funkcji. 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 lm
wywoł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?
źródło
boot
pakietu w R. Następnie możesz pozwolić regresji liniowej na zestaw danych ładowania początkowego.lm
uż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.Odpowiedzi:
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
ROOT
dopasowanie 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ćσ
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ólnyfa= m a fa= m a + ϵ σ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
lm
lm
interpretacji 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σ
lm
skaluje 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,lm
aglm
nie robić). Oto krótki przykład tego, jak to zrobićROOT
. Po pierwsze,ROOT
moż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.Wprowadziłem pierwiastki kwadratowe jako niepewności dotyczące wartości . Moc wyjściowa dopasowania toy
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:
rma
narzędzie zmetafor
pakietu (pierwotnie zinterpretowałem tekst w tej odpowiedzi, aby nie obliczyć przechwytywania, ale tak nie jest). Przyjmując wariancje w pomiarach y po prostu y:To zdecydowanie najlepsze czyste narzędzie R dla tego typu regresji, jakie znalazłem.
źródło
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 !!!lm
celu uzyskania prawidłowego wyniku. (Jeśli ktoś jest ciekawy, pokażę, jak to zrobićROOT
).