Lineares Modell, bei dem die Daten unsicher sind, unter Verwendung von R.

Angenommen, ich habe Daten mit Unsicherheiten. Beispielsweise:

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

Die Art der Unsicherheit kann beispielsweise Wiederholungsmessungen oder -experimente oder Messinstrumentenunsicherheit sein.

Ich möchte mit R eine Kurve daran anpassen, was ich normalerweise tun würde lm. Dies berücksichtigt jedoch nicht die Unsicherheit in den Daten, wenn es mir die Unsicherheit in den Anpassungskoeffizienten und folglich die Vorhersageintervalle gibt. In der Dokumentation finden Sie auf der lmSeite Folgendes:

... Gewichte können verwendet werden, um anzuzeigen, dass unterschiedliche Beobachtungen unterschiedliche Varianzen aufweisen ...

Ich denke also, dass das vielleicht etwas damit zu tun hat. Ich kenne die Theorie, es manuell zu machen, aber ich habe mich gefragt, ob es möglich ist, dies mit der lmFunktion zu tun . Wenn nicht, gibt es eine andere Funktion (oder ein anderes Paket), die dazu in der Lage ist?

BEARBEITEN

Nach einigen Kommentaren finden Sie hier einige Erläuterungen. Nehmen Sie dieses Beispiel:

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)

Gibt mir:

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

Grundsätzlich sind meine Koeffizienten a = 39,8 ± 22,3, b = 92,0 ± 9,3, c = -4,3 ± 0,8. Nehmen wir nun an, dass für jeden Datenpunkt der Fehler 20 ist. Ich werde ihn weights = rep(20,10)im lmAufruf verwenden und erhalte stattdessen Folgendes:

Residual standard error: 84.87 on 7 degrees of freedom

Die Standardfehler der Koeffizienten ändern sich jedoch nicht.

Manuell weiß ich, wie man das macht, indem man die Kovarianzmatrix unter Verwendung der Matrixalgebra berechnet und die Gewichte / Fehler dort einfügt und die Konfidenzintervalle daraus ableitet. Gibt es eine Möglichkeit, dies in der lm-Funktion selbst oder in einer anderen Funktion zu tun?

Diese Art von Modell ist in bestimmten Bereichen der Wissenschaft (z. B. Physik) und Ingenieurwissenschaften weitaus häufiger als die "normale" lineare Regression. In Physik-Tools wie ROOTist diese Art der Anpassung trivial, während die lineare Regression nicht nativ implementiert ist! Physiker neigen dazu, dies nur als "Anpassung" oder als Chi-Quadrat-Minimierungsanpassung zu bezeichnen.

$\sigma$

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

lm Gewichte und der Standardfehler

In den Antworten gibt es einige mögliche Lösungen. Insbesondere schlägt eine anonyme Antwort dort die Verwendung vor

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

lm$\sigma$

BEARBEITEN

Wenn Sie so etwas häufig tun, sollten Sie es in Betracht ziehen ROOT(was dies nativ zu tun scheint lmund glmnicht). Hier ist ein kurzes Beispiel dafür ROOT. Zunächst einmal ROOTkann über C ++ oder Python verwendet werden, und es ist ein riesiger Download und Installation. Sie können es im Browser mit einem Jupiter-Notizbuch versuchen, indem Sie dem Link hier folgen und rechts "Binder" und links "Python" auswählen.

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    

und eine schöne Handlung wird produziert:

$x$lm

ZWEITE BEARBEITUNG

Die andere Antwort aus derselben vorherigen Frage von @Wolfgang bietet eine noch bessere Lösung: das rmaTool aus dem metaforPaket (ich habe den Text in dieser Antwort ursprünglich so interpretiert, dass er den Achsenabschnitt nicht berechnet hat, aber das ist nicht der Fall). Nehmen Sie die Varianzen in den Messungen y einfach als 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

Dies ist definitiv das beste reine R-Tool für diese Art von Regression, die ich gefunden habe.