Segmentierte Regression
Neue Frage »

28.11.2011, 19:04 vogelstrauss Auf diesen Beitrag antworten »
Segmentierte Regression
Hallo,

ich habe leider einige Probleme mit der Interpretation der Ergebnisse meiner Regressionsanalyse, ich hoffe mir kann jemand helfen.
Ich analysisere Zeitreihendaten von zwei verschiedenen Gruppen (control und test) und zwei Interventionen innerhalb der Zeitreihe.
Meine abhaengige Variable heisst "Prices" und zu jedem Zeitpunkt T habe ich zwei Werte (control und test).

Mein Modell sieht so aus:



Die zwei interventionen a und b sind zum Zeitpunkt T=13 und T=21. Die variablen Da und Db sind "dummy"-variablen fuer die Interventionen, Da=1 wenn T>=13, Db=1 wenn T>=21, ansonsten Da,Db=0.
Die Variable G ist ebenfalls eine dummy variable, sie ist 0 fuer die Kontroll- und 1 fuer die Testgruppe.
Die Variablen Pa und Pb beschreiben die Zeit, die nach Intervention a bzw. b vergangen ist.

Dann habe ich mit R eine lineare Regression durchgefuehrt, mit forlgendem Output:

code: 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:

Call:
lm(formula = Acquisition.Prices ~ T + Da + Db + Pa + Pb + G + 
    GT + GDa + GDb + GPa + GPb, data = d)

Residuals:
     Min       1Q   Median       3Q      Max 
-2605.74   -21.02     3.01   124.61  2036.50 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2435.364    440.325   5.531 1.05e-06 ***
T             21.675     59.828   0.362    0.719    
Da            15.538    679.489   0.023    0.982    
Db             9.894    638.092   0.016    0.988    
Pa             9.325    125.565   0.074    0.941    
Pb            -7.304    125.565  -0.058    0.954    
G           2686.773    622.714   4.315 7.17e-05 ***
GT            41.612     84.610   0.492    0.625    
GDa          136.385    960.942   0.142    0.888    
GDb         -371.773    902.398  -0.412    0.682    
GPa          169.388    177.576   0.954    0.345    
GPb          -53.843    177.576  -0.303    0.763    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 715.4 on 52 degrees of freedom
Multiple R-squared: 0.9327,     Adjusted R-squared: 0.9185 
F-statistic: 65.51 on 11 and 52 DF,  p-value: < 2.2e-16


Wenn ich diese Koeffizienten benutze, kann ich meine Daten ganz gut beschreiben (siehe Plot im Dateianhang).

Mein Problem ist: Die errechneten Koeffizienten sind alle nicht signifikant (bis auf den y-intercept fuer test und control)!!
Die Koeffizienten fuer Da und Db zum Beispiel muessten mir den Anstieg in Price gleich nach Interventionen a und b geben. Die Koeffizienten fuer Pa und Pb sollten mir die Aenderung im Trend der Kurve nach Interventionen a und b geben... Ich kann mir nur nicht vorstellen, dass das alles nicht signifikant ist!
Wenn Fast alle Variablen nicht signifikant sind, dann muesste ich sie ja aus dem Modell entfernen koennen und trotzdem meine Daten beschreiben koennen.

Wenn ich mir die Vorhersage des Modells anschaue (siehe Anhang, blau=control, rot=test), dann sieht mir die Aenderung nach den Interventionen schon signifikant aus!

Koennte jemand mir helfen zu verstehen, wie eigentlich die P-Werte fuer die Signifikanz errechnet werden? Habe ich etwas vergessen oder ist vielleicht mein Modell falsch?

Entschuldigung fuer den langen Beitrag, aber ich hoffe mein Problem ist verstaendlich und ich hoffe mir kann jemand helfen!

Vielen Dank im Voraus,

vogelstrauss
28.11.2011, 19:33 Black Auf diesen Beitrag antworten »

Die p Werte für die Signifikanz ergeben sich jeweils durch einen Test auf

Der Output macht mir den Anschein dass du viel zu viele abhängige Variablen drin hast.
Du kannst versuchen per gemischter Rückwärtsinduktion ein reduziertes Modell zu bestimmen, R hat dafür glaub ich sogar eine vorgefertigte Funktion.
29.11.2011, 14:29 vogelstrauss Auf diesen Beitrag antworten »

Hallo,

danke fuer die schnelle Antwort!

Es scheinen in der Tat sehr viele Variablen in dem Modell, die nicht benoetigt werden um meine Daten zu beschreiben...

Ich habe leider noch keine R-Funktion fuer die Rueckwaertsinduktion gefunden,
aber ich konnte mein Modell 'von Hand' reduzieren, indem ich einfach jeweils die Variablen mit den hoechsten P-Werten aus dem Modell entfernt habe.

Jetzt sieht das Modell und der Ouput so aus:

code: 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:

Call:
lm(formula = Acquisition.Prices ~ T + G + GPa + GPb, data = d)

Residuals:
     Min       1Q   Median       3Q      Max 
-2678.40   -39.95     3.82    65.50  2202.99 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2352.87     240.00   9.804 5.37e-14 ***
T              30.12      12.59   2.393   0.0199 *  
G            3045.89     247.65  12.299  < 2e-16 ***
GPa           233.25      49.95   4.670 1.79e-05 ***
PGb          -113.66      81.71  -1.391   0.1694    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 680.2 on 59 degrees of freedom
Multiple R-squared: 0.931,      Adjusted R-squared: 0.9263 
F-statistic: 198.9 on 4 and 59 DF,  p-value: < 2.2e-16


Der R-squared Wert ist aehnlich bei diesem Modell, jedoch sind die Variablen signifikanter.

Vielen Dank,

vogelstrauss
29.11.2011, 15:52 Black Auf diesen Beitrag antworten »

Einfach die Variablen mit dem höchsten p-Wert rauszuwerfen ist im Allgemeinen keine gute Idee.
Besser ist man verwendet den AIC um die Modelle zu vergleichen, genau das macht man bei der Rückwärtsinduktion.
In R gibt es dafür die Funktion step(), (step direkt auf das lm Objekt mit allen Variablen anwenden) diese wirft so lange Variablen raus bis sich der AIC nicht mehr verbessert.
29.11.2011, 19:54 vogelstrauss Auf diesen Beitrag antworten »

Mit der 'step' methode bekomme ich in der Tat ein anderes Modell:
T + G + GT + GPa
mit einem AIC von 839.39 und R-quadrat von 0.9313.

Zum Vergleich, mein von Hand reduziertes Modell hat einen AIC von 839.67 und R-quadrat: 0.93.

Das Probelm ist, dass mit diesem Modell die Aenderungen im Slope nach den Interventionen (z.B. beschrieben durch GPa) nicht signifikant sind:

code: 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:

Call:
lm(formula = Acquisition.Prices ~ T + G + GT + GPa, data = d)

Residuals:
     Min       1Q   Median       3Q      Max 
-2645.27   -21.83     2.68    86.49  2228.17 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2405.92     245.72   9.791 5.63e-14 ***
T              26.90      13.00   2.070   0.0428 *  
G            2582.44     458.08   5.638 5.12e-07 ***
GT             67.26      45.34   1.483   0.1433    
GPa            92.56      58.75   1.576   0.1205    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 678.8 on 59 degrees of freedom
Multiple R-squared: 0.9313,     Adjusted R-squared: 0.9266 
F-statistic: 199.8 on 4 and 59 DF,  p-value: < 2.2e-16


War das andere Modell (T + G + GPa + GPb) wirklich so falsch?
Dieses Modell zeigt, was ich eigentlich urspruenglich sehen wollte, naemlich dass GPa signifikant ist.

Heisst das nun, dass die Aenderungen in meinen Daten mit den gegebenen Interventionen wirklich nicht signifikant sind?
Bin ich damit am Ende und muss mir eingestehen, dass die Daten schlecht sind?

Oder kann ich noch andere Sachen probieren z. B. Outlier korrigieren oder fuer eventuelle Korrelationen der Error- Terme korrigieren?

Vielen Dank,

vogelstrauss
30.11.2011, 17:07 Black Auf diesen Beitrag antworten »

Also optimal sind sicher beide Modelle nicht.
Könntest du vllt mal den exakten R Befehl mit dem du die Regression erstellt hast posten?

Es gibt auch noch einige andere Methoden mit denen man reduzierte Modelle bestimmen kann, sagt dir resettest etwas?

Wie ich gerade erfahren musste macht R bei step() nicht automatisch gemischte Rückwärtsinduktion, schau mal ob er dir mit step(deinlm,direction="both") was anderes ausspuckt.

Edit: Ich glaube nicht dass Ausreißer das Problem verursachen, dazu müsste das Bestimmtheitsmaß deutlich schlechter sein.
 Anzeige 
 
30.11.2011, 17:36 vogelstrauss Auf diesen Beitrag antworten »

Wenn ich step(deinlm,direction="both") verwende, bekomme ich das gleiche Ergebnis wie nur mit step(meinlm).

Wenn ich das richtig vertehe, testet resettest auf H0: Modell ist linear. Mir ist leider noch nicht klar, was das argument 'type' zu bedeuten hat und welche Potenzen ich fuer das Argument 'power' nehmen sollte...

Hier sind meine R-commands fuer die Regression:

code: 
1:
2:
3:
4:
5:

table = read.table('data.txt', sep=';', header=TRUE)
d <- as.data.frame(table)
o = lm(Acquisition.Prices ~ T + Da + Db + Pa + Pb + G + GT + GDa + GDb + GPa + GPb, data=d)
summary(o)
step(o)

Die Daten habe ich angehaengt.

glm() gibt uebrigens das selbe Ergebnis.

Vielen Dank fuer Deine Hilfe!!
30.11.2011, 22:45 Black Auf diesen Beitrag antworten »

So hab mir jetzt das Ganz mal genauer angeschaut.

Das Hauptproblem ist dass GDa, GDb, GPa, GPb allesamt extrem korreliert sind.
Ein zweites Problem ist dass du im Verhältnis zu der dir vorliegenden Datenmenge viel zu viele Variablen gewählt hast.
Zu der inhaltlichen Interpretation der Variablen kann ich natürlich nichts sagen, aber du hast dir bei der Wahl der Variablen bestimmt etwas gedacht. Es ist möglich dass viele der gewählten Variablen in Wahrheit tatsächlich einen signifikanten Einfluss haben, aber um das stichhaltig zu belegen reicht hier die Datenmenge nicht.


Tatsächlich wird hier nämlich die Anpassung im Wesentlichen von 2 Variablen beinflusst, nämlich von G und GPa:

code: 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:

Call:
lm(formula = Acquisition.Prices ~ GPa + G, data = d)

Residuals:
     Min       1Q   Median       3Q      Max 
-2654.38  -276.77    -2.13   256.02  2465.08 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2849.81     126.67   22.50   <2e-16 ***
GPa           214.05      18.55   11.54   <2e-16 ***
G            2894.85     216.60   13.37   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 716.6 on 61 degrees of freedom
Multiple R-squared: 0.9208,     Adjusted R-squared: 0.9182 
F-statistic: 354.6 on 2 and 61 DF,  p-value: < 2.2e-16


Man könnte GPa in obiger Regression auch durch GPb,GDa oder GDb ersetzen - die Anpassung wäre immernoch gut.

Das obige Modell ist das einzige das ich persönlich als zufriedenstellend bezeichnen würde, ob es inhaltlich sinnvoll ist musst du entscheiden.


Falls das nicht der Fall sein sollte, würde ich an deiner Stelle nochmal die Variablen (insbesondere die Dummies) überdenken: Wirf nur mal einen Blick auf die Spalten Da, Pa, GDa, GPa, GDb, GPb deiner Modellmatrix, dann siehst du dass du dir einen fast perfekten (und mehrfachen) Spaltendefekt konstruiert hast.
06.12.2011, 18:58 vogelstrauss Auf diesen Beitrag antworten »

Hallo,

ich hatte ein paar Tage keine Zeit mich mit der Regression zu befassen, hab jetzt aber noch ein paar Fragen...

Zitat:
im Verhältnis zu der dir vorliegenden Datenmenge viel zu viele Variablen gewählt hast.


Das sehe ich genauso, daher versuche ich immer noch, das Modell zu reduzieren.

Sehr wichtige Variablen in meinem Modell sind GPa und GDa, welche den Slope und die Aenderungen des Preises nach der ersten Intervention beschreiben. Diese Beiden Variablen sollten eigentlich signifikant sein. Deshalb irritiert es mich, dass man GPa und GDa im reduzierten Modell austauschen kann...

Ich habe noch ein wenig weitergeschaut, warum mein Modell nicht so ganz gut ist und habe gemerkt, dass ich keine Saisonbereinigung auf meinen Daten gemacht habe. Die Daten sind halbjaehrlich gesampled, deshalb koennten vielleicht versteckte Muster in den Daten sein.

Deshalb habe ich mir die Reste vom fit angeschaut und die Autokorrelationsfunktion der Reste berechnet:

code: 
1:
2:
3:
4:

a = acf(o$residuals)
par(mfrow=c(1,2))
plot(o$residuals, ylab="Residual value", xlab="Time", main="Residuals", type='b')
plot(a$acf, ylab="Autocorrelation", xlab="lag", main="Residual autocorrelation", type='b')


Der Plot ist im Anhang.
Beim Plot der Reste (links) kann ich kein besonderes Saison-Muster erkennen.

Trotzdem sind manche Autokorrelations-Werte recht hoch...

Ich habe hier
http://www.duke.edu/~rnau/testing.htm
gelesen, dass die Autokorellations-Werte ungefaehr bei sein sollten (0+-0.25 bei meinen Daten).

Wie man im Plot jedoch sieht sind einige Werte bei ungefaehr -0.5.

Koennte das ein Grund sein, warum mein Modell nicht gut ist?

Auf der selben site wurde gesagt, dass man mit einer weiteren Dummy-Variable fuer die Saison korrigieren kann.
(Ich habe in den Daten eine weitere Spalte 'Summer' angehaengt, die 1 ist wenn im Sommer gesampled werden und die Variable 'Summer' im Modell beruecksichtigt... )
Dadurch das mein Modell schon recht viele Variablen hat, hat das jedoch nichts gebracht.

Koennte es troztdem mit einer moeglichen Saisonabhaengigkeit in meinen Daten zusammenhaengen, dass ich kein gutes Modell hab?


Vielen Dank fuer die viele Hilfe,

Gruss, vogelstrauss
Neue Frage »
Antworten »



Verwandte Themen

Die Beliebtesten »
Die Größten »
Die Neuesten »


MatheBoard » Hochschulmathematik » Stochastik & Kombinatorik » Segmentierte Regression