Matematisk statistik allmän kurs, MASA01:B, HT-14 Laboration 2
Rapporten till den här laborationen skall lämnas in senast den 19e December 2014. Det går bra att mejla rapporten till mig.
Obs: Alla rapporter som jag får efter deadlinen kommer inte att godkännas.
Jag kommer rätta rapporterna innan jul och ni får titta på dem i Januari. Om ni har frågor gällande mina kommentarer så skicka gärna ett mejl så kan vi bestämma tid då vi träffas (därför att jag inte längre är anställd vid Lunds Universitet från och med början av Januari och pga detta inte längre kommer vara här varje vardag).
God jul till er alla och - om vi inte ses - lycka till med tentan!
Ifall ni har frågor eller något är oklart så är ni välkommna att höra av er.
Mareile Große Ruse, mareile@maths.lth.se, MH:325
1 ICKE PARAMETRISKA TEST
Laboration 4
I den här laborationen ägnar vi oss åt icke parametriska test och linjär regression.
1 Icke parametriska test
1.1 Chitvåtest och Fisher’s exakta test
I den första delen tittar vi på en kontingenstabell av följande form:
grupp/klass K1 K2
A OA,K1 OA,K2 nA B OB,K1 OB,K2 nB
n1 n2 n
Här är
nA= OA,K1 + OA,K2, n1= OA,K1+ OB,K1
(nB, n2 defineras på liknande sätt) och n = nA+ nB. Vi har lärt oss att statistikan
Q = (OA,K1 − eA,K1)2
eA,K1 +(OB,K1 − eB,K1)2
eB,K1 +(OA,K2− eA,K2)2
eA,K2 +(OB,K2 − eB,K2)2 eB,K2
har - under H0 - en approximativ χ2(1)-fördeling och vid testet av H0 använder man sig av χ2(1)-fördelningens p-värde.
Problem* 1.1.
• Vilken tumregel gäller för användbarhet av det approximativa chitvåtestet?
• Om tabellen ovan inte uppfyller tumregeln så kan man testa H0 med Fisher’s exakta test. På vilken (exakt) fördeling grundar sig Fisher’s exakta test?
I nästa uppgift undersöker vi kvaliteten av χ2-approximationen i det fall då vår 2 × 2 tabell inte uppfyller tumregeln.
I R kan en 2 × 2 tabell skapas på följande sätt:
grA <- c(111,876) # o_{A,K_1} = 111, o_{A,K_2} = 867 grB <- c(131,777) # o_{B,K_1} = 131, o_{B,K_2} = 777
# sammanfattar vektorerna till en tabell min.data <- as.table(rbind(grA, grB))
# ger namn åt rad och kolonn:
dimnames(min.data) <- list( grupp = c("A", "B"), klass = c("klass1", "klass2"))
# tittar på tabellen min.data
Fisher’s exakta test kan erhållas med hjälp av funktionen fisher.exact() och ett chitvåtest kan genomföras med chisq.test(). Vi får fram p-värdena genom fisher.exact()$p.value såsom chisq.test()$p.value. Om vi kör
fisher.test(min.data)$p.value # det exakta p-värdet
chisq.test(min.data)$p.value # det approximativa p-värdet så kan vi jämföra det exakta och det approximativa p-värdet.
Problem 1.2. Skaffa två olika 2×2 kontingenstabeller där du för vardera tabell väljer värdena så att tumbregeln inte är uppfylld. För varje tabell beräknas p-värdet av Fisher’s exakta test och p-värdet av chitvåtestet varefter värdena jämförs. Hur bra fungerar χ2-approximationen i sådana extrema fall?
1.2 Chitvåtest vid verkliga data
Den övriga delen av labben är baserad på verkliga data som finns i filen data.csv på kursens hemsida. Ladda ner filen, spara den på din dator och läs in den i R1. Kör
dat <- read.table("data.csv", header = TRUE, sep = ",")
Denna data är en del av data Pima Indian Diabetes som finns i UCI Machine Learning Repos- itory databas (http://www.ics.uci.edu/˜mlearn/MLRepository.html).
År 1990 testade den amerikanska National institute of Diabetes and Digestive and Kidney dis- eases 392 Pima indiska kvinnor för diabetes enligt Världshälsoorganisationens kriterier. När studien genomfördes var samtliga kvinnor äldre än 21 år och bosatta i närheten av Phoenix, Arizona.
Data består av n = 392 observationer av 11 variabler:
1
1 ICKE PARAMETRISKA TEST
name description scale
diastolic diastolic blood pressure integer-valued
glucose plasma glucose concentration in an oral integer-valued glucose tolerance test
insulin serum insulin concentration integer-valued
bmi Body Mass Index real-valued
diabetes diabetes pedigree function real-valued
age age in years integer-valued
test test whether the patient shows sign of diabetes positive, negative (coded "‘positive"’ if patient shows sign of diabetes)
dia.categ dichotomized diastolic blood pressure, high encodes low, high values > 70, low encodes values ≤ 70
age.categ dichotomized age ≤ 27, > 27
gluc.categ dichotomized glucose, high low, high
encodes values > 122, low encodes values ≤ 122
ins.categ dichotomized insulin, high low, high
encodes values > 156, low encodes values ≤ 156
Vi börjar med att skaffa oss ett överblick över data. Först tittar vi på de första raderna. Detta ger oss en uppfattning om variablernas värden.
dat[c(1:10), ] # tittar på de första 10 raderna I R finns många möjligheter för att plotta data.
Problem 1.3. Vi vill reda ut hur de stokastiska variablerna glucose och insulin hänger ihop.
Därför plottar vi glucose som en funktion av insulin (dvs vi plottar glucose på y-axeln och insulin på x-axeln):
plot(dat$insulin, dat$glucose) Finns något samband mellan variablerna?
En sådan plot kan vi får för alla par av slumpvariabler samtidigt (vi tittar här bara på de första 6 variablerna)
plot(dat[ , c(1:6)])
I föreläsningen har vi lärt oss om hur man använder ett χ2 oberoende test för att undersöka om två kontinuerliga slumpvariabler är oberoende eller inte.
I uppgiften ovan tittade vi på dem kontinuerliga slumpvariablerna glucose och insulin och plottade varje värdepar (insulini, glucosei), i = 1, . . . , n.
Vi kan också undersöka deras diskreta motsvarigheter ins.categ och glu.categ, som ger oss ett ytterligare sätt för att studera sambandet mellan båda variabler. Härför är hjälpsam funktion mosaicplot.
Funktionen mosaicplot kräver en värdestabell som input. En sådan kan skaffas genom glu.ins <- table(glucose = dat$gluc.categ, insulin = dat$ins.categ)
glu.ins
Problem 1.4. Gör en mosaicplot av glu.ins. Bekräftar mosaicplotten ditt intryck från förra uppgiften?
Testa (med ett lämpligt test) om sambandet mellan ins.categ och gluc.categ är signifikant.
Problem 1.5. Man kan tänka sig diabetes förekommer oftare bland äldre personer jämfört med yngre. Särskilt skulle detta betyder att ålder och diagnos av diabetes är beroende. Kan det här antagendet bekräftas statistiskt?
2 Enkel linjär regression
I den här andra delen av laborationen ägnar vi oss åt linjär regression.
Antagandet för en vanlig enkel linjär regression är att responsvariabeln Y har formen Y = β0+ β1x + ,
där x är en (känd) förklarande variabel och ∼ N (0, σ2). Parametrarna β0, β1, σ2 > 0 är okända och skall skattas.
Dessutom antar man att man har ett stickprov av n oberoende värden Y1, . . . , Yn av Y , Yi = β0+ β1xi+ i,
där 1, . . . , n är oberoende och har samma fördelning som . x1, . . . , xn är (kända) värden av den förklarande variabeln x, t.ex. är Yi den uppmätta längden av holländare i, xi är skostorleken av holländare i och i är mätfelet vid längdens mätning.
Problem* 2.1. Vad har Yi, i = 1, . . . , n för fördelning?
2 ENKEL LINJÄR REGRESSION
Problem* 2.2. Ange ML-skattaren ˆβ0, ˆβ1 av β0, β1.
Problem 2.3. Skatta täthetsfunktionen av variabeln glucose och lägg till täthetsfunktionen av en normalfördelning med lämpliga värden på väntevärde och varians.
Verkar glucose vara normalfördelad?
Ovan har vi sett att glucosvärdet och insulinvärdet har ett starkt samband. Nu vill vi göra en linjär regression med glucose som responsvariabel och insulin som förklarande variabel.
Dvs vi modellerar
glucosei = β0+ β1 insulini+ i, i = 1, . . . , n (1) där i∼ N (0, σ2).
Problem 2.4. Beräkna ML-skattarna ˆβ0 och ˆβ1.
I R finns en funktion, lm(), som räknar ut skattarna automatiskt. Det första argumentet till lm() är responsvariabeln, därefter följer tecknet ˜ som skiljer av ekvations (1) vänstra sida från den högra sidan. Sist anger man namnet insulin av den förklarande variabeln.
linreg <- lm(glucose ~ insulin, data = dat)
Jämför dina skattade värden (uppgift 2.4) med de skattade värdena som ges av lm funktionen linreg. Ifall värdena inte stämmer överens gör om uppgift 2.4.
Ofta vill man testa om det skattade värdet ˆβ1 skiljer sig signifikant från 0. Ett hypotespar H0: β1 = 0 vs H1 : β1 6= 0
är ekvivalent med påståenderna
H0 : Det finns inget linjärt samband mellan glucose och insulin vs H1 : Det finns ett linjärt samband mellan glucose och insulin
För att kunna testa H0 mot H1 behöver vi fördelningen av skattaren ˆβ1 under H0: β1= 0.
Man kan visa att - under H0 -
βˆ1 q
V(βˆ1)
∼ N (0, 1).
Om σ är okänd så är också V( ˆβ1) okänd (variansen av ˆβ1 beror på σ) och behöver skattas. I så fall är
T (Y ) = βˆ1
q V(\βˆ1)
∼ t(n − 2).
Därför kan vi testa H0 mot H1 med ett t-test. Vi erhåller den skattade standardavvikelsen q
V(\βˆ1) genom
summary(linreg)$coefficients
Den skattade standardavvikelsen av ˆβ1 finns i rad två , andra kolonnen. t-statistikans värde står i den tredje kolonnen och p-värdet i den sista kolonnen.
För att få en bättre förståelse gör vi följande beräkningar:
beta1_hat <- summary(linreg)$coefficients[2,1]
sd_beta1_hat <- summary(linreg)$coefficients[2,2]
beta1_hat sd_beta1_hat
teststat <- beta1_hat/sd_beta1_hat teststat
n <- length(y)
pvärde <- 1-pt(teststat, n-2) pvärde
och ser att våra värden stämmer överens med värdena vi fick av summary(linreg)$coefficients.
Vi får även ett 95%-konfidensintervall för β1 (och för β0, men vi är bara intresserade av β1, dvs av den andra raden):
confint(linreg)
eller om vi beräknar på egen hand:
c(beta1_hat - qt(0.975, n-2)*sd_beta1_hat, beta1_hat + qt(0.975, n-2)*sd_beta1_hat)
