1. No category

Matematikcentrum
Matematisk Statistik
Lunds Universitet
MASB11
VT15, lp3
Laboration 5
Regressionsanalys
2015-03-13
Syftet med laborationen är att vi skall bekanta oss med lite av de funktioner som finns i R för att utföra
korrelations- och regressionsanalys. När ni arbetat er igenom laborationshandledningen kan ni fortsätta
med projektet.
Introduktion – Regressionsanalys i R
Genom hela introduktionen ställs frågor i anslutning till analyserna. Kortfattade svar på en del av
dessa finns i slutet på introduktionen. Vi skall börja med att göra en regressionsanalys på följande
datamaterial:
Datamaterial: För 6 torskar har vi värden på variablerna Längd (cm) och Ålder (år).
Längd (cm)
Ålder (år)
Längd(cm)
Ålder(år)
15
1
60
6
30
2
58
5
35
3
25
2
50
4
12
1
55
5
43
4
Kan vi påvisa något samband (linjärt) mellan längd och ålder? Går det alltså avgöra åldern med hjälp
av fiskens längd.
Inläsning av data.
Börja med att mata in data till R. Lägg in värdena i en data.frame, som du döper till torskar, med
värdena i två kolumner: Längd och Ålder.
torskar <- data.frame(Längd=c(15,30,35,50,55,60,58,25,12,43), Ålder=c(1,2,3,4,5,6,5,2,1,4))
Datamaterialet skall alltså innehålla 2 kolumner med 6 värden i varje kolumn.
a) Beskrivning av data.
Börja nu med en grafisk beskrivning av sambandet genom att rita ett spridningsdiagram med Ålder på
x-axeln och Längd på y-axeln:
plot(torskar$Ålder, torskar$Längd)
Ser det ut som om det kan finnas ett linjärt samband?
b) Korrelationer.
Vi börjar med att beskriva sambandet mellan variablerna med hjälp av korrelationskoefficienten.
Beräkna Pearson's korrelationskoefficient med R's inbyggda funktion cor(...) och testa om den är
skild från noll med cor.test(...):
cor(torskar$Ålder, torskar$Längd)
cor.test(torskar$Ålder, torskar$Längd)
Tyder resultaten på att det finns något linjärt samband mellan Längd och Ålder?
c) Enkel linjär regression. Vi skall nu undersöka hur sambandet mellan variablerna ser ut genom
att anpassa en rät linje till data. Kommandot lm(y ~ x) (lm står för linear model) anpassar en linjär
modell för den beroende variabeln y som funktion av en eller flera förklarande variabler x. Sedan kan
vi få ut olika egenskaper hos modellen och skattningarna med ytterligare kommandon:
modell <- lm(Längd ~ Ålder, data=torskar)
modell # skattningarna av beta0 och beta1
summary(modell) # ger mer information, t.ex. signifikanser för skattnigarna
confint(modell) # Konfidensintervallet för skattnigarna av beta0 och beta1
Genomför analysen. Identifiera följande mått i utskriften: r – korrelationskoefficienten, r 2 –
förklaringsgraden, s – residualspridningen, de skattade koefficienterna med standardfel, t-test och
konfidensintervall.
För att få den skattade regressionslinjen utritad i figuren ni skapade i 1 a) kan du använda kommandot:
abline(modell)
d) Prognoser och konfidensintervall. Om man vill använda sin regressionsmodell för att göra
prognoser så kan detta enkelt göras efter att man skattat modellen. I R kan man prediktera i
samma datamaterial som man använde för att skatta modellen. Då får man en prediktion för varje
individ=rad. Man kan också mycket enkelt ange en helt annan uppsättning individer som man
vill prediktera för istället. Det är praktiskt när man t.ex. vill rita ut konfidensintervall och
prediktionsintervall snyggt.
# vi vill prediktera för en ålderssekvens i steg om halvår: 0.5, 1.0, 1.5, …, 7.0, 7.5:
x0 <- data.frame(Ålder=seq(0.5,7.5,0.5))
mu0konf <- predict(modell, x0, interval=”confidence”) # konfidensintervall
mu0pred <- predict(modell, x0, interval=”prediction”) # prediktionsintervall
cbind(x0, mu0pred)
Vad blir prognosen för längden för en sju år gammal torsk? Vad blir prognosintervallet?
För att få intervallen utritade i figur ritar vi linjer med våra prediktionsåldrar på x-axeln och
tillhörande intervallgränser på y-axeln. Vi vill dessutom rita konfidensintervallet som streckade
blå linjer och prediktionsintervallet som prickade röda:
lines(x0$Ålder, mu0konf[,”lwr”], col=”blue”, lty=2) # undre (lower) gränsen
lines(x0$Ålder, mu0konf[,”upr”], col=”blue”, lty=2) # övre (upper) gränsen
lines(x0$Ålder, mu0pred[,”lwr”], col=”red”, lty=3)
lines(x0$Ålder, mu0pred[,”upr”], col=”red”, lty=3)
e) Kontroll av förutsättningar. Vi skall nu kontrollera två av de antagandensom finns i analysen.
För det första antagandet om att slumpfelen skall vara normalfördelade och för det andra antagandet
om lika varians.
Residualerna beräkna med kommandot
res <- residuals(modell)
Undersök nu om residualerna är normalfördelade genom plotta residualerna i ett normalfördelnigspapper, använd R-funktionen normplot(res). Du kan behöva installera ett paket som heter nsRFA för
att få tillgång till funktionen normplot(). Ni kan också en göra en QQ-plot för att kompletera
undersökningen av residualernas fördelning (i R används kommandot qqnorm(res) ).
Antagandet om lika varianser (konstant spridning kring linjen) kan vi undersöka genom att plotta
residualerna mot Ålder. Gör det!
plot(torskar$Ålder, res)
Kan du se några strukturer i figuren?
Datamaterial: Bradfordmetoden.
I laborationen ”Proteinbestämning enligt Bradfordmetoden” i en kurs i cellbiologi undersöktes absorbansen hos prov med olika spädningar av Bovint Serum Albumin (BSA)-standard. Prov med 0-10 μg
protein spädes till 100 μl med vatten och två prover förberedes per koncentration.
Data för enav laborationsgrupperna finns i filen Labbdata.RData som du hittar på kursens hemsida.
Modell: Enligt Lambert-Beers lag gäller att absorbansen (A) kan beskrivas som en linjär funktion av
koncentrationen (c): A=k·c där konstanten k beror på ämnets molära absorptionskoefficient vid en viss
våglängd samt kyvettens längd. Vid mätningar får man naturligtvis räkna med en viss slumpmässig
variation, en rimlig modell är att absorbansen vid mätning nr i, A i, beskrivs linjärt av koncentrationen
ci plus ett slumpmässigt fel ei:
Ai = β0 + β1·ci + ei
där ei är oberoende och normalfördelad slumpfel med väntevärdet 0 och standardavvikelsen σ. Här
motsvaras konstanten β1 av den tidigare k medan β0 är absorbansen i den lösning som BSA:n är löst i.
a) Undersök på labbdata om den linjära regressionsmodellen ovan är rimlig att anpassa.
b) Hur mycket ökar absorbansen då man ökar koncentrationen en enhet? Ange ett 95 %
konfidensintervall för denna storhet
c) Vad är genomsnittlig absorbans för prov med koncentration 50 (mg/l). Ange ett 95 %
konfidensintervall för denna storhet. Tips: skapa en ny data.frame(…) och använd predict(…).
x50<-data.frame(koncentration=c(50))
mu50konf<-predict(modell2,x50,interval="confidence")
d) Vi har ett prov med koncentration 50 (mg/l). Ange ett 95 % prediktionsintervall för absorbansen
i just detta prov.
e) Huvudsyftet med mätningarna var att erhålla en standardkurva för hur absorbansen påverkas av
koncentrationen. Anta att vi på ett prov med okänd koncentration c 0 uppmätte absorbansen 0.43.
En skattning av c0 kan vi få fram genom att lösa ut x ur sambandet 0.43 = β0 + β1·x så här (om
den anpassade modellen sparats i variabeln modell2):
beta0 <- modell2$coefficients[1]
beta1 <- modell2$coefficients[2]
c0 <- (0.43 – beta0) / beta1
Vad blev den skattade koncentrationen?
f) Man skulle också vilja ha ett intervall som uppskattar inom vilka gränser c 0 kan ligga – ett sådant
intervall kallas kalibreringsintervall. Tyvärr kan man inte få det direkt i R. Däremot kan man få
en uppfattning om hur brett intervallet är eftersom kalibreringsintervallet är omvändningen till
prediktionsintervallet. Om vi beräknar prediktionsintervallet för ett antal koncentrationer kan vi
sedan undersöka för vilka koncentrationer som intervallet gränser blir den uppmätta absorbansen
0.43. Vi börjar med att beräkna prediktionsintervallet för koncentrationerna 0, 1, …, 100 och rita
upp dem samt den linje, y=0.43, som vi vill hitta skärningarna till:
x0 <- data.frame(koncentration=seq(0,100))
mu0pred <- predict(modell2, x0, interval=”prediction”)
plot(Labbdata$koncentration, Labbdata$absorbans, ylim=c(0.35, 0.5))
abline(modell2)
lines(x0$koncetration, mu0pred[,”lwr”], col=”red”, lty=3)
lines(x0$koncetration, mu0pred[,”upr”], col=”red”, lty=3)
abline(h=0.43) # horisontell linje på y=0.43
Vi kan försöka hitta skärningarna i figuren men eftersom vi har sparat beräkningarna av
prediktionsgränserna kan vi leta där istället. Vi skriver ut koncentrationerna och tillhörande
prediktionsintervall:
cbind(x0, mu0pred)
Leta upp det koncentrationsvärde där övre prediktionsgränsen upr är lika med 0.43 (vi kan kalla
det c_lwr eftersom det kommer att bli undre gränsen i kalibreringsintervallet). Leta också upp det
koncentrationsvärde där undre prediktionsgränsen lwr är lika med 0.43 (vi kallar denna koncentration
för c_upr).
c_lwr <- ? # skriv dit ditt avlästa värde istället för ?
c_upr <- ? # skriv dit ditt avlästa värde istället för ?
För att se hur det fungerade lägger vi till kalibreringsintervallet i figuren:
abline(v=c_lwr) # vertikal linje på x=c_lwr
abline(v=c_upr) # vertikal linje på x=c_upr
Vad blev kalibreringsintervallet? Är det användbart?
Datamaterial: Kids
Under ett givet år samlade man in vikten och längden hos nästan samtliga nyfödda vid Malmö
Allmänna Sjukhus, MAS. Materialet finns i datafilen kids.RData. Hämta datafilen från
hemsidan och undersök ifall det finns något beroende mellan variabeln kids$Langd och
variabeln kids$Vikt.
Tips:Börja med att plotta variablerna mot varandra. Finns det några anmärkningsvärda
datapunkter – så kallade outliers – ska dessa vara med vid modellbygget eller ej?
Hur stark är korrelationen mellan variablerna? Är den signifikant?
Om man bygger en linjär modell mellan de två variablerna. Ska man teckna sambandet som
att Vikten ska vara en linjär funktion av Längden eller ska man teckna sambandet som att
Längden är en linjär funktion av Vikten? Anpassa nu en linjär modell till datamaterialet:
Yi = β0 + β1·Xi + ei
Är parametrarna signifikanta? Vad blev förklaringsgraden R2? Har residualerna rätt fördelning
är de oberoende och med lika varians?
Är beroendet linjärt? Är en linjär modell bra? Väg ihop informationen om signifikansen av
parametrarna i modellen samt en enkel residualanalys. Tips: Kolla fördelningen hos
residualerna samt plotta dem både mot Vikt och Langd och leta efter strukturer?
Hur kan man göra modellen bättre? Kan tex sambandet vara kvadratiskt?
Om man vet att en nyfödd väger 3400 gram. Bestäm ett prediktionsintervall för barnets längd
med konfidensgraden 0.95.
Om man vet att längden hos en nyfödd är 48 cm. Bestäm ett prediktionsintervall för barnets
vikt med konfidensgraden 0.95.
Svar:
Torsk:
1. b) Ja! rp = 0.9828, t=15.0709 => p=0.000 – Man kan förkasta nollhypotesen om inget samband.
c) Förklaringsgrad R^2 = 0.966; r = 0.9828; s = 3.443; β0=5.993 (2.4044); β1=9.790 (0.6496);
t=15.071; p=3.72e-07.
d) Prognos för längd vid åldern 7 år = 74.52; PI, prediktionsintervall = [64.52, 84.53].
e) NF: Njä;Konstant varians: Ej helt lätt att bedöma (få värden)
Kids:
Det finns ett tydligt beroende mellan variablerna Vikt och Langd hos de nyfödda. Modelleringen kan
lika gärna göras som Vikt(Langd) som Langd(Vikt). Beroendet är inte linjärt. Leta efter strukturer i
residualplottarna. Plocka bort outliers och kontrollera ifall modellerna är bättre med eller utan outliers.
Bradfordmetoden:
2. b) 0.0008 intervall: 0.00063−0.0011
c) KI 0.425 – 0.441
d) PI 0.407 – 0.459
e) c0 skattas till 47 mg/l
f) Intervallet skattas till 15 – 78; tämligen brett.
Sammanfattning R
cor(x, y) Korrelationskoefficient
cor.test(x, y) Test för korrelationskoefficient
lm(y ~ x) regression av y som funktion av x
lm(y ~ x, data=dataframen som innehåller x och y)
… i ett visst datamaterial
summary(modell) skattningar, signifikanser, etc,
confint(modell) konfidensintervall för parametrarna
predict(modell, x0) prediction av förväntat värde när x=x0
predict(modell, x0, interval=”confidence”)
… med konfidensintervall
predict(modell, x0, interval=”prediction”)
… med prediktionsintervall
residuals(modell) residualer