Wie werden Gewichte für die WLS-Regression in R bestimmt?

8

Ich versuche, das Alter als Funktion einer Reihe von DNA-Methylierungsmarkern vorherzusagen. Diese Prädiktoren sind kontinuierlich zwischen 0 und 100. Wenn ich eine OLS-Regression durchführe, kann ich sehen, dass die Varianz mit dem Alter zunimmt.

Daher habe ich mich für ein gewichtetes Regressionsmodell entschieden. Ich habe jedoch Probleme bei der Entscheidung, wie die Gewichte für mein Modell definiert werden sollen. Ich habe die fGLS-Methode wie folgt verwendet:

OLSressq <- OLSres^2                 # Square residuals
lnOLSressq <- log(OLSressq)          # Take natural log of squared residuals
aux <- lm(lnOLSressq~X)              # Run auxillary model
ghat <- fitted(aux)                  # Predict g^
hhat <- exp(ghat)                    # Create h^
fGLS <- lm(Y~X, weights = 1/hhat)    # Weight is 1/h^

Und das waren meine Ergebnisse:

Call:
lm(formula = Y ~ X, weights = 1/hhat)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-4.9288 -1.2491 -0.1325  1.2626  5.1452 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 23.1009494  5.2299867   4.417 1.64e-05 ***
XASPA       -0.1441404  0.0474738  -3.036  0.00271 ** 
XPDE4C       0.6421385  0.0812891   7.899 1.83e-13 ***
XELOVL2     -0.2040382  0.0866564  -2.355  0.01951 *  
XELOVL2sq    0.0088532  0.0009381   9.438  < 2e-16 ***
XEDARADD    -0.1965472  0.0348989  -5.632 5.98e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 1.762 on 200 degrees of freedom
Multiple R-squared:  0.9687,    Adjusted R-squared:  0.9679 
F-statistic:  1239 on 5 and 200 DF,  p-value: < 2.2e-16

Bevor ich jedoch herausfand, wie man die fGLS-Methode durchführt, spielte ich mit verschiedenen Gewichten herum, nur um zu sehen, was passieren würde. Ich habe 1 / (quadratische Residuen des OLS-Modells) als Gewichte verwendet und am Ende Folgendes erhalten:

Call:
lm(formula = Y ~ X, weights = 1/OLSressq)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.0893 -0.9916 -0.7855  0.9998  2.0238 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 30.8756737  1.1355861   27.19   <2e-16 ***
XASPA       -0.1956188  0.0116329  -16.82   <2e-16 ***
XPDE4C       0.6168490  0.0102149   60.39   <2e-16 ***
XELOVL2     -0.1596969  0.0116723  -13.68   <2e-16 ***
XELOVL2sq    0.0078459  0.0001593   49.26   <2e-16 ***
XEDARADD    -0.2492048  0.0068751  -36.25   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 1 on 200 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.133e+06 on 5 and 200 DF,  p-value: < 2.2e-16

Da der verbleibende Standardfehler kleiner ist, R² gleich 1 ist (ist das überhaupt möglich?) Und die F-Statistik viel höher ist, bin ich versucht anzunehmen, dass dieses Modell besser ist als das, was ich mit der fGLS-Methode erreicht habe. Es scheint mir jedoch, dass die zufällige Auswahl von Gewichten durch Versuch und Irrtum immer schlechtere Ergebnisse liefern sollte, als wenn Sie tatsächlich mathematisch versuchen, die richtigen Gewichte zu schätzen.

Kann mir jemand einen Rat geben, welche Gewichte ich für mein Modell verwenden soll? Ich habe auch hier und da gelesen, dass Sie R² nicht so interpretieren können, wie Sie es bei der Durchführung einer OLS-Regression tun würden. Aber wie sollte es dann interpretiert werden und kann ich es trotzdem verwenden, um mein WLS-Modell irgendwie mit meinem OLS-Modell zu vergleichen?

I. Smeers
quelle
2
Da wäre ich sehr vorsichtig R-squared = 1. Haben Sie eine Vorstellung davon, wie hoch die möglichen Gewichte sein sollten? Klingt wie 1/(squared residuals of OLS model)war nur eine Vermutung. Wenn Sie über die Gewichte im Dunkeln sind, empfehle ich die Verwendung von GLS oder iterativ gewichteten kleinsten Quadraten.
Jon
Ich muss hinzufügen, dass beim Anpassen des gleichen Modells an einen Trainingssatz (die Hälfte meiner ursprünglichen Daten) das R-Quadrat von 1 auf 0,9983 gesunken ist. Es war in der Tat nur eine Vermutung, weshalb ich schließlich fGLS verwendet habe, wie oben beschrieben. Ist es das, was Sie mit "Ich schlage vor, GLS zu verwenden" meinen? Ich bin nur verwirrt darüber, warum das Modell, das ich durch einfaches Erraten der Gewichte erstellt habe, besser passt als das Modell, das ich durch Schätzen der Gewichte durch fGLS erstellt habe. Ich habe noch nichts von iterativ gewichteten kleinsten Quadraten gehört, aber ich werde es untersuchen. Vielen Dank.
I. Smeers
@Jon, machbare GLS erfordert die Angabe der Gewichte (während nicht realisierbare GLS, die theoretisch optimale Gewichte verwenden, kein machbarer Schätzer sind, dh in der Praxis nicht verwendet werden können).
Richard Hardy
Ja das ist richtig. Sie könnten jedoch die Korrelationsstruktur in der nlme::glsFunktion angeben . Das nlme::corClassesbietet eine Liste verschiedener Korrelationsstrukturen.
Jon

Antworten:

0

Warum benutzt du FLGS? Haben Sie Heteroskedastizität und Korrelation zwischen den Residuen? Und ist die Matrix-Var-Cov-Matrix unbekannt? Versuchen Sie es bptest(your_model)und wenn der p-Wert kleiner als das Alpha ist (z. B. 0,05), liegt eine Heteroskedastizität vor. Und dann sollten Sie versuchen zu verstehen, ob es eine Korrelation zwischen den Residuen mit einem Durbin Watson-Test gibt: dwtest(your_model)Wenn die Statistik W zwischen 1 und 3 liegt, gibt es keine Korrelation. Wenn Sie also nur Heteroskedastizität haben, sollten Sie WLS wie folgt verwenden:

mod_lin <- lm(Price~Weight+HP+Disp., data=df)
wts     <- 1/fitted( lm(abs(residuals(mod_lin))~fitted(mod_lin)) )^2
mod2    <- lm(Price~Weight+HP+Disp., data=df, weights=wts)

So mod2ist es auch mit dem alten Modell, jetzt mit WLS.

R-square = 1Es ist zu komisch. Vielleicht gibt es Kollinearität.

Lorenzo Famiglini
quelle
1
Warum sollte ein DW-Test angemessen sein? Ich denke, es wird nur für die Autokorrelation verwendet und ich sehe nicht, wie das in diesem Fall zutreffen würde.
meh
1
Willkommen bei xvalidated! Bitte angeben , aus dem Paket Funktionen bptestund dwtestkommen aus als sie sind nicht Teil der Standard - R - Verteilung.
Helix123
Weil Sie verstehen müssen, welcher Schätzer der beste ist: wie wls, fgls, ols ect ..
Lorenzo Famiglini