FMS012/MASB03: M ATEMATISK STATISTIK 9 HP , HT-16
Datorövning 4: Regression
Syftet med den här laborationen är att du skall bli mer förtrogen med
• Enkel linjär regression
• Multipel linjär regression
Specialrutiner och datamaterial finns att hämta på kursens hemsida:
http://www.maths.lu.se/kurshemsida/fms012masb03/
1 Förberedelseuppgifter
1. Läs igenom denna handledning och lös uppgift ST35 i övningshäftet.
2. Förvissa dig om att du förstår hur regression fungerar, speciellt skillnaden mellan konfidensin- tervall för linjen och prediktionsintervall.
3. Redovisas vid laborationens start!
(a) Om vi skattat lutningen β i en enkel linjär regression till β ∗ = 0.75 och beräknat ett 95 % konfidensintervall I β = (−0.95, 2.45), kan vi då påstå att det finns ett statistiskt signifikant linjärt samband mellan x och y?
(b) Skriv upp hur residualerna i en enkel linjär regression beräknas.
(c) Hur skattar man σ med hjälp av residualerna?
(d) Hur bör residualerna se ut, enligt modellen? Hur kan man kontrollera det?
2 Enkel linjär regression
Vid enkel linjär regression söker man anpassa en rät linje till datamaterialet, dvs modellen är y i = α + β x i + ε i , i = 1, . . . , n,
där ε i är oberoende likafördelade störningar med väntevärdet 0 och variansen σ 2 .
Vi kommer i den följande framställningen att arbeta med matrisformuleringen av modellen, Y = Xβ + ε,
där de ingående matriserna har följande form:
Y =
y 1
y 2 .. . y n
, X =
1 x 1
1 x 2 .. . .. . 1 x n
, β =
α β
och ε =
ε 1 ε 2 .. . ε n
.
Som uppvärmning använder vi MATLAB -funktionen regress för att lösa uppgift ST35:
>> help regress
>> x=[1:8]’; % En kolumn med x-värden.
>> y=[1.5 2.3 1.7 2.0 2.5 1.9 2.2 2.4]’; % En kolumn med y-värden.
>> X=[ones(size(x)) x]; % En kolumn med ettor och en med
% x-värdena.
>> [b,bint]=regress(y,X) % b = skattade beta-parametrar,
% bint = konf.int. för beta-parametrarna.
>> mu=X*b; % Den skattade linjen.
>> plot(x,y,’*’,x,mu,’-’) % Rita obervationer och skattad linje.
>> xlabel(’x: inställning’)
>> ylabel(’y: dimension’)
Uppgift: Vad blev α ∗ och β ∗ ? Jämför β ∗ med din lösning.
Uppgift: Vad blev I α och I β ? Jämför I β med din lösning.
På kurshemsidan finns en specialskriven funktion, reggui.m (med ett antal hjälpfiler). Den både skattar modellparametrarna, och ritar upp data, skattad linje, residualer och konfidens- och predik- tionsintervall. Spara ner den, och dess hjälpfiler, från kurshemsidan och lös uppgift ST35 igen:
>> help reggui
>> reggui(x,y)
Uppgift: Identifiera α ∗ och β ∗ i figuren och jämför med dina tidigare beräkningar.
Uppgift: Identifiera I α och I β i figuren och jämför med dina tidigare beräkningar.
Uppgift: Den inritade skattade linjen är röd. Det betyder att minst en av parametrarna inte är signifikant skild från noll. Vilken/vilka?
Uppgift: Vi kan passa på att titta på residualerna också, dvs e i = y i − α ∗ − β ∗ x i . I den nedre vänstra figuren är e i ritade mot x i . Tänk efter att det stämmer. Ser residualerna ut som de ska?
Residualerna ska vara normalfördelade också. I den nedre högra figuren är de utritade i ett normal- fördelningspapper. På x-axeln finns residualerna i storleksordning. På y-axeln finns deras empiriska fördelningsfunktion, se laboration 1, men skalan är satt så att en normalfördelning blir en rät linje medan alla andra fördelningar blir krokiga.
Uppgift: Ser det ut som om residualerna kan vara normalfördelade?
3 Kalibrering av flödesmätare
Vi ska nu studera ett lite större problem med hjälp av reggui.
Bakgrund
Kalibrering av en flödesmätare genomförs oftast i en speciell kalibreringsrigg. Här finns en refe- rensmätare eller referensmetod för att mäta flödet. För att erhålla en god bild av hur den testade flödesmätaren fungerar utförs kalibreringen vid ett stort antal flöden. Tyvärr kan man även vid kalibrering råka ut för situationer där den testade mätaren störs av testförhållandena.
Om t.ex. pulsationer uppträder i flödet kommer detta att negativt påverka resultaten för den testade mätaren. Detta visar sig oftast vid låga flödeshastigheter, då ultraljudsmätare tenderar att överskatta flödeshastigheten. Detta orsakas av att vi erhåller en laminär flödesprofil i röret, vilket medför att en ultraljudsmätare kan överskatta flödet med upp till 33 % vid fullt utbildad laminär strömning.
Vid låga flöden ser vi även att vi har stora fluktuationer i resultaten. Detta beror troligen på att vi har flödespulsationer i flödesriggen vilka kommer att orsaka fluktuerande resultat för ultraljuds- flödesmätaren, bland annat orsakat av s.k. aliasproblem (d.v.s. mätsystemet arbetar med en för låg sampelfrekvens i förhållande till frekvenserna hos det uppmätta).
Vid höga flöden uppträder troligen kavitation (ett slags bubbelbildning) inne i ultraljudsflödesmä- taren vilket kan förklara de positiva felen och den ökade spridningen för strömningshastigheter över 6.3 m/s.
Metod
Vi har nu tillgång till data från en kalibrering av en ultraljudsflödesmätare. Datamaterialet, som kommer från institutionen för värme- och kraftteknik, omfattar 71 mätningar och är lagrat i matri- sen flow, där varje rad innehåller data från en mätning, variabeln fx2 avser referensflödesmätningar från kalibreringsriggen och fy2 avser respektive flöden uppmätta med den testade ultraljudsflödes- mätaren (flödeshastigheterna givna i enheten m/s).
Den använda kalibreringsriggen använder kontinuerlig vägning av det genomströmmande vattnet för att bestämma ett massflöde som sedan kan räknas om till medelhastighet i röret, vilket är vad ultraljudsmätaren mäter.
Tanken är här att vi med hjälp av de gjorda mätningarna med givare och referens skall skatta para- metrarna i en enkel linjär regressionsmodell. Vi antar då att referensmätningarnas fel kan försummas i jämförelse med ultraljudsgivarens (varför måste vi bekymra oss om detta?) och att ultraljudsgiva- rens fel är oberoende, likafördelade och har väntevärdet noll.
För att studera detta datamaterial ska vi använda funktionen reggui vars finesser du förhoppnings- vis redan bekantat dig med. Observera att du t.ex. automatiskt kan rita ut konfidensintervall och prediktionsintervall genom att markera i tillämplig ruta. Genom att kryssa i ”Mark ints” och sedan klicka (och dra) i figuren kan man låta reggui beräkna och rita intervallen för ett visst x 0 -värde.
För att bilden skall bli tydligare börjar vi med att studera en liten delmängd av materialet, 10 talpar av flödesmätningar som ges i variablerna fx1 och fy1.
>> load flow.mat
>> reggui(fx1,fy1)
Använd nu funktionen interaktivt för att göra följande beräkningar:
Uppgift: Avläs det förväntade värdet enligt ultraljudsmätaren, då flödet enligt kalibreringsriggen är 0.40 m/s.
Uppgift: Avläs ett 95 %-igt konfidensintervall för detta förväntade värde.
Uppgift: Avläs ett 95 %-igt prediktionsintervall för en framtida observation från ultraljudsmätaren, då kalibreringsriggen ger mätvärdet 0.40 m/s.
Uppgift: Tänk efter vad det är som skiljer de två intervallen åt.
När vi sedan skall använda den kalibrerade ultraljudsmätaren innebär det i princip att vi ”läser baklänges” i kalibreringskurvan. Antag att vi med ultraljudsmätaren får mätvärdet 0.48 m/s.
Uppgift: Avläs ett 95 %-igt kalibreringsintervall för den ”sanna” flödeshastigheten (dvs det värde som kalibreringsriggen skulle ge). Identifiera i figuren de kurvor som används vid den grafiska bestämningen av detta konfidensintervall och tänk efter varför det är just dem, man skall använda.
Uppgift: I vilka av intervallen ovan (konfidensintervallet för linjen, prediktionsintervallet oxh ka- libreringsintervallet) är det viktigt att mätfelen är normalfördelade och har konstant varians?
Dvs, vilka av intervallen har en bredd som inte kan göras hur liten som helst genom att öka antalet observationer?
Intervallen, särskilt kalibreringsintervallet, är oanvändbart breda men 10 observationer är för litet för att kunna ge bra skattningar. Låt oss istället använda hela datamaterialet:
>> reggui(fx2,fy2)
Upppgift: Avläs det förväntade värdet enligt ultraljudsmätaren, då flödet enligt kalibreringsriggen är 0.40 m/s.
Uppgift: Avläs ett 95 %-igt konfidensintervall för detta förväntade värde.
Uppgift: Avläs ett 95 %-igt prediktionsintervall för en framtida observation från ultraljudsmätaren, då kalibreringsriggen ger mätvärdet 0.40 m/s.
Uppgift: Antag att vi med ultraljudsmätaren får mätvärdet 0.48 m/s och avläs ett 95 %-igt kalibre- ringsintervall för den ”sanna” flödeshastigheten (dvs det värde som kalibreringsriggen skulle ge).
Uppgift: Jämför intervallbredderna baserade på de 10 mätningarna med motsvarande intervall-
bredder för den modell som är anpassad till alla de 71 mätpunkterna. Skiljer de två konfi-
densintervallen väsentligt åt (eller inte)? Hur kan det förklaras?
Uppgift: Jämför de två prediktionsintervallen. Skiljer de sig väsentligt åt (eller inte)? Hur kan det förklaras?
Uppgift: Jämför de två kalibreringsintervallen. Skiljer de sig väsentligt åt (eller inte)? Hur kan det förklaras?
Uppgift: Innan vi törs använda den skattade regressionslinjen för prediktion, måste vi naturligtvis förvissa oss om att modellen är adekvat. Ger residualplottarna anledning att förkasta mo- dellen eller anser du att du på goda grunder kan använda den skattade regressionslinjen för kalibrering av ultraljudsmätaren?
4 Polynomregression
Datamaterialet som du skall arbeta med i detta avsnitt är koldioxidhalter uppmätta över en vulkan varje månad under en period av 32 år, dvs totalt finns 32 · 12 = 384 mätvärden. Materialet finns i filen co2.mat. Mätvärdena ligger i variabeln co2data. Plotta mätvärdena:
>> load co2
>> plot(co2data)
Det finns uppenbarligen en kraftig periodicitet (årsvariation) i mätningarna, och en sådan låter sig inte så lätt fångas med en polynomiell regressionsfunktion. Detta problem kan lösas på flera sätt.
Ett är att införa en sinus-funktion som modellerar variationen, ett annat är att differentiera datase- kvensen, dvs undersöka z i = y i − y i−1 i stället för y-värdena själva. Vi skall dock välja den mycket enkla lösningen att medelvärdesbilda över varje år. Årsmedelvärdena finns i variabeln co2medel. Vi skapar även en vektor med den förklarande variabeln (årtalet, räknat från lämplig nollpunkt).
>> x=(1:32)’;
>> plot(x,co2medel,’*’)
Uppenbarligen är den periodiska variationen borta, vilket också var syftet med medelvärdesbild- ningen. Vi skall nu göra polynomregression på materialet, dvs vår modell är
y i = β 0 + β 1 x i + β 2 x i 2 + . . . + β k x i k + ε i , i = 1, . . . , n,
där ε i är oberoende likafördelade störningar med väntevärdet 0 och variansen σ 2 . Som modellen är skriven ovan är den olinjär, ty ett polynom är inte en linjär funktion, men vi kan göra den linjär genom att införa de nya förklarande variablerna x ji = x i j för j = 1, . . . , k, i = 1, . . . , n, och skriva
y i = β 0 + β 1 x 1i + β 2 x 2i + . . . + β k x ki + ε i , i = 1, . . . , n.
Vi börjar med att anpassa en linjär funktion till datamaterialet, dvs polynomets ordningsgrad k = 1.
Skattningarna av β 0 och β 1 erhålles med hjälp av funktionen reggui:
>> reggui(x,co2medel)
Detta är den modell vi skall arbeta vidare med.
Uppgift: Verkar en rät linje vara en tillfredsställande regressionsmodell?
Diagrammet visar att residualerna i mitten av mätserien tycks komma från en annan fördelning är residualerna i början och slutet av densamma. Alternativt finns en stark korrelation mellan stör- ningarna vilket strider mot oberoendeantagandet. Vi drar alltså slutsatsen att en enkel linjär regres- sionsmodell inte passar det aktuella datamaterialet. Nästa steg är att försöka anpassa en kvadratisk funktion till mätvärdena, dvs vi använder ordningstalet k = 2 för regressionspolynomet (använd knappen ”degree” i reggui).
Uppgift: Verkar den kvadratiska modellen vara bättre än den linjära? Kan residualerna tänkas kom- ma från samma fördelning? Kan de tänkas vara normalfördelade?
Uppgift: Avsluta med att studera de 95 %-iga konfidensintervallen för β 0 , β 1 och β 2 . Är β 2 signi- fikant skild från 0, dvs, om H 0 : β 2 = 0 och H 1 : β 2 6= 0, kan vi då förkasta H 0 (på nivån 5 %)? I så fall kan vi med gott samvete anta den kvadratiska modellen före den linjära.
Uppgift: På samma sätt kan man gå vidare och testa om en tredjegradsterm i regressionsfunktionen är relevant. Undersök de olika möjligheterna reggui ger dig att studera en regressionsmodell och välj olika gradtal i modellen.
Fick du några varningsmeddelanden i kommandofönstret? Vad kan det i så fall bero på?
Uppgift: Gör en bedömning av figurerna och utskriften med de skattade parametrarna och konfidens- intervallen och avgör vilken polynommodell som är mest adekvat.
5 Multipel regression — Frost
Hur beror antalet frostdagar i en ort på höjd och latitud? Om man känner höjden och latituden hos en ort kan man då förutsäga (prediktera) antalet frostdagar? Man noterade det genomsnittliga antalet frostdagar vid 20 olika väderstationer i West Virginia. Detta tillsammans med höjden över havet (feet) och stationens latitud finns i filen frost.mat. Väderstationernas namn finns i variabeln namn, höjden över havet i x1, latituden i x2 och antalet frostdagar i y.
5.1 Vilken typ av modell?
För att få en första uppfattning om eventuella samband, rita punktdiagram över den beroende variabeln gentemot de oberoende variablerna var och en för sig, dvs rita dels Frostduration mot Höjd över havet och dels Frostduration mot Latitud.
>> load frost
>> figure(1)
>> plot(x1,y,’*’)
>> xlabel(’x1: höjd över havet’)
>> ylabel(’y: antal frostdagar’)
>> figure(2)
>> plot(x2,y,’*’)
>> xlabel(’x2: latitud’)
>> ylabel(’y: antal frostdagar’)
Uppgift: Ger bilderna någon antydan om att ett linjärt samband skulle kunna råda mellan förkla- rande och beroende variabler?
Använd reggui för att göra två enkla linjära modeller, en där frostdagarna beror bara på höjden över havet och en där de istället beror på latituden.
Uppgift: Om man bara tänker använda en av x-variablerna, vilken ska det vara? Verkar det rimligt att latituden inte verkar spela någon roll?
I en enkel linjär regression kan man inte ta hänsyn till flera faktorer samtidigt. Man riskerar då att missa inverkan av vissa x-variabler när de drunknar i ett mycket starkare samband med andra x-variabler. I multipel regression tar man hänsyn till alla x-variablerna samtidigt och har större möjligheter att få med även svagare, men viktiga, samband. 1 Då undersöker man hur k uppmätta variabler, x 1 , . . . , x k , gemensamt påverkar en responsvariabel, y. I vårt problem har vi två förklarande variabler, dvs k = 2. Regressionsmodellen är alltså av formen
y i = β 0 + β 1 x 1i + β 2 x 2i + ε i , i = 1, . . . , n
där man tänker sig att alla ε i är oberoende och normalfördelade med väntevärde 0 och varians σ 2 . Uppgift: Identifiera de olika variablerna i modellen ovan i vårt frostproblem. Vad är alltså er ansatta
modell i just detta exemplet?
5.2 Regression med regress
Reggui är en specialskriven funktion för våra grundkurser i matematisk statistik. Den fungerar emel- lertid bara för enkel linjär regression och polynomregression, den går alltså inte att använda här.
Om man vill arbeta med multipel linjär regression (flera x-variabler) måste man använda Matlabs inbyggda funktion för regressionsanalys, regress. Börja med att bygga upp matrisen X som, i det här fallet, är en (20×3)-matris (det finns 20 st observationer av x1-värden respektive x2-värden) med första kolumnen enbart ettor, andra kolumnen bestående av x1-värdena och tredje av x2-värdena.
Vektorn b ger skattningarna av parametrarna β 0 , β 1 och β 2 medan deras konfidensintervall, med konfidensgraden 95 %, finns i matrisen Ib. Vektorn r ger residualerna.
>> X = [ones(size(x1)) x1 x2]
>> [b Ib r] = regress(y,X)
Uppgift: Vad blev skattningarna av β 0 , β 1 och β 2 ? Ange också motsvarande tre konfidensintervall.
Uppgift: Vilka av modellparametrarna är signifikant skilda från noll (på 5 %-nivån)? Kan vi för- enkla modellen och ta bort någon av variablerna?
Uppgift: Stämmer detta med era slutsatser från de två enkla regressionerna? Försök förklara skill- naden.
1
